We model growing bacteria as two-dimensional expanding rods, whose diameter d ₀ is fixed and whose length l _i expands with a mean growth rate g up to length l ₀ (Figure 1A). When a rod reaches length l ₀ it splits into two rods of length l ₀/2. The growth of each rod with rate g _i is linear in time, and for each rod we choose g _i to be a random variable in the interval [g/2, 3g/2] to avoid synchronization effects [28]. Following mechanics-focused approaches used in recent computational studies of bacterial colonies [28, 34, 48–50], the rods’ dynamics are governed by an overdamped motion for the positions r _i and orientations e i = ( cos θ i , sin θ i ) T expressed by the angle θ _i . Explicitly, they are updated according to. d r i d t = 1 ζ l i ∑ j F i j , (1) d θ i d t = 12 ζ l i 3 ∑ j r i j × F i j ⋅ e z , (2) where ζ is a friction coefficient, e _z is the unit vector that is perpendicular to the plane of motion of the rods, r _ij = r _i − r _j is the separation vector between two rods, and F _ij are the forces between rods. The forces are modeled using a Hertzian repulsion, which reads F i j = F 0 d 0 1 / 2 h i j 3 / 2 n i j , (3) where F ₀ is the repulsion strength, h _ij is the overlap length of particle i and j, and n _ij is the normal vector of the particles’ closest point of contact (Figure 1B). In the following we use l ₀ as a unit of length, and g as a unit of velocity, implying a natural unit of time l ₀/g. Our model has two dimensionless parameters, which are the aspect ratio l ₀/d ₀ and gζ/(F ₀ l ₀); the latter we keep fixed to gζ/(F ₀ l ₀) = 5 × 10^–7. We simulate up to 25,000 rods from an initial state of a single rod.

To analyze the simulated colonies, we compute a nematic tensorial orientation field from our particle-based simulations. Following Ref. [32], we apply a smoothing function that transforms the ith particle into h i r = 1 4 tanh e i ⋅ Δ r i + l i / 2 σ − tanh e i ⋅ Δ r i − l i / 2 σ × tanh e i ⊥ ⋅ Δ r i + d 0 / 2 σ − tanh e i ⊥ ⋅ Δ r i − d 0 / 2 σ , (4) where Δ r _i = r − r _i is the separation from the particle’s position r _i ; e i ⊥ is the vector perpendicular to the particle’s orientation e _i ; l _i is the particle’s current length; and σ = 2d ₀ is the smoothing length. The nematic field is then given by Q μ ν r = ∑ i N h i r 2 e i , μ e i , ν − δ μ ν , (5) where the sum runs over particle indices i and μ, ν refer to the vector components, and the nematic order parameter is S = Q x x 2 + Q x y 2 . Similarly, we calculate a velocity field as v r = ∑ i N h i r v i , (6) where v _i = d r _i /dt is the velocity of an individual particle approximated from its displacement in the most recent timestep.

Before studying colonies containing many defects, we examine the early-time behavior of isolated nematic defects programmed into the colony’s initial condition. In hydrodynamic field theories of fixed-size active nematics, director fields containing isolated defects display very distinct vorticity patterns depending on whether the defect has positive or negative charge: a +1/2 defect generates two regions of opposite vorticity nearby, resulting in a net active force in the transverse direction, whereas a −1/2 defect is surrounded by a threefold-symmetric pattern of alternating vorticity that produces no net force at the defect core [14, 39].

To investigate whether this picture applies in the case of growing active nematics, we initialize our simulation with rod orientations approximating an individual defect. We then allow the colony to grow until each cell has doubled at least once, resulting in the configurations shown in Figures 2A,D for a −1/2 and +1/2 defect respectively. Figures 2B,E display the corresponding nematic field, exhibiting the typical structures of +1/2 and −1/2 defects.

The velocity field reveals a flow generally oriented radially outward from the center of the colony for both +1/2 and −1/2 defects (see Figures 2C,F), which is induced by the radial outward growth of the colony. We also computed the flow fields’ vorticity ω = ∂ v y ∂ x − ∂ v x ∂ y (color code in Figures 2C,F). In line with nematic field theories [14, 51] the +1/2 defect displays two regions of vorticity with opposite sign and the −1/2 defect shows a combination of six regions with alternating sign. These results for isolated defects lead us to expect that the +1/2 defect self-propulsion dynamics familiar from fixed-size active nematics will compete with the radial expansion flow in growing active nematics.

Typical simulation snapshots resulting from a single-cell initial condition are shown in Figures 3A–C, displaying the growth of a colony. Consistent with prior work [28], we find microdomains of rods with similar orientation. Neighboring microdomains meet at grain boundaries, many of which contain nematic point-disclination defects. Figure 3D shows a typical nematic field for a colony at the end of a simulation, where we mark nematic defects with +1/2 (magenta arrows) and −1/2 (red circles) charge. The defects are identified by integrating changes in orientation on a small circuit around every point in the nematic field, leading to a number of defect candidates around each (true) defect position, which is then found using a clustering algorithm.

As shown in Figure 3E, this method identifies defects reliably such that we can extract defect trajectories. These are shown in Figure 3G for +1/2 defects and in Figure 3H for −1/2 defects. For both defect types, we observe a predominant radial outward trend, implying that the expansion of the colony has a strong influence on the defect motion.

The absence of significant braiding among the +1/2 defect trajectories (Figure 3G) represents a notable distinction from fixed-size active nematics, where such braiding corresponds strongly to the active material’s chaotic flow field [37]. Relatedly, we find that Lyapunov exponents calculated for our system, using the velocity field of Eq. 6 (see also [37]), are statistically indistinguishable from zero, regardless of particle aspect ratio. Chaotic dynamics are characterized by nonzero Lyapunov exponents, as these measure the separation rate of two infinitesimally close trajectories.

As a measure of nematic order, we compute the orientational correlation function, commonly used for characterizing nematic liquid crystals [51], which reads g N r = ⟨ P 2 e i ⋅ e j ⟩ , (7) where P 2 ( x ) = 1 2 ( 3 x 2 − 1 ) is the second Legendre polynomial and r = | r _i − r _j | is the distance between particle i and j. As we would expect, g _N (r) decays slowly for large-aspect ratio particles, indicating nematic orientational correlations over some significant correlation length, whereas for smaller aspect ratios, a sharp decay of g _N near r = 0 reveals an absence of nematic order (Figure 3F).

We now examine the statistical properties of the nematic defects’ motion. First, we compute their velocities v and their speed distribution, shown in Figure 4A. We find that both types of defects have a peak at intermediate speeds, which is more prominent for +1/2 defects. Conversely, and surprisingly, higher speeds are somewhat more common for −1/2 defects. A similar picture arises when we project the defects’ velocities onto the radially (from the center of mass of the colony) outward direction at the defect’s position (Figure 4B). Again, both +1/2 and −1/2 defects show a peak at a positive value, and the −1/2 defects’ distribution is enhanced towards high radial velocities. That is, while colony growth generally pushes all defects farther from the colony center, outward radial motion is greater for −1/2 defects than for +1/2 defects.

However, for +1/2 defects, we have a second contribution to their motion: the asymmetry in the nematic field around the defect leads to an active motion or “self-propulsion”, a well-known phenomenon in fixed-size active nematics [7, 13, 14]. To determine whether, and how strongly, + 1/2 defects self-propel in growing active nematics, we compute each defect’s heading direction p , the orientation of the “comet-head” arrangement in the nearby director field, which is computed using the polarization p _μ = −∂ _ν Q _μν /|∂ _ν Q _μν | at the defect core. We then examine the projection of the defect’s velocity onto its heading direction. The defect’s velocity is computed with two independent methods, first from its trajectory and second from the velocity field v ( r ) at the defect’s location. The resulting distributions, shown in Figure 4C, both display a broad range of positive and negative values, but with a bias toward positive values. This implies that +1/2 defects move preferentially along their heading direction, the same defect self-propulsion trend found in fixed-size, extensile active nematics. We additionally find that the defect heading directions are isotropically distributed in the lab frame and have no significant correlation with the radial direction (Figure 4D).

The mean-square-displacement of defects (Figure 4E) shows a quadratic scaling with time, signifying ballistic motion, that stems from the radial expansion of the colony. Ballistic motion implies that the radial expansion dominates over the diffusive defect stirring dynamics typical of fixed-size active nematics [37, 39].

From the distributions in Figure 4, along with the defect trajectories in Figures 3G,H, the following conclusion emerges: While +1/2 defect self-propulsion is detectable in growing active nematics, the motion of +1/2 defects (as well as that of −1/2 defects) is dominated by the radial expansion flow.

We now turn to a population genetics perspective on mixing properties of our growing colonies, i.e. the degree of relatedness for spatially nearby cells. In order to quantify mixing of particles, we compute τ, the number of generations separating two cells phylogenetically, which is twice the number of generations since their most recent common ancestor (Figure 5A). Here, we first find the local neighbors of a particle i within a radius r cut = λ d 0 l 0 / π , where λ = 125 controls the size of the local neighborhood. We note that the results presented in the following do not qualitatively change upon varying λ within reasonable bounds. For each particle i we compute τ _i as the average of its τ measure with each of its neighbors. In Figures 5B,C we show τ _i for colonies growing up to aspect ratios l ₀/d ₀ of 2 and 10, respectively. Visual inspection reveals that rods with a higher aspect ratio have a higher τ _i and therefore mix more than at a lower aspect ratio.

On a global level, we determined the average phylogenetic distance τ ̄ = 1 N ∑ i τ i . As expected, τ ̄ increases with time (Figure 5D). Additionally, τ ̄ increases with aspect ratio. In the fully grown colony, we show τ ̄ as a function of aspect ratio in Figure 5E. We find that τ ̄ increases for small aspect ratios, implying an increase of genetic mixing as aspect ratio increases, until around l ₀/d ₀ = 7, where τ ̄ plateaus.

To make a connection to experiments [35, 40, 41] we performed simulations with two “alleles” (yellow and blue colors in Figures 6A,B) assigned randomly once the colony has reached N = 100 particles. Subsequently, cell division preserves the mother cell’s allele in both daughter cells, giving rise to the allele distributions shown in Figures 6A,B. As with the τ ̄ measure, the alleles show that a colony with aspect ratio 10 (Figure 6B) is better mixed than a colony with aspect ratio 2.

As a measure of genetic mixing in our two-allele representation, we compute the heterozygosity H, which, most generally, is the probability that two randomly selected organisms in the population share the same allele, but which can depend on distance r in spatially structured populations [41]. H is accessible to measurement in an experimental setting [35, 40, 41]. We study the small-distance limit of H(r) averaged over all cells in the colony, which we call H ₀, by considering all pairs of neighboring particles. With this restriction, the probability that two particles have different alleles equals the probability that both cells are on a line of contact between blue and yellow sectors, which recasts H ₀ as a geometric measurement. For our simulations we compute H ₀ = C _BY /A _B from configurations like Figures 6A,B, where C _BY is the contour length of contact between blue and yellow sectors (including single-particle “sectors”), and A _B is the total area occupied by the blue particles. Similarly to τ ̄ , larger values of H ₀ imply better mixing of the colony.

For a given Brownian dynamics simulation, we measure H ₀ for an ensemble of random blue/yellow color assignments when the colony has size N = 100 particles, propagating the colors forward in time by following the simulation’s phylogenetic tree, obtaining additional configurations similar to Figures 6A,B. These are used to compute the colony-averaged value of H ₀ and its associated standard deviation, shown in Figure 6C. We find that for smaller aspect ratios, heterozygosity increases with increasing aspect ratio, implying an enhanced genetic mixing. Comparing our findings for τ ̄ and H ₀ to the measurement of nematic correlation length (Figure 3F), we see that at low aspect ratios, colony self-mixing improves with increasing aspect ratio just as nematic order begins to emerge.

It is possible that simulations of larger systems would observe active nematic self-mixing more akin to the fixed-size case. The exponentially increasing computational expense with number of generations favors continuum models at these larger scales. Agent-based modeling such as the results we have presented can help to constrain and validate continuum hydrodynamic models for these purposes.