Mechanics and morphology of proliferating cell collectives with self-inhibiting growth

The morphology of a cell collective is affected by factors across scales, from the mechanics and sensing of individuals to system-wide environmental conditions. At the local level, cell elasticity and geometry influence packing configurations Farhadifar et al. (2007); Staple et al. (2010); You et al. (2018), intermittent phases of growth and motility drive concentric density variations Shimada et al. (2004); Liu et al. (2011), and chemical cues give rise to internal patterning and interfacial instabilities Adler (1966); Budrene and Berg (1991, 1995); Lu et al. (2022); Martínez-Calvo et al. (2022). Global factors such as fluid flows Atis et al. (2019); Pearce et al. (2019); Benzi et al. (2022) or confinement Helmlinger et al. (1997); Yan et al. (2021); Amchin et al. (2022) can also modify growth rates and interfacial dynamics. Essential to all of these features is cell proliferation, in which individuals produce biomass from the consumption of chemical nutrients Hallatschek et al. (2023).

In contrast to conventional active matter systems, such as particle suspensions whose configuration is specified by spatial and orientational degrees of freedom, systems undergoing proliferation possess additional degrees of freedom that might include, for example, age or size Rubinow (1968); Lebowitz and Rubinow (1974); Wittmann et al. (2023). Another distinguishing feature of proliferating systems are growth-induced stresses, which are particularly significant in cell collectives due to the inevitable many-body collisions between cells. These stresses not only cause local cellular rearrangements You et al. (2021), but may also impact growth rates Delarue et al. (2016); Alric et al. (2022).

The influence of stress on growth and morphology has been analyzed in both discrete and continuum settings Shraiman (2005); Gniewek et al. (2019); Schnyder et al. (2020); Carpenter et al. (2024); Ye and Lin (2024). Less attention has been paid to the coupled dynamics of growth-induced stress and microstructure. In this work, we study how stress-dependent growth affects the spatiotemporal dynamics and structure of proliferating colonies of elongated particles, such as bacteria or yeast. Many-particle simulations reveal concentric wave-like patterns in cell size whose space and time scales are determined by the sensitivity of individual cell growth to local mechanical stress. Based on this discrete model, we derive a continuum theory that couples the microscopic size degrees of freedom to the macroscopic mechanics, and find exact solutions in radially symmetric geometries. These solutions identify the origin of cell size patterning as the result of cumulative resistance to growth across generations.

Discrete models of cell proliferation are invaluable for their precise control of environmental parameters and mechanical responses. In such models, contact resolution is a common challenge owing to the typically dense packing of cell bodies. Here, we employ a contact resolution algorithm that enforces a no-overlap constraint by means of a nonlinear convex optimization problem. This method builds on contact resolution algorithms originally developed for granular matter and computer graphics Stewart and Trinkle (1996); Anitescu and Potra (2002); Mazhar et al. (2015); Macklin et al. (2019) that have recently been adapted to inertialess active matter systems Yan et al. (2019, 2022). This constraint-based formulation avoids numerical stiffness associated with repulsive potentials, allowing for larger time steps at the cost of an iterative solution.

We consider a collection of $N(t)$ sphero-cylindrical particles of constant cross-sectional diameter $b$ that grow from an initial length $\ell_{0}$ to length $2\ell_{0}$, at which point they divide into two particles of length $\ell_{0}$. For simplicity, we neglect variations in cell size after division Amir (2014); Taheri-Araghi et al. (2015). Each particle is characterized by its center of mass $\mathbf{x}_{n}(t)$, orientation angle $\theta_{n}(t)$, and length $\ell_{n}(t)$, shown schematically in Fig. 1. Particles within the ensemble move and rotate under the contact forces created by growth. Considering particles atop a two-dimensional substrate with constant friction factor $\xi$, the translational and angular velocities induced by the total contact forces $\mathbf{F}_{n}^{c}$ and torques $T_{n}^{c}$ are $d\mathbf{x}_{n}/dt=(1/\xi\ell_{n})\mathbf{F}_{n}^{c}$ and $d\theta_{n}/dt=(12/\xi\ell_{n}^{3})T_{n}^{c}$, respectively. Here $\mathbf{F}_{n}^{c}$ and $T_{n}^{c}$ are expressed in terms of non-negative Lagrange multipliers determined from the constraints that particles cannot overlap, non-contacting particles do not generate force, and collisions minimize the rate of frictional dissipation (see Supplementary Material). Between divisions, each particle grows according to $d\ell_{n}/dt=(\ell_{n}/\tau)\exp(-\lambda\sigma_{n}^{c})$, where $\tau$ is a growth timescale, $\sigma_{n}^{c}$ is the compressive stress acting along the long axis of the particle, and $\lambda$ is the stress sensitivity parameter. For uninhibited growth ($\lambda=0$ or $\sigma_{n}^{c}=0$), individuals grow exponentially fast, which has been reported experimentally for several types of rod-like bacteria, though this is still debated Godin et al. (2010); Campos et al. (2014); Kar et al. (2021). For stress-sensitive growth ($\lambda>0$), the growth rate decreases with the compressive stress, an effect which is also observed in live cells Minc et al. (2009); Delarue et al. (2016); Alric et al. (2022). Like the forces and torques, the compressive stress $\sigma_{n}^{c}$ is expressed in terms of the contact Lagrange multipliers, which themselves depend on growth, resulting in nonlinear coupling between contact and growth when $\lambda\neq 0$. Further details on the formulation and its numerical implementation can be found in the Supplementary Material.

We non-dimensionalize the equations of motion by the cell length $\ell_{0}$, growth time scale $\tau$, and stress scale $\tau/\xi\ell_{0}^{2}$. The sole free parameter is the dimensionless growth sensitivity $\lambda^{\prime}=(\tau/\xi\ell_{0}^{2})\lambda$. The simulations begin with a single cell of unit length and diameter $b=\ell_{0}/2$, and run until the colony radius reaches $R=150\ell_{0}$, corresponding to $\mathcal{O}(10^{5})$ particles. Each panel in Fig. 2(A-C) shows five snapshots, equally spaced in time, of simulations for three different values of $\lambda^{\prime}$, with the particles colored by their length. As $\lambda^{\prime}$ increases, concentric waves form about the colony center. Their characteristic wavelength decreases with $\lambda^{\prime}$ up to a point where it approaches the particle scale and is difficult to distinguish by eye.

To gain insight into these patterns and their dynamic development, we formulate a continuum theory that describes both the collective macroscopic mechanics and the microscopic distribution of particle sizes. Let $\Psi(\mathbf{x},\ell,t)$ denote the number of particles with length $\ell\in[\ell_{0},2\ell_{0}]$ in a control volume centered at $\mathbf{x}\in\Omega(t)$, where $\Omega(t)$ is the colony domain. This distribution function satisfies a Smoluchowski equation in the material frame,

where $D/Dt=\partial/\partial t+\mathbf{u}\cdot\mathbf{\nabla}$ is the material derivative with the cell velocity field $\mathbf{u}(\mathbf{x},t)$, and $\dot{\ell}(\mathbf{x},\ell,t)=(\ell/\tau)\exp(-\lambda p)$ is a growth flux, identical to the single particle growth rate in the discrete model with the compressive stress $\sigma^{c}$ replaced by an isotropic pressure $p(\mathbf{x},t)$. In this model, cell division is described by a boundary condition on $\Psi$ with respect to $\ell$. Specifically, the boundary condition balances the number density at the initial and division lengths, $(\dot{\ell}\Psi)\big{\rvert}_{\ell=\ell_{0}}=2(\dot{\ell}\Psi)\big{\rvert}_{\ell=2\ell_{0}}$, where the factor of two comes from the assumption that cells divide into two daughter cells. A similar distribution with a division boundary condition was originally formulated to describe cell maturation Rubinow (1968), and was recently used to describe particle length distributions in the context of dynamic density functional theory Wittmann et al. (2023).

The macroscopic mechanics are described by a Darcy-like model, which has been applied to collective cell growth in many contexts Greenspan (1976); Roose et al. (2007); Lowengrub et al. (2009); Wu et al. (2013); Martínez-Calvo et al. (2022); Ye and Lin (2024). In the present context, this model arises as the minimizer of frictional dissipation between the particles and the substrate subject to a volumetric growth constraint, analogous to the principle satisfied by the particle simulations (see Supplementary Material). Specifically, the constitutive equations are

Here $\xi$ is a constant friction factor as in the discrete model and $g=[\dot{\ell}]_{\ell}/[\ell]_{\ell}=(1/\tau)e^{-\lambda p}$ is the volumetric growth rate with $[f]_{\ell}:=\int_{\ell_{0}}^{2\ell_{0}}f\Psi~d\ell$. We non-dimensionalize the continuum equations using the same characteristic scales as before. The sole free parameter is again the dimensionless growth sensitivity $\lambda^{\prime}=(\tau/\xi\ell_{0}^{2})\lambda$, which exactly corresponds to that in the discrete model. In the remainder of this work, and in all figures, we assume all variables are dimensionless and omit primes. Further details on the model and its derivation can be found in the Supplementary Material.

Because $g$ does not depend on $\Psi$, the system of Eqs. (2)-(3) can be solved independent of the Smoluchowski equation (1), after which Eq. (1) can be solved along spatial and length characteristics. From this solution, given a Dirac delta initial condition ${\Psi_{0}(\mathbf{x},\ell)=\delta(\ell-1)}$, the average cell length $\overline{\ell}:=[\ell]_{\ell}/[1]_{\ell}$ can be computed explicitly as

where ${\rm frac}[z]\in[0,1)$ is the non-integer part of $z\geq 0$ and $\bm{\Phi}(\mathbf{x}_{0};t):\Omega(0)\mapsto\Omega(t)$ is the Lagrangian flow map. For uninhibited growth, this gives the homogeneous solution $\overline{\ell}(\mathbf{x},t)=2^{{\rm frac}[t/\ln 2]}$, implying cells divide synchronously at time $t=n\ln 2$ for integer $n$, as expected. In contrast, for stress-sensitive growth inhomogeneities may arise through anisotropically accumulated resistance along characteristics. Notably, inhomogeneity is not necessarily a result of instability.

Motivated by the particle simulations, we consider radially symmetric colonies of radius $R(t)$ with $R(0)=\sqrt{b/\ell_{0}\pi}$, the initial size of which is chosen so that the initial area of the colony is approximately that of a single particle. In this geometry, Eqs. (2)-(3) can be solved analytically in $r=|\mathbf{x}|\leq R(t)$ for the pressure $p$ and radial velocity $u$. For $\lambda=0$ the solution is $p(r)=(R^{2}-r^{2})/4$ and $u(r)=r/2$, while for $\lambda\neq 0$ the solution is

where $c(R)=[(1+\lambda R^{2}/2)^{1/2}-1]/2\lambda R^{2}$. Figure 2(D-F) shows the average cell length $\overline{\ell}$ from the continuum theory for the same parameters as the discrete simulations, without fitting. Similar concentric waves form whose wavelength closely follows that of the discrete simulations. At the largest sensitivity $\lambda=10^{-2}$ [Fig. 2(F)], the wavelength is roughly $5\ell_{0}$, suggesting a wave-like structure underlies the seemingly random distribution in the discrete simulation of (C).

The analytical solution given in Eq. (5) provides various predictions for the macroscopic dynamics and their dependence on $\lambda$. Figure 3(A) shows the radially averaged stress from the discrete simulations $\langle\sigma_{n}^{c}\rangle$ (solid), binned in annuli of equal area, along with the analytical solution for $p$ (dashed). To reduce noise, the discrete stress is averaged over a time interval during which the colony grows by less than one particle length. As $\lambda$ increases, we find the overall stress decreases, with an order of magnitude difference between the $\lambda=0$ and $\lambda=10^{-2}$ cases. The continuum stress is typically lower than the discrete stress, likely due to the coarse-graining assumption that particles are isotropic in orientation. As a consequence, we do not expect exact agreement; however, we do anticipate equivalent scaling with $\lambda$. Indeed, the maximum pressure, shown in Fig. 3(C), which is always located at the center of the colony, shows great agreement with the scaling of the analytical solution.

Figure 3(B) shows the colony radius $R(t)$ over time for discrete (solid) and continuum (dashed) simulations for various $\lambda$. At early times the growth rates are comparable and show an exponential trend. The curves begin to diverge around $t\approx 5$, with larger values of $\lambda$ corresponding to slower growth. The continuum theory is in excellent agreement with the discrete model, showing that, while the model may differ at instantaneous points in time, as in Fig. 3(A), it provides an accurate description of the system’s dynamical evolution. From the continuum theory, for $\lambda=0$ we have $R(t)=R_{0}e^{t/2}$. For any non-zero value of $\lambda$, however, taking $\lambda^{1/2}R\gg 0$ we instead find linear asymptotic scaling with time ${R(t)\sim(2/\lambda)^{1/2}t}$. Figure 3(D) shows the rate of change $R^{\prime}(t_{f})$ at the final time, determined by a linear fit. At larger values of $\lambda$ we see an approach to the asymptotic slope, though it does not appear to be in the asymptotic regime for the lower values shown here.

To characterize the emergent waves, we consider the radially averaged length distribution $\langle\ell_{n}\rangle$ from the particle simulations, shown in Fig. 4(A). For all non-zero sensitivities a wave-like signal is identified. The wavelength appears to decrease with the distance from the center of the colony $r$, and generally has a stronger signal at large distances. A quantitative characterization is given by the wavelength between the last two peaks, whose evolution in both discrete and continuum simulations is shown in Fig. 4(B). The time of appearance of the first wave tends to decrease with $\lambda$, as indicated by the lack of data at early times. Once a well-defined wave has formed, its wavelength decreases over time, approaching an asymptotic value which scales with $\lambda^{-1/2}$ as shown in Fig. 4(C). This scaling can be understood through the solution for the average length (4). Rescaling space as $\tilde{r}=\lambda^{1/2}r$ makes the growth rate $g(\tilde{r})$ independent of $\lambda$. Because the colony grows linearly with rate $\propto\lambda^{-1/2}$ at late times, this suggests the integral in Eq. (4) approaches the same value for characteristics $\hat{r}_{\lambda}(t)$ with initial position $\hat{r}_{\lambda}(0)\propto\lambda^{1/2}$.

This work shows through large-scale particle simulations and a corresponding continuum theory how stress-sensitive growth affects the spatiotemporal dynamics of cell size in proliferating cell collectives. Concentric ring morphologies like those reported here are common and often form through chemical means and/or periodic processes Shimada et al. (2004); Liu et al. (2011, 2021). In contrast, the patterns we observe arise from a constant process of mechanical resistance. In natural settings, such effects would be compounded by other growth-limiting factors such as access to nutrients, and controlled measurements may be needed to assess their relative importance to specific organisms.

The continuum theory is generic and its predictions might apply to other cell collectives such as tissues, though material elasticity is likely an important factor that we do not address here Goriely (2017). Nonetheless, similar patterns have been observed in simulations of a vertex model of growing epithelial layers Carpenter et al. (2024). A prediction of our model is that colonies of the same size can reduce internal pressure through stress-sensitivity, and organisms might exploit this property to optimize the balance of growth and stress. The predicted linear growth of the colony for non-zero sensitivities is consistent with experimental observations Hallatschek et al. (2007) and can be used to estimate the model parameter $\lambda$. Using data from Hallatschek et al. (2007), for E. Coli the growth timescale is $\tau\approx 54{\rm min}$, the initial cell length is approximately $\ell_{0}\approx 1\mu{\rm m}$, and the colony radial velocity is $dR/dt\approx 1.3\times 10^{-2}{\mu\rm m}/{\rm s}$, yielding the dimensionless sensitivity $\lambda\approx 1.1\times 10^{-3}$, which is in the range of values where the mechanical effects discussed here are relevant.

The Smoluchowski equation provides an explicit coupling between the micro- and macro-scales, and its solution in terms of characteristics provides an analytical description of inheritance, or how the cumulative stress of a cell line affects the size and growth of subsequent generations. Within this context, other inherited variables, such as genes and mutations thereof, could be incorporated into the model and coupled to growth. Finally, while not discussed here, the discrete simulations exhibit rich orientational structure, and future work will explore the dependence of this structure on growth sensitivity.