Stress anisotropy in confined populations of growing rods

A central feature of living matter is its ability to grow and multiply. The mechanical activity associated with growth produces both macroscopic flows shaped by confinement, and striking self-organization phenomena, such as orientational order and alignment, which are particularly prominent in populations of rod-shaped bacteria due to their nematic properties. However, how active stresses, passive mechanical interactions and flow-induced effects interact to give rise to the observed global alignment patterns remains elusive. Here, we study in silico colonies of growing rod-shaped particles of different aspect ratios confined in channel-like geometries. A spatially resolved analysis of the stress tensor reveals a strong relationship between near-perfect alignment and an inversion of stress anisotropy for particles with large length-to-width ratios. We show that, in quantitative agreement with an asymptotic theory, strong alignment can lead to a decoupling of active and passive stresses parallel and perpendicular to the direction of growth, respectively. We demonstrate the robustness of these effects in a geometry that provides less restrictive confinement and introduces natural perturbations in alignment. Our results illustrate the complexity arising from the inherent coupling between nematic order and active stresses in growing active matter, which is modulated by geometric and configurational constraints due to confinement.

Self-organization in multicellular biological systems is driven by inherent cellular activity.This non-equilibrium activity can take a variety of forms, including selfpropulsion from cell motility [1,2], active adhesion [3] or chemical activity from metabolism or signaling [4,5], which also occur in non-biological active matter [6,7].Growth is one of the hallmarks of life.It constitutes another process by which energy can be injected into the system at the microscopic scale, and has been shown to able to balance out chemical interactions at the level of large-scale behaviour [8].The mechanisms by which growing active matter self-organizes to form multicellular communities such as biofilms [9,10] and functioning tissues [11][12][13] are complex and often involve other forms of activity or internal regulation [14][15][16][17].Here, we focus on the mechanical aspects of growth, mediated by steric interactions between individual rod-shaped cells and confinement.They are sufficient to reproduce alignment and large-scale flow patterns observed in the initial stages of bacterial colony formation [18][19][20][21] and represent a prime example for the class of systems known as active nematics [22].
The consequences of these ingredients depend critically on the mechanical environment: For freely expanding colonies in two dimensions, the overall colony shows no preferred direction but generates locally ordered microdomains arising from a competition between passive elastic properties and active extensile stresses generated by growth, which are themselves functions of cell properties such as their length-to-width aspect ratio and growth rate [19,23].Topological defects typical for active nematics, which arise naturally at the domain boundaries [23,24], have also been shown to be important for the transition into the third dimension [25,26].
In contrast, both in nature and lab experiments, colonies often grow in confined environments.A prominent example are rectangular channels, where populations of growing rod-shaped cells develop global nematic order, with their elongation direction oriented towards the exits (see Fig. 1A).Despite the simple setup, these systems show a plethora of phenomena which have served to elucidate fundamental processes governing growing active nematics in a number of studies.In the course of these investigations, the emergence and maintenance of global orientational order were attributed to the response of the elongated particles to growth-induced expansion flow [27] as well as to globally anisotropic stresses which build up due the distinct boundary conditions in different directions [28].Similar to the local order characterized by the microdomain size in freely expanding colonies, the orientational order parameter for this global alignment was found to depend on the mechanical parameters of the system, such as the length-to-width ratios of the growing rods and the growth rate, the latter also confirmed in experiments [29].For near-perfect alignment, the system also becomes inhomogeneous in time and is characterized by intermittent breakdowns of order due to a buckling instability [30,31].
The picture of anisotropic stress maintaining global order is appealing, as it can be directly related to the geometry of the system: In the direction of the openings, stress can be dissipated more easily than in the confined direction perpendicular to it and this perpendicular stress can help stabilize the emergent order, leading to a self-consistent, dynamic steady state.While this phenomenology was reproduced by a visco-elastic continuum theory and was also observed without the perpendicular confinement [28], the buckling analogy [30] (which is also related to anisotropic stress) and resulting dynamical complexity [31] suggest that the interaction between active and passive stresses becomes more complex near perfect order.
In this work, we further investigate this limit.Using particle-based simulations and coarse-graining the resulting spatially varying stress tensor fields, we show that the picture of an excess stress perpendicular to the growth direction required to stabilize alignment breaks down.Instead, we find that the system can spontaneously organize into states with a long-term average stress anisotropy which is inverted, i.e., showing excess stress in the direction of growth.Starting from asymptotic theories for perfect alignment, we show that the resulting stress tensor is consistent with a decoupling of active and passive stresses, which persists even in the case of open channels, where confinement is not provided by perfectly straight walls but by the fluctuating growing colony itself.

A. Agent-Based Modeling
We use an agent-based model of compressible rod-like cells that grow in length and divide.Each rod is comprised of a rectangular fixed-width body with half-circle caps at both ends.Over time, the total cell length increases linearly to l max , twice its starting value, and then divides into two daughter cells with identical orientation ϕ and randomly drawn growth rate from a uniform distribution γ ∈ 3 4 , 5  4 .Neighbouring rods interact mechanically only and repel one another according to Hertzian repulsion law [32] as illustrated in Fig. 1B, where d is the shortest connection between the two center lines.The simulations are set in the over-damped limit, such that the velocity of the ith particle v i = µ µ µ µ µ µ µ µ µ µ µ µ µ µ µ µ µ i F i , where F i is the total force on particle i and µ µ µ µ µ µ µ µ µ µ µ µ µ µ µ µ µ i is an anisotropic mobility tensor.The model is very similar to other models of dividing rods in the literature [27,28,30,31] and builds on physical intuition, hence, we refer the reader to the Supplementary Material for additional details.

Confinement
We consider a two-dimensional channel with outlets on the two opposing sides in the y direction and confinement on two sides in the x direction.Particles are removed once their center point crosses the outlet boundaries.Confinement in x direction is modeled by walls that exert Hertzian forces as in Eq. (1) on the cells, where d is now the shortest connection from the center line of the cell to the wall.Two examples of such a channel simula- tions on a 50×50 unit domain are shown in Fig. 1C.In our actual investigation, we will use domains of size 200×200.All lengths are given in multiples of the rod width, which is kept constant throughout this work and measurements are taken after discarding transient dynamics.

B. Extraction of Continuous Fields
To study the emergent dynamics at greater length scales, we compute approximate continuum fields of the orientational order parameter and stress tensor on a fine rectangular grid.We write for the former where a (k,∆V ) is the overlapping area of particle k and small sub volume ∆V around point r and the normalization factor Z(r) = k a (k,∆V ) corresponds to the total overlap area of all cells with ∆V .The average orientation is given by where arg refers to the angle in the complex plane and we divide by two to account for the nematic symmetry.The scalar order parameter within the same region is given by where ξ = 1 corresponds to complete alignment and ξ → 0 to no order.The stress tensor field is computed as [33] where F kl are the cell-cell interaction forces and 0 ≤ λ kl ≤ 1 is fractional length of the line segment r kl = r k − r l that lies within ∆V .Averages across larger areas or along an entire axis can then be computed as required.Throughout this study, absolute stress values are reported in units of S 0 , for which we arbitrarily choose the central |σ yy | for l max = 6 in a domain with our standard height of 200 length units (see Supplementary Material for an explicit definition of S 0 ).A useful quantity introduced by You et al. is the normalized stress anisotropy [28] ∆Σ = which quantifies the excess stress across the channel, i.e., perpendicular to the walls, in the direction of confinement [34].

A. Spatially Varying Fields
We set up the numerical experiments as described above by placing four randomly oriented rods in the center of a 200 × 200 unit channel with confining walls at the sides and let the dynamics evolve for many cell division times.Examples of this process are shown for a smaller domain and division lengths l max = 2 and l max = 6 in Fig. 1C with snapshots taken after 7, 9, and 15 generations, respectively.An important observation is that preferential alignment is always directed towards the channel outlets.However, the final, highly ordered, columnar structure with ξ ≈ 1 as in the second example is only attained for l max 4, whereas for l max 4, the steady-state is characterized by imperfect order with ξ < 1.
For rods with l max = 2 (resulting in ξ ≈ 0.5), we resolve the stress field of the final steady state along the channel.We find that both |σ yy | and |σ xx | take on a parabolic profile with |σ xx | being consistently larger as shown in Fig. 1D (top).The system is therefore, both locally and globally, characterized by ∆Σ > 0, which has previously been observed as typical [28], along with the intuitive explanation that stress is more difficult to dissipate in the confined direction and needed to generate and maintain a preferential alignment towards the channel outlets.
However, cells dividing at a greater maximal length l max = 6 (resulting in ξ ≈ 1) exhibit a different behaviour as shown in Fig. 1D (bottom): Here, measurements of the stress tensor along a central cut through the channel show that, while the vertical stress |σ yy | follows the familiar parabolic profile, the confined horizontal stress locks onto a constant value in large parts of the domain.Surprisingly, this even leads to an inversion of the stress anisotropy ∆Σ in the colony center, with the horizontal stress |σ xx | remaining at lower values than the vertical stress |σ yy |.
Observing near uniform values of |σ xx | while |σ yy | changes continuously indicates a decoupling of the stresses made possible by perfect ordering inside the channel.We will approach the study of this (de-)coupling from two angles: We begin by deriving limiting theories which reveal the physical origin of |σ yy | and |σ xx | by predicting them separately.We then continue numerically by modifying the geometry, primarily impacting constraints on the horizontal stress.

B. Column Theory for Active Stress
A colony of our model rods has a well-defined and spatially homogeneous distribution of cell age g ∈ [0, 1) and growth rate γ ∈ 3 4 , 5 4 described by with a normalization parameter χ and averages γ ≈ 0.985 and g ≈ 0.441.
With this distribution at hand, we predict the parabolic stress profile using the emergent and effectively onedimensional columnar structures.In this limit, we neglect all transverse (horizontal) dynamics and derive an expression for the |σ yy (y)| stress profile within the colony.Also assuming incompressibility, a generic dry continuum ∂ t ρ + ∇ • (vρ) = αρ for the cell density ρ with an effective growth rate α immediately leads to the condition ∇ • v = α for the steady state, or ∂ y v y = α for our onedimensional consideration.The effective rate α at which line density is produced in our model can be calculated from Eq. ( 5), such that where ∆l = lmax 2 is the length (including caps) by which a cell grows between divisions and l = ∆l 1 + g is the average length occupied by a cell.By defining y = 0 as the center of the colony with v y (y = 0) = 0, the velocity profile therefore becomes where the second equality represents force balance and µ is an effective mobility that accounts for substrate friction, which, in general, could be a function of density and therefore of y.Again using the incompressibility assumption, it is computed as where, in the agent-based model, the mobility parallel to the symmetry axis µ || is a function of the cell aspect ratio l 2R with the cell radius R, as defined in the appendix model description, and the average can be computed numerically.With the boundary condition that the stress vanishes at the outlets, σ yy (y max ) = 0, we integrate Eq. ( 6) and obtain a parabolic stress profile σ yy (y) = γ µ 1 + g (y max − y) 2  2 with the maximum stress at the center.This result is used in Fig. 1D to predict the stress profile and matches closely for the case of l max = 6 while significantly underestimating the central stress of l max = 2 due to compression effects neglected in the theory.A more detailed comparison is shown in Fig. 3A, where the maximal stress |σ yy | is computed for a range of division length values l max .The theory matches numerical results for long cells, the limit for which the theory was derived, as it yields the near-perfectly ordered quasi-columnar structure.For shorter cells, disorder and compression becomes more relevant, causing the theory to underestimate the measured stresses.We conducted additional simulations with a single one-dimensional colum of cells while freezing the orientational degree of freedom (grey line in Fig. 3A).The fact that this measurement matches the theory even more closely confirms that a large part of the deviation is caused by transverse compression, which effectively increases the density and thereby also changes the mobility.
It may initially seem surprising that the theory approximately captures |σ yy | even at small division lengths l max 4, where the system self-organizes into a weakly ordered state far away from a columnar structure assumed above.However, on a mean-field level, the start- ing point of our incompressible theory, including v x = 0, ∂ y v y = α and the distribution of cell ages, holds in the disordered system as well.Therefore, the theory can also be viewed as an approximation to the disordered 2D scenario.An additional inaccuracy, besides the violations of the incompressibility assumption explained above, then arises from the anisotropic mobility of the particles, which changes the estimate of µ in Eq. ( 7).

C. Passive Stress Theory
To determine whether active and passive stresses indeed decouple in the highly ordered state, we calculate the expected |σ xx | arising only from the passive repulsion between neighbouring columns.This is possible in such a state since the overlap ∆ < 2R between neighboring columns can be calculated from the number of columns, the channel width and the width 2R of a single cell.Parallel cells at distances ∆ from one another then exert Hertzian repulsion forces of the form on each other.
We then consider the line density of these force vectors along the vertical axis: Each cell pair has at most one interaction, but depending on the cell length, which ranges between l max /2 and l max , and configuration, some cells may interact with up to three cells on each side, as illustrated in Fig. 2.
Considering two adjacent columns, we can count the number of unique touching cell pairs.Every new cell in the left column has one interaction partner on the right and likewise every new cell on the right interacts with the one to its left as is indicated by red arrows in Fig. 2A.This scheme of counting captures all interactions and we write 4R l since on average there is a new cell every l in each column.The rightmost term represents a correction factor: As illustrated in Fig. 2B, whenever the cell caps are too close to those in the other column, one of the two potential interactions is lost.Assuming there are no longrange correlations between the neighbouring columns, two cell cap regions are expected to overlap with probability p = 4R l −1 .A first comparison with numerical estimates is shown in Fig. 1D (bottom) where a dashed line indicates the expected configurational stress |σ xx | for 205 columns in the 200 unit wide channel using l max = 6.The measured value itself is matched well by our theory, again with a slight underestimation due to neglected compression along the columns, which increases the line density of cells and therefore of the interactions.

D. Statistically negative stress anisotropy
An important input for the prediction is the exact number of columns in the configuration, which is not predetermined by the parameters of the system.Instead, it is an an emergent result of the self-organization process that leads to perfect alignment and, in principle, can also fluctuate over long time scales [31].Therefore, the question arises whether the inverted anisotropy |σ xx | < |σ yy | in Fig. 1D (bottom) is typical, or, more generally, what the distribution of the number of columns is that the system naturally attains.
To evaluate this in detail, we generated channel simulations with l max = 6 for slightly varying domain width.The result is shown in Fig. 3C.The allowed discrete configurational stress values for a fixed number of columns change continuously with increasing domain width (dashed lines).For each value of the domain width, 20 independent initial conditions lead to a spectrum of observed of column numbers, whose stress values show excellent agreement with the corresponding discrete levels.While a few simulation runs generated a state with |σ xx | > |σ yy |, the vast majority exhibit ∆Σ < 0. This is further emphasized by the marginal distributions shown in Fig. 3C.

E. Open Domains
The channel geometry used so far is in many ways special: Not only does it impose a preferential direction for growth (or, in fact, allows for only one possible flow direction), it also has perfectly straight walls which ensure that the same amount of space is available to each column as cells are pushed along the channel.Furthermore, the walls do not introduce any perturbations.As we have seen, the passive horizontal stress is purely determined by these boundary conditions and is therefore also constant along the channel.Given that the growth-induced active stress nevertheless drives the self-organization process that selects a certain number of columns, it may seem as if the special properties of the channel confinement alone are responsible for the decoupling of the two stresses and other unusual features of this state, such as the negative stress anisotropy ∆Σ < 0. To find out whether this is the case and better understand the origin of emergent effects, we carry out simulations in rectangular domains that have outlets on all four sides.Open domains provide a weaker constraint on the dynamics as they permit, a priori, flow in all directions.Still, for non-square domains, we find preferential alignment between the long opposing sides, as shown in Fig. 4A.From the initial unordered phase, the dynamics produce a highly ordered central region supported by disordered sides that generate the required compression.This is a non-trivial result in itself and highlights the system's natural tendency to generate order.However, due to the complex dynamics in the side regions, the ordered phase may spontaneously become unstable leading to a macroscopic buckling event before the ordered state is eventually reestablished.An example of this is shown in Supp.Fig. 4.
Measurements of the stress tensor fields are shown in Fig. 4B-4E.Analogous to Fig. 1D, these plots display stress profiles as measured along a vertical cut through the domain center (Fig. 4B and 4C) and, in addition, the same along the horizontal direction (Fig. 4D and 4E).For short cells, the behaviour along the vertical shown in Fig. 4B is the same as in channels and along the horizontal axis (Fig. 4D) we find a profile that flattens towards the center.In the case of longer cells, as displayed in Fig. 4C and 4E, the behaviour deviates from the channel simulations.Most importantly, the open domains do not allow for spatially homogeneous passive horizontal stress.Instead, horizontal compression is actively generated near the domain sides which can be seen from the steep increase in stress along the horizontal cut in Fig. 4E.Still, the dynamics in the center are quite similar to the channel, as no horizontal motion is apparent and no active stress is generated horizontally due to vertical nematic ordering (compare snapshot in Fig. 4A and the flat part of the horizontal profile Fig. 4E).Therefore, the horizontal stress in the central region can again be expected to be mostly of passive origin, due to the compression of vertical cell columns.However, the fact that the horizontal stress varies vertically in Fig. 4C (in contrast to Fig. 1D bottom) means that the space that columns occupy is not constant as in the channel, but increases from the middle towards top and bottom.Therefore, the open geometry successfully removed this artificial constraint imposed by the channel walls and, at the same time, introduced perturbations from the disordered sides, lending additional weight to the observation, that, nevertheless, clearly a negative anisotropy |σ xx | < |σ yy | is observed on average.

F. Length dependency
So far, we have only considered exemplary cell division lengths l max to understand the qualitatively different dynamics of short and long cells.It is therefore natural to ask how the results compare quantitatively across the entire l max range and between channels and open domains.A finely sampled length scan is summarized in Fig. 5, displaying the average order parameter ξ (Fig. 5A) and stress anisotropy ∆Σ (Fig. 5B) for both geometries.These values are computed in a central 100 × 20 box around the center of the domains.Averaging is done over 150 instantaneous measurements sampled from 15 generations of time evolution and additionally over 10 independent initial conditions.While always being slightly more disordered, the dynamics in the open domain centre closely resembles that of the channels across the full parameter range.In both cases, the order approaches unity for division lengths l max 4 and this is accompanied by a transition to negative anisotropy ∆Σ.Interestingly, even though the order is slightly reduced in the open domains for all l max , they exhibit a cleaner relation of l max and ∆Σ as shown in Fig. 5B.Both effects mostly likely have the same explanation, namely, that the lack of hard boundaries in combination with the larger open domain reduces finite-size effects in the column number.

III. CONCLUSION
We employed agent-based simulations to study the growth of bacterial colonies in confinement.This topic has seen a range of publications in recent years, but a careful analysis of the interaction forces and the corresponding spatially resolved stress fields has revealed previously unreported phenomena.
To understand the dynamics from a theoretical perspective, we focused on the limiting case of perfect ordering.In this limit, it was possible to predict the (active) vertical stress field.It gives reasonable predictions even for short cells that do not generate the column structure used in the derivation.
The column structure also allowed for an estimation of the configurational horizontal stress.Comparison with simulation data confirmed that, in the column structure, the horizontal stress is indeed largely passive with the theoretical prediction matching measurements up to a few percent.
Many of the qualitative features of these states were also found in simulations on open rectangular domains, which revealed that, even without hard and flat walls, near-perfect order and strong negative stress anisotropy may emerge.This is an important result, as it indicates that neither strong ordering nor the negative stress anisotropy were an artifact of the particular confining channel, but rather a natural emergent property of the dynamics.It provides a compelling argument that the ordered state can indeed be stable in the presence of a (limited) stress anisotropy inversion.Away from perfect order, there is a tight coupling of the stress field components and order, as any imbalance in the interaction forces causes microscopic reorientation away from the direction of strongest compression.This feeds back to the stresses, both through relaxation of the existing stress, as well as by redirecting the active stress of growth.However, as the order parameter approaches unity, compression along the principal axes of the orientation field ceases to generate torques on the rods involved and the components of the stress tensor decouple.
An interesting connection of our work to that of You et al. [19] is that, in our open domains, we observe the characteristic formation and breakup of microdomains.The rectangular geometry formed by outlets is sufficient to drive the described ordering both instantaneously and in the long-time average.However, at intermediate timescales, stochastic buckling events are observed that break up even the macroscopically ordered domain in the middle region into smaller fragments.The frequency of these events depends strongly on the relative scales of the total domain size versus typical microdomains which themselves depend on the division length l max .This effect required the use of many independent initial conditions to produce Fig. 4E without significant fluctuations in the |σ xx | profile.The statistics of these events in open domains and the connection to the smectic properties of the column structure seem to warrant more detailed future investigations (as done for the channel geometry [31]).
Dell'Arciprete et al. [23] explain the underlying alignment mechanism of the emergent dynamics using torques experienced by rods in shear flow, similar to recent observations in other geometries [21].You et al. [28] on the other hand attribute this to an anisotropy in the stress tensor.One could, of course, argue that these two are not very different at all, given that the stress field creates flows in the first place.However, our work adds to this discussion.We, unexpectedly, found an inversion of stress anisotropy in the strongly ordered limit, even in open domains.This is not consistent with a purely stress based coupling.On the other hand, the observed stochastic collapsing of the ordered structures in the open domains can not be explained using flow based arguments only, as it lacks a destabilizing mechanism for this case.A theory combining both effects may be a promising avenue for capturing all observed phenomena.
In this work, we presented a detailed analysis of a computational model of growing rod-shaped particles and discovered states in which activity-induced stresses, passive nematic properties, volume exclusion and confinement can interact in novel ways and yield strongly anisotropic systems.Despite these emergent intricacies, our model merely captures one specific physical aspect of reality when it comes to the growth of bacteria in colonies or biofilms.It would be interesting to elucidate what role the effects discovered here can play among the complex processes that characterize natural systems.This includes nutrient distribution [16,35,36], which leads to gradients of activity and thereby represents another opportunity for the passive properties of the system to become important, but also extends to structure formation in three dimensions [20,35,36], interactions with other types of activity such as motility and chemical signal production [37] and different kinds of confinement [38].In addition, it may prove rewarding to look for even more general physical principles by exploring the similarities the system exhibits to dense systems of active polymers and filaments [39,40] due to its columnar structure.

CONFLICTS OF INTEREST
There are no conflicts to declare.

FIG. 1 .
FIG. 1. (A) Experimental observation of growing E. coli exhibiting alignment when confined in a quasi-2D microfluidic trap.The openings of the trap are at the top and bottom (bright refracting areas), the confining walls are towards the left and the right (outside the area shown here).(B) Illustration of Hertz-based repulsion force law.(C) System configurations in a 50 × 50 channel after 7, 9, and 15 cell generations for parameter lmax = 2 (top row) and lmax = 6 (bottom row).Colouring indicates orientation as defined in the colour wheel.(D) Cuts along the length of a simulation channel with height/length 200 and filled with lmax = 2 (top) and lmax = 6 (bottom) growing rods.Theoretical predictions (dashed) accompany simulation results.

FIG. 2 .
FIG. 2. (A) Illustration of counting argument for passive stress theory.(B) Detailed view of missing interaction.

FIG. 3 .
FIG. 3. Evaluation of theoretical predictions.(A) Central ystress in channels (blue) with theory (dashed) and a simplified 1D column simulation (grey).(B) Stationary central stress field components for cells with lmax = 6.|σxx| and |σyy| are measured in a channel-wide 10 unit tall slice (as illustrated in Supp.Fig. 3) and averaged over five generations.Predictions for various column numbers are shown as dashed lines.(C) Marginal distributions of points in Figure 3B.

FIG. 4 .
FIG. 4. Open rectangular domains.(A) A snapshot of a 600 × 200 unit open domain filled with lmax = 6 growing rods.Colour indicates rod orientation as illustrated in the colour wheel.(B-E) Central vertical and horizontal cuts through the domain measuring the local stress tensor coefficients for lmax = 3 (top) and longer lmax = 6 (bottom) cells.Approximations were computed as illustrated in Supp.Fig. 3 using 10 initial conditions for lmax = 3 and 210 initial conditions for lmax = 6 to better capture increased fluctuations.

FIG. 5 .
FIG. 5. Direct comparisons of channels and open domains under varying cell division length lmax.(A) Order parameter in a central 100 × 20 region as illustrated in Supp.Fig. 3 averaged over time and initial conditions.(B) The similarly averaged central anisotropy for channels and open domains.