Research ArticleBIOPHYSICS

Confinement-induced self-organization in growing bacterial colonies

See allHide authors and affiliations

Science Advances  22 Jan 2021:
Vol. 7, no. 4, eabc8685
DOI: 10.1126/sciadv.abc8685


We investigate the emergence of global alignment in colonies of dividing rod-shaped cells under confinement. Using molecular dynamics simulations and continuous modeling, we demonstrate that geometrical anisotropies in the confining environment give rise to an imbalance in the normal stresses, which, in turn, drives a collective rearrangement of the cells. This behavior crucially relies on the colony’s solid-like mechanical response at short time scales and can be recovered within the framework of active hydrodynamics upon modeling bacterial colonies as growing viscoelastic gels characterized by Maxwell-like stress relaxation.


The ability to exploit the physical properties of the environment, to achieve biological functionality at the large scale, is the hallmark of self-organization in multicellular systems (14). Both in eukaryotes (1, 2) and prokaryotes (3, 4), cells can take advantage of the geometrical and mechanochemical features of their surroundings, such as confinement, friction, compliance, and degradability, to migrate and proliferate. In prokaryotes, a spectacular example of this behavior is offered by growing colonies of sessile bacteria. Although lacking of a self-propulsion machinery, these cells are still capable of taking advantage of the extensile forces resulting from their individual growth to achieve collective migration while relying on the environmental topography to navigate (5).

Volfson et al. (6) investigated this mechanism using a monolayer of dividing nonmotile Escherichia coli cells cultured in microfluidic channel with open ends that allow the cells to escape. Such a confined bacterial layer is initially structureless, but, as the density progressively increases, global nematic order starts to develop within the system, with most of the cells oriented along the longitudinal direction of the channel. This alignment mechanism, resulting from the combination of growth and confinement, allows the cells to efficiently relieve the internal extensile stresses and prevents the colony from overcrowding by facilitating the emergence of an expansion flow directed toward the channel outlets. In long channels, however, such a globally ordered state is unstable, and domains of nonlongitudinally oriented cells arise throughout the colony, rendering it disordered at the large scale (7).

Whereas the latter instability, triggering the transition from global alignment to disorder, is well understood, the nature of the mechanism leading to global alignment is still debated, despite markedly affecting the overall fitness of the colony. Volfson et al. (6) proposed that global longitudinal alignment arises in response to the bacterial expansion flow, because of the propensity of nematic liquid crystals to align with respect to the flow direction. Thus, the channel geometry dictates the direction of the expansion flow, which, in turn, acts as an ordering field for the bacteria themselves. A different mechanism, recently proposed by Karamched et al. (8), revolves instead around the role of parallel anchoring at the channel walls. Using a stochastic model based on the spatial Moran process, they demonstrate that global alignment can emerge even in the absence of flow alignment. The latter finding is, however, built upon specifically designed growth dynamics, which, based on previous experiments and simulations (6, 7, 911), are by no means indispensable for the occurrence of global alignment in vitro.

Here, we propose an alternative explanation, rooted in the emergence of globally anisotropic stresses in the confined bacterial population. Using computer simulations of a hard-rod model, we find that cell growth gives rise to a persistent accumulation of mechanical stress in the colony. While longitudinal stresses can be efficiently relieved via the expansion flow toward the channel outlets, lateral confinement prohibits cell motion in the transverse direction, leading to a buildup of transversal stress. The resulting mechanical anisotropy determines a net torque that reorients the cells in the direction of minimal stress, which, in the case of a straight channel, coincides with the longitudinal direction. This interpretation is additionally supported by a continuum model of the expanding colony as an active viscoelastic nematic gel, characterized by Maxwell-like stress relaxation. These discoveries not only deepen our understanding on the roles of confinement and mechanical stresses on the self-organization of growing bacterial colonies, providing a potential way to control growing bacterial colonies, but also shed new light on the interplay between orientation and stress in active viscoelastic liquid crystals.


Growing colonies under confinement

We use a minimal model of proliferating hard rods to study growing bacteria under confinement (1214). Each cell is modeled as a spherocylinder, which can grow, divide, and push away its neighbors as it elongates. Details of the model can be found in the “Hard-rod model” section. Figure 1 (A to R) illustrates the typical dynamics of our in silico colonies, subject to boundary conditions of three different types. When the colony is confined by rigid walls (Fig. 1, A to F), our results reproduce previous experiments and simulations (611, 15, 16): The long time dynamics consists of a steady state characterized by global nematic order, with most of the cells aligned along the horizontal direction. A possible explanation of this phenomenon, proposed in (8), relies on the intuitive idea that steric repulsion from the rigid walls generates torques on the peripheral cells, thereby aligning these horizontally. Such an alignment then propagates to the interior of the channel as a consequence of the cell-cell steric repulsion, resulting into uniform horizontal alignment across the entire colony. The latter mechanism occurs in thermotropic nematic liquid crystals [see, e.g., (17)] and was recently observed in confined monolayers of eukaryotic spindle-like cells (18). However, our results on periodically confined colony (Fig. 1, G to L) readily disprove this hypothesis. In this instance, the colony is not confined by rigid horizontal walls; however, global horizontal alignment is still prominent. Specifically, once the top and bottom boundaries of the colony merge, disconnected domains of horizontally aligned cells start to emerge uniformly throughout the colony (Fig. 1J). Over time, these domains expand and eventually form a percolating cluster spanning the entire length of the colony (Fig. 1, K and L). Global alignment emerges even when vertical confinement is removed (Fig. 1, M to R). Here, a region of horizontally aligned cells appears in the center of the channel once the colony height exceeds a given threshold and then gradually expands vertically, leaving two disordered caps of constant thickness at the top and bottom of the colony. Last, increasing the channel width eventually disrupts global horizontal alignment (see fig. S1). In particular, the wider the channel, the larger the height of the colony at which the horizontally aligned region first appears. In addition, nematic order in wide channels is lower on average and the disordered caps are thicker than those in narrow channels. The latter phenomenon is well suited to be tested with in vitro experiments, for instance, by growing multiple bacterial colonies within confining channels of varying width.

Fig. 1 Colonization under confinement.

Snapshots of growing colonies at different time points. The boundary conditions are (A to F) rigid-wall confinement, (G to L) periodic confinement, and (M to R) no confinement in the y direction. Cells are color-coded by their orientation, as indicated by the color wheel in (A). In all three colonies, Lx = 70 μm, while Ly = 50 μm in (A) to (L). The colony height in (R) is about 300 μm. In all three scenarios, the direction of the confining channel is orthogonal to the red outlets, hence horizontal.

Anisotropic stress drives cell reorientation

Volfson et al. (6) proposed that the longitudinal alignment of E. coli bacteria in a channel resulted from the combination of the expansion flow, caused by cell division, and the tendency of nematic liquid crystals to reorient in the presence of a velocity gradient. Whereas plausible to explain the observed longitudinal alignment within a straight channel, however, this mechanism would also determine preferential radial alignment in freely expanding colonies, as a consequence of prominent radial flow characterizing these systems [see, e.g., (14)]. By contrast, such a behavior has never been observed in either experiments or simulations (19). Furthermore, longitudinal alignment can also emerge in the absence of flow, provided that the cells are subject to anisotropic normal stresses.

To demonstrate this latter statement, we have simulated an initially disordered distribution of nongrowing cells, confined in a rectangular doubly periodic box and subject to a uniform contraction along the y direction (Fig. 2). Initially, each cell has a random position and orientation (Fig. 2A). Then, the box height is uniformly shortened by means of the scaling transformation: y(t + Δt) = ky(t), where y is the abscissa of the cells and the horizontal boundaries, k = (LyVyΔt)/Ly is the shrinking factor, and Vy is the relative speed of the two horizontal boundaries. Note that cell orientations remain unchanged during the scaling. During the process, we measure the stresses experienced by the cells using the virial construction, namelyσi=1aiΣjrij Fijc(1)where ai=ai/ϕ is the effective area occupied by the ith cell, with ai = d0(li + πd0/4) the exact cell area and ϕ the local packing fraction, and rij is the position of point of contact between the ith and the jth cell with respect to the center of mass of the ith cell. To characterize the amount of global orientation, we define the following alignment parameterΦ=2px21(2)where 〈 · 〉 denotes an average over all cells in the box. By construction, −1 ≤ Φ ≤ 1, with Φ = 0 corresponding to a completely disordered configuration and Φ = ± 1 to perfect alignment in x and y, respectively. Figure 2 (A to C) summarizes the dynamics of the shrinking colony. The system is initially isotropic and the orientation of the cells is uniformly distributed (Fig. 2A and its inset). Shrinking increases the cells’ overlap, hence the stress across the colony. However, because of the anisotropy in the scaling transformation, the vertical normal stress σyy increases more quickly than the horizontal normal stress σxx (Fig. 2D). Simultaneously, longitudinal alignment develops throughout the colony (Fig. 2, C and D).

Fig. 2 Dynamics of the shrinking colony.

(A to C) Snapshots at different time points. Cells are color-coded by their orientations according to the color wheel in (C). The insets show the histograms of cell orientation at the corresponding time points. (D) Average stresses and the alignment parameter Φ as functions of time. The system contains 5000 cells in total, with a width Lx = 150 μm and an initial height Ly = 200 μm. The shrinking speed, i.e., the relative speed between the top and the bottom boundaries, is Vy = 10 μm/hour.

These observations suggest a causal relation between the occurrence of transient anisotropy in the normal stresses and the emergence of longitudinal alignment within the bacterial population. Here, we postulate the following mechanism. Passive spherocylinders in proximity of the isotropic-nematic phase transition are known to organize in clusters consisting of highly aligned cells, sometimes referred to as cybotactic clusters (20, 21). This effect is further enhanced in colonies of growing bacteria, as a consequence of the long-wavelength instability of the nematic ground state driven by the extensile active stresses (14). These clusters are not held together by attractive interactions and eventually break up or merge into larger domains. At short time scales, however, clusters may exhibit solid-like mechanical response to environmental forces and, in particular, undergo internal rearrangements while subject to anisotropic stresses, resulting in a redistribution of the stress that eases the local stress anisotropy. Figure 3 illustrates this process in the case of uniformly shrunk colonies. The cells in the central part of Fig. 3A are initially loosely packed and do not exhibit a preferential orientation, but, as time progresses, the shrinkage causes them to form a stack of nine tightly packed cells roughly oriented at 45 with respect to the horizontal direction (Fig. 3B). Last, as σyy > σxx, this cluster undergoes a collective rearrangement, whose effect is to redistribute the stress among the normal components by reorienting the cells along the horizontal direction (Fig. 3C).

Fig. 3 Alignment mechanism in uniformly shrunk colonies.

(A) The cells in the central part of the panel are initially loosely packed and do not exhibit a preferential orientation. (B) As time progresses, the shrinkage causes them to form a stack of nine tightly packed cells roughly oriented at 45° with respect to the horizontal direction. (C) Last, as σyy > σxx, this cluster undergoes a collective rearrangement, whose effect is to redistribute the stress among the normal components by reorienting the cells along the horizontal direction. Snapshots correspond to different times in the simulations of a shrunk colony.

Some remarks are in order. First, the aligning mechanism illustrated above in the case of uniformly shrunk colonies results from the orchestrated action of multiple effects. These include the entropic alignment of the cells originating from their excluded volume interactions, further enhanced by the progressive increase in density, and the cluster’s collective rotation in the presence of anisotropic normal stresses. The latter process, in turn, crucially relies on the fact that cells can slide along their longitudinal direction p, while their motion in the transversal direction p is obstructed. Second, because of the absence of attractive interactions among our in silico cells, the same effect cannot be produced solely by shear stresses, as these would mainly slide the cells with respect to each other, resulting in a spreading of the domain along the horizontal direction. Real bacteria do exhibit attractive interactions, mediated by proteins such as the adhesin FimH (22), but these are inevitably weaker than the repulsive excluded volume interactions among cells. Last, such an alignment mechanism evidently relies on the viscoelastic nature of bacterial layers, and its performance depends on how the time required for the stress anisotropic to build up compares with the viscoelastic crossover time. In the “Continuum description” section, we directly test this mechanism by implementing it in a continuum viscoelastic model of growing colonies characterized by Maxwell-like stress relaxation.

Stress anisotropy and alignment in growing colonies

To test whether the mechanism proposed in the previous section carries over to colonies of growing bacteria, we look back at the expanding colonies depicted in Fig. 1 and track the alignment parameter Φ as well as the normal stresses σxx and σyy during expansion. Generally speaking, the normal stress decreases with the distance from the colony center and increases monotonically in time. To more precisely capture such spatial-temporal dynamics, let us focus on a region ℛ0, defined as the area within 20 μm from the colony center: i.e., x2+y2<20 μm. Figure 4 (A to C) shows the average normal stresses and the alignment parameter within ℛ0 of the three colonies displayed in Fig. 1, subject to different boundary conditions. As the fronts of the colonies approach the boundaries, stress is globally isotropic and the colonies have the typical structure observed in free space (14). This behavior drastically changes once the colonies meet the boundaries at t = tc. In particular, in colonies subject to hard-wall (Fig. 1, A to F) or periodic confinement (Fig. 1, G to L), both normal stresses undergo a marked increase; hence, σyy > σxx as a consequence of the inhibition of the bacterial motion in the y direction (Fig. 4, A and B). A similar behavior is also found in unconfined colonies (Fig. 1, M to R). Here, stress builds up more gently than in the presence of hard walls, and yet, σyy > σxx at any time t > tc (Fig. 4C). This results from the fact that, while the removal of cells in the x direction relieves the associated normal stress σxx, σyy increases with time as the colony elongates along the y direction. Intuitively, such an implicit confinement is much weaker than the explicit confinement provided by a hard wall or a periodic boundary; thus, stresses increase more gradually. In all three colonies, the occurrence of stress anisotropy, represented by the nonvanishing normal stress difference ∣σyy − σxx∣, continuously drives cells to align with the horizontal direction. Eventually, a steady state is reached where most of the cells are aligned along the x direction; thus, Φ ≳ 0.5 and σyy ≳ σxx (see, e.g., Fig. 1, F, L, and R).

Fig. 4 Stress dynamics in growing colonies.

Normal stresses and the alignment parameter in region ℛ0 of the colonies. The boundary conditions are (A) hard-wall confinement (see, e.g., Fig. 1, A to F), (B) periodic confinement (Fig. 1, G to L), and (C) no confinement (Fig. 1, M to R). In (A) and (B), the black dashed line indicates tc, the time at which the colonies reach the horizontal boundary. In (C), tc is set as the time where the normal stress difference in region ℛ0 becomes nonzero. The results at each boundary condition are obtained by averaging over 200 simulations of growing colonies.

Despite the aligning mechanisms having the same origin in both shrinking and growing colonies, the latter exhibits a number of system-specific properties that are not found in passive systems. First, shrinking colonies undergo a jamming transition for sufficiently high densities; this effectively freezes the cells’ orientation. By contrast, in growing colonies, the perpetual cell growth and division prevent the system from jamming, thus allowing the colony to further improve its global longitudinal alignment. This is exemplified by the fact that Φ is systematically larger in growing colonies than in shrinking ones. Second, because of the aligning mechanism described in the previous section, domains of horizontally aligned cells are more stable than others. As cells duplicate, these horizontal domains can then expand without being reoriented and eventually take over the entire colony. Last, as we detail in the “Continuum description” section, cell growth results in a continuous injection of active stress σa ∼ 〈pp〉. This gives rise to a feedback mechanism that directly converts orientational anisotropy into stress anisotropy, in such a way as to enhance the redistribution of stress from the transverse to the longitudinal component. This effect stabilizes the horizontally aligned configuration, where the normal stress difference ∣σyy − σxx∣ is minimal at a given cell density.

To further check our hypothesis, we explore the onset of stress anisotropy and longitudinal alignment in the unconfined colony displayed in Fig. 1 (O to R). For this purpose, we focus on a 20-μm-width region located at the center of the colony and spanning its entire height. Figure 5 shows the spatial distribution of the stresses and the alignment parameter Φ at different times, or, equivalently, colony heights. At early times and for small colony heights, stress is isotropic and the cells have random orientation (Fig. 5A). Once the colony becomes sufficiently elongated in the y direction, however, the normal stresses in the center of the colony start to differentiate as a consequence of the effective confinement originating from the top and bottom caps (Fig. 5B). The anisotropy increases with time as the colony progressively elongates (Fig. 5, C and D). The alignment parameter Φ and the normal stress difference ∣σyy − σxx∣ transition from vanishing to positive simultaneously and continue evolving in parallel at any later time, thus confirming the correlation between these two quantities.

Fig. 5 Stress distribution in growing colonies.

Spatial distribution of stresses and alignment parameter Φ in a 20-μm-width region located at the center of the colony and spanning its entire height at different times. All the quantities are measured upon dividing the region into smaller boxes of height 10 μm and calculating their averages within each box. The channel is Lx = 70 μm wide, whereas the colony’s height is (A) 70 μm, (B) 150 μm, (C) 250 μm, and (D) 390 μm. The results are obtained by averaging over 100 simulations of growing colonies.

To quantify such a correlation, we construct a two-dimensional map of the alignment parameter Φ as a function of σxx and σyy for the case of unconfined colonies (Fig. 6A). This is achieved by varying the width of the channel in the interval 50 μm ≤ Lx ≤ 150 μm in steps of 10 μm and performing 100 simulations for each Lx value, with the goal of generating a large amount of (σxx, σyy, Φ) combinations. As it is evident from Fig. 6A, the full dataset only occupies a narrow wedge-shaped region of the (σxx, σyy) plane, straddling the ∣σxx∣ = ∣σyy∣ line, as a consequence of the stress redistribution resulting from the rotation of the cells. Such a region is slightly asymmetric toward the ∣σyy∣ > ∣σxx∣ half-plane because of the elongated morphology of the colony and the continual release of σxx by the absorbing boundaries. This corresponds to a region of increased horizontal alignment, marked in red in Fig. 6A, which is more likely to occur in the portions of the region characterized by large stress anisotropy. The latter property is even more evident in Fig. 6B, where the alignment parameter Φ is plotted directly against the stress anisotropy parameter, defined asΔΣ=σyyσxxσxx+σyy(3)thus confirming our hypothesis that bacterial alignment is ultimately related with the occurrence of a confinement-induced anisotropy in the normal stresses.

Fig. 6 Stress anisotropy in growing colonies.

Relation between the stresses and the alignment parameter in colonies unconfined in the y direction. (A) Two-dimensional map of the anisotropy parameter Φ as a function of the normal stresses σxx and σyy in channels of increasing width, in the range Lx = 50 to 150 μm with ΔLx = 10 μm. The cyan line corresponds to the bisectrix ∣σxx∣ = ∣σyy∣. The maxima of Φ correspond to the maxima of ΔΣ for a given Lx value, and their periodicity results solely from the combination of multiple Lx values, each associated with a maximum in Φ and ΔΣ. (B) Stress anisotropy parameter, Eq. 3, versus the alignment parameter. The error bars show the SDs of data samples about the average values.

Last, topological defects are observed in the early stages of the system for both discrete and continuous models. However, their presence is transient and does not affect the long-time behavior of the colony.

Continuum description

As a final test of the alignment mechanism proposed in the previous sections, we introduce a continuum model where most of the emergent properties found in our hard-rod model can be directly implemented and controlled. For this purpose, we describe monolayers of growing bacteria as two-dimensional viscoelastic active nematic gels, characterized by three continuum fields: the cell density ρ, the local mechanical stress σ, and the nematic tensor Q indicating the local cell orientation. More details of the continuum model can be found in the “Continuum model” section.

Figure 7 (B to D) summarizes the main outcome of the model when κ > 0. After an initial transient, our continuum colony relaxes on a steady state characterized by perfect longitudinal alignment throughout the channel (Fig. 7B, where Qxy = S/2 sin 2θ = 0 and Qxx = S/2 cos 2θ > 0) and a stationary flow in the x direction (Fig. 7C). By contrast, the density ρ − ρ0 is maximal at x = 0 and vanishes at the outlets, as a consequence of the accumulation of cells in the center of the channel. A comparison between Figs. 7D and 6B, in particular, corroborates our hypothesis that the faster buildup of transverse stress, originating from the inhibition of the flow along the y direction, drives the longitudinal alignment of the nematic director. In Fig. 7D, different combinations of (σxx, σyy, Φ) are sampled during the relaxation dynamics of the colony, and the error bars indicate the SD.

Fig. 7 Continuum model.

(A) Schematic representation of a one-dimensional active Maxwell material. Activity is incorporated into the constitutive equation in the form of a generator supplying a constant stress σa. (B) Steady-state configurations of nematic tensor Q for κ > 0. After an initial transient, the colony relaxes toward a state characterized by perfect longitudinal alignment. (C) Steady-state configuration of the density (black), vx (blue), and vy (red) in a colony with κ > 0. (D) Stress anisotropy parameter ΔΣ, Eq. 3, versus the alignment parameter Φ for a colony with κ > 0, compare with Fig. 6B. (E) Steady-state configurations of Qxx (or Φ/2), Qyy, and the stress anisotropy versus κ. As κ changes in sign from positive to negative, the colony switches from longitudinal to transverse orientation. The results are obtained from a numerical integration of Eqs. 8A to 8C on a rectangular domain with adsorbing boundaries in the x direction and periodic boundaries in the y direction.

Figure 7E illustrates the effect of the coupling parameter κ on the performance of global alignment. For each κ value, we run 50 simulations with different random initial configuration of the Q field to generate statistics. For κ > 0, the colony rapidly evolves toward a longitudinally aligned configuration. The latter, in turn, is insensitive to the specific κ value as exemplified by the plateau in Φ shown in Fig. 7E. By contrast, for κ < 0, the colony exhibits global vertical alignment (i.e., Φ ≈ −1), but the performance of the alignment mechanism is reduced compared to the κ > 0 case, as demonstrated by the large error bars. We expect this behavior to originate from the incompatibility between vertical alignment and the release of active stress into the bacterial flow. While in the case of horizontally aligned cells the active stress σaαx̂x̂ originating from bacterial growth can be instantaneously released in the x direction by sourcing a longitudinal expansion flow sinking at the absorbing boundaries, the active stress σaαŷŷ generated by vertically aligned cells must first be redistributed to the other stress component before being released through the flow.

Here, we have investigated the interplay among orientational order, geometrical confinement, and growth in multicellular systems of sessile bacteria. Using molecular dynamics and continuous modeling, we have demonstrated that geometrical anisotropies in the confining environment give rise to an imbalance in the normal stresses, which, in turn, drives a collective reorientation of the cells. Cell reorientation not only eases the mechanical anisotropy by redistributing the existing stress among the two normal components, but also redirects the growth-induced flow in such a way to facilitate the escape of cells and reduce further stress buildup. The combination of these mechanisms results in a three-way regulatory loop that, although having a purely mechanical origin, allows bacterial colonies under confinement to control internal stresses and navigate through their environment in such a way to maximize their fitness. As an example of this behavior, here, we discussed the case of bacterial colonies in a channel, where the fast buildup of transverse stress drives the cells to reorient along the longitudinal direction and proliferate throughout the channel. A chief prediction of our analysis is the relation between the stress anisotropy ΔΣ and alignment parameter Φ (see Figs. 6B and 7D), which could be, in principle, experimentally tested upon confining bacterial colonies with embedded stress sensors in a microchannel.

Aside from the microfluidic setting, our work sheds light on special circumstances where confinement is induced by the colony itself, even in the absence of physical barriers. An interesting illustration of such a self-confinement effect is provided by the merger of bacterial communities as they independently colonize the surrounding environment (Fig. 8). In the merging region (enclosed by colored boxes), cells collectively align tangentially with respect to the colonies’ front. This observation can be interpreted on the basis of our results by noticing that the cells at the interface between the merging colonies are confined by their respective bulks, thus subject to a normal stress σNN = N · σ · N that eventually exceeds the tangential stress σTT = T · σ · T, where N and T are, respectively, the tangent and normal directions at the colony’s front. Like colonies growing within an open channel, these interfacial cells are then reoriented in the tangential direction. This reorganization favors the merging of the bacterial communities and allows them to cooperatively colonize the surrounding space, rather than compete for the available resources.

Fig. 8 Bacteria alignment in merging colonies.

Experimental snapshots of the merger of three Vibrio cholerae colonies at different times. In the merging regions, marked with colored boxes, cells collectively align with the common tangent at the colonies’ front. Courtesy of A. Sengupta.

Our work additionally suggests a possible strategy to control the expansion trajectory of growing bacterial colonies with engineered anisotropic stresses. This could be achieved, for instance, by designing apposite confinement geometries or by introducing cell sinks at specific locations, in such a way to adjust the direction of minimal stress. These sinks could be, in principle, opened and closed on demand to achieve a fully controllable collective bacterial flow.

Last, our analysis advances current knowledge about hydrodynamic stability of active fluids, by highlighting the role of confinement and viscoelasticity as stabilizing features. It has been extensively shown, both theoretically and experimentally, that confined active fluids depart from the lowest free-energy state once the system size L exceeds the characteristic length scale 𝓁a at which active and passive forces and torques balance (2328). For L ≫ 𝓁a, the system becomes chaotic (2933), often after having experienced an intermediate regime characterized by the appearance of coherent structures and periodic flow patterns (3437). Here, we showed that the combination of stress anisotropy and viscoelasticity enhances the stability of the lowest energy state, even in the presence of strong active flows, thus suggesting a new angle to explore the role of confinement in active matter.


Hard-rod model

Cells are modeled as spherocylinders of a fixed diameter d0 constrained to lie in the xy plane. Each cell is characterized by three time-dependent degrees of freedom—namely, the position ri; the orientation pi=cos θi x̂+sin θi ŷ, with θi the angle with respect to the x axis of a Cartesian frame; and the length li (excluding the caps on both ends)—whose dynamics is governed by the following set of overdamped ordinary differential equationsdridt=1ζli (Σj=1NiFijc+Σα=±1Fiαw+ηi)(4A)dpidt=12ζli3 Mi×pi(4B)

The quantity ζ represents the effective drag per unit length originating from the cell-substrate adhesive interactions [see, e.g., (38)]. Adhesive molecules expressed on the bacterial surface dynamically bind and unbind to the substrate and can be considered uniformly distributed along the cellular wall on average. The resulting drag force is, consequently, proportional to the cell area and is independent of the cell orientation. The first term on the right-hand side of Eq. 4A represents the steric force due to the contact between the ith and jth cell, with 𝒩i the total number of cells in contact with the ith. The points of contact have positions rij with respect to the center of mass of the ith cell and apply a Hertzian forceFijc=Yd01/2hij3/2Nij(5)where Y is proportional to the Young’s modulus of the cell, hij is the overlap distance between the two cells, and Nij is their common normal. The second term on the right-hand side of Eq. 4A, on the other hand, represents the force associated with the interaction between the cell caps, at positions riα = αlipi/2 with respect to the center of mass, and the confining walls and is calculated analogously to Fijc. The vector ηi corresponds to a random force, whose components are sampled from the uniform distribution in the range −10−6 N ≤ ηiβ ≤ 10−6 N, with β = x, y. This force embodies the inevitable randomization of the bacterial dynamics resulting from small irregularities in the cell walls and the substrate, preventing the dividing cells from forming long one-dimensional chains that are never observed in experiments. Analogously, Eq. 4B describes the rotation of the cell axis in response to the torqueMi=Σj=1Ni(rij×Fijc)+Σα=±1(riα×Fiαw)(6)

Last, the length li increases linearly in time with rate gi and, after having reached the division length ld, the cell divides into two identical daughter cells. To avoid synchronous divisions, the growth rate is randomly sampled from the uniform distribution in the range g/2 ≤ gi ≤ 3g/2, with g the average growth rate. Immediately after duplication, the daughter cells have the same orientation as the mother cell, but independent growth rates. Equation 4 has been numerically integrated using the following set of parameter values, corresponding to the typical values of E. coli microcolonies under physiological conditions [see, e.g., (12)]: d0 = 1 μm, ld = 2 μm, g = 2 μm/h, Y = 4 MPa, and ζ = 200 Pa h. The integration is performed via the Euler scheme with a time step Δt = 0.5 × 10−6 hours.

The proliferating colony is confined in a rectangular Lx × Ly box, whose horizontal (top and bottom) boundaries are either consisting of impenetrable rigid walls, periodic, or placed at an infinite distance from the colony. The vertical (left and right) boundaries, on the other hand, are always set as absorbing, in such a way to emulate the outlets of a microfluidic channel (6, 7, 911). Cells crossing the absorbing boundaries are then removed from the system. For the sake of brevity, in the remainder of the paper, we refer to Lx as the “width” and to Ly as the “height” of the channel. All simulations start with a cell at the origin, with random orientation. We change only the boundary condition and the channel geometry and fix all the other parameters.

Continuum model

We describe monolayers of growing bacteria as two-dimensional viscoelastic active nematic gels, characterized by a Maxwell-like stress relaxation (Fig. 7A), namely2ηu=(1+τDt)(σσa)(7)where uij = (∂ivj + jvi)/2 is the strain rate tensor and the constants η and τ represent, respectively, the shear viscosity and the viscoelastic crossover time: i.e., τ = η/E, with E the Young’s modulus of the colony. The operator Dt is the corotational material derivative: i.e., Dtσ = tσ + v · ∇ σ + ω · σσ · ω, with ωij = (ivjjvi)/2 the vorticity tensor. The tensor σa = αQ represents the extensile active stress resulting from the growth of the cells (23, 39), whereas Q=pp1/2=S(nn1/2) is the nematic tensor, with 𝟙 the two-dimensional identity tensor, n the nematic director, representing the average orientation of the cells, and S = 〈2 ∣n · p∣ −1〉 the local order parameter. At long times, i.e., t ≫ τ, Eq. 7 yields the usual constitutive equation of active fluids: i.e., σ = 2ηu + σa; while at short times, Eq. 7 describes the characteristic response of a solid material: i.e., tσu.

The hydrodynamic equations for the cell density ρ, momentum density ρv, and nematic tensor Q are given bytρ+·(ρv)=kgρ(8A)t(ρv)+·(ρvv)=·σ˜Pξρv(8B)DtQ=γ1H˜(8C)where the tilde indicates the deviatoric (i.e., traceless) component: e.g., σ=σtr (σ)1/2.

Equation 8A accounts for the exponential growth in the number of bacteria resulting from cell division, with kg the growth rate: i.e., N = ∫ dA ρ/m = N0 exp kgt, with m the cell mass. On the other hand, Eq. 8B describes the redistribution of momentum among the stresses, with P the pressure, as well as its dissipation via the drag force −ξρv, with ξ a drag coefficient, resulting from the interaction with the substrate. Last, Eq. 8C governs the evolution of the bacterial orientational degrees of freedom, with γ the orientational viscosity and H = − δFQ the molecular tensor describing the relaxation of the free energy F = ∫ dA f. The latter must account for both the entropic contribution, resulting from a departure from the uniformly aligned configuration, and the elastic contribution associated with the colony solid-like behavior and can be expressed in the classic de Gennes’ form (40)f=12KQ2+12Atr Q2+14C(tr Q2)2Γtr (Q·ϵ)(9)with K the orientational stiffness of the nematic phase, in one elastic constant approximation, and A and C being mean-field coefficients controlling the isotropic-nematic phase transition. As bacterial colonies are generally not subject to thermal fluctuations and nematic order arises as a consequence of crowding and activity [see, e.g., (19)], one can assume the order parameter S=2tr Q2 to be nonzero even in two dimensions. The constant Γ expresses the strength of the coupling between the nematic tensor and the strain tensor ϵ. Consistently with Eq. 7, we assume this to be related to the stress tensor via Hooke’s law: i.e., ϵ=(1+ν)/E σν/E tr (σ)1, with ν the Poisson ratio. ThusH˜=K2Q+(A+12CS2)Qκσ˜(10)where κ = Γ(1 + ν)/E. All the coefficients appearing in Eqs. 8B and 8C depend, in principle, on the cell density ρ. In particular, here, we assume the following simple linear dependenceS0=2AC=1ρcρ, α=α0(1ρρ0)

Here, ρc is the critical density associated with the isotropic-nematic phase transition, whereas ρ0 > ρc represents the density at which the cells are closely packed and able to transfer stresses on their surroundings. In addition, we assume the equation of state P = P0(ρ/ρ0 − 1).

The continuum model outlined here is closely related with the active gel model introduced by Kruse et al. (41) to explain the formation of defective patterns in polar actomyosin gels but further accounts for the remodeling caused by elastic deformations of the bacterial domains at short time scales. Experimental estimates of the viscoelastic crossover time τ are available in the literature on biofilms [see, e.g., (42, 43) for recent publications] but are less common for microcolonies, as the small size and the monolayer structure of these bacterial communities render standard rheological techniques ineffective. Yet, from (6, 19), one can infer τ to be of order 100 min for E. coli under physiological conditions in confined and freely expanding colonies respectively, thus comparable with our molecular dynamics simulations (see Fig. 4). Furthermore, Eqs. 7 and 8 do not account for rotational viscoelasticity, such as that encountered in certain type of polymeric liquid crystals, such as dense suspensions of actin filaments, and originating from the entanglement of these semiflexible polymers in three dimensions (44).

Equations 8A to 8C are numerically integrated using finite differences on a rectangular domain endowed with stress-free absorbing boundaries in the x direction and periodic boundaries in the y direction. The system is initiated with ρ=ρ0, v = 0, and a random configuration of the Q field and evolved until a steady state is found. More details about the numerical integration of Eqs. 8A to 8C, including all parameter values, can be found in the Supplementary Materials.


Supplementary material for this article is available at

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is not for commercial advantage and provided the original work is properly cited.


Acknowledgments: Funding: We are indebted to A. Sengupta for several illuminating discussions and for sharing with us the images displayed in Fig. 8. This work is partially supported by The Netherlands Organization for Scientific Research (NWO/OCW) as part of the Frontiers of Nanoscience program and the Vidi Scheme (Z.Y., D.J.G.P., and L.G.). Author contributions: Z.Y., D.J.G.P., and L.G. developed the models, performed the simulations, analyzed the data, and wrote the manuscript. Competing interests: The authors declare that they have no competing interests. Data and materials availability: All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. Additional data related to this paper may be requested from the authors.

Stay Connected to Science Advances

Navigate This Article