A competitive advantage through fast dead matter elimination in confined cellular aggregates

Collective competition is ubiquitous in proliferating cellular systems, with implications ranging from maintenance of stable homeostasis in tissues [1], to robust morphogenesis [2], tumor development and suppression [3, 4] and survival of bacterial communities [5, 6], for example during range expansion [7]. Competitive advantages can be achieved either via direct adversarial interactions such as induction of elimination [8–10] (i.e. one cell inducing death in another) or due to competitive exclusion [11–15] where sharing a common resource or limited space leads to competition even without any specific response depending on the cell type of rival cells. The current work is aimed at the latter, non-adversarial, competition scenario.

While many past works focused on species' typical division and death rates, arguing that fitness increases as the difference between the former and the latter [2], this intuitive argument cannot be used when cell aggregates do not expand but instead compete for free space in confinement, as validated numerically e.g. in [16]. In this case, division is, on average, balanced by death, creating a dynamic steady state known as homeostasis [4, 17, 18]. Then, the question arises whether an emergent property of the homeostatic state itself (in which competing cell types are absent) can instead serve as a suitable fitness indicator for competition outcomes. We ignore the role of chemical signalling, which can also participate in the competition that controls homeostasis [19]. In systems with pressure-regulated growth [20] and an equation of state relating pressure to the density of (living) cells, one property that has already been proposed is the homeostatic pressure [4, 16]. This makes intuitive sense, since, in a mixture of cell types with different homeostatic pressures, there is a range of intermediate pressures at which growth can still exceed death in one species while death exceeds growth in the other. However, an equation of state implicitly assumes that the active (growing) cells are the only mechanically significant component in the system. In contrast, the inclusion of any other continuously produced mechanical components such as dead cells or extracellular matrix may decouple the density of growing (active) cells from pressure and lead to a fitness that is not captured by homeostatic pressure.

Here, we demonstrate that the presence of passive by-products of the processes maintaining homeostasis can indeed significantly affect a species' ability to compete. Specifically, we focus on systems where these by-products are cells that have already been fated to death and lost their ability to grow and divide, but still retain full or partial structural integrity before they are eventually removed. In real biological systems, the time scale for the removal of dead or dying cells can have a wide range of values relative to the cell's life cycle depending on the mechanism of removal. Prominent examples are extrusion [17, 21, 22], structural disintegration [23], phagocytosis [9] or 'explosive' bacterial lysis [24]. Until now, models of growth-driven dense cellular matter have often implicitly employed instantaneous cell removal processes both in continuum theories and in most simulation works [25, 26] (and see also the discussion section). Taking into account this additional time scale is therefore not only theoretically interesting, but also of immediate biological relevance. However, we would like to note that remnants of dead cells or cells destined for apoptosis do not necessarily sterically hinder growth of other cells in all biological systems. In fact, they can lead to active force generation which can be important for development [27] or extrusion [28], something we do not consider here.

We consider a prototypical cellular system with only a few essential 'building blocks' (figure 1), namely steric repulsion, cell growth and proliferation, cell death, and, for the sake of cell number homeostasis, also pressure regulation of growth [29–31]. A final ingredient of our model, which is central to our investigation, is the explicit, finite-time removal of dead cells, which translates to a continuous supply of (space-filling) dead passive matter and thus breaks the equation of state between pressure and the density of active, growing cells (figures 2(a)–(d)). In practice, the chosen removal mechanism consists of a decaying steric interaction for dead cells, as described in greater detail below. The relevant intuitive picture is the degradation of dead cells' structural integrity and load supporting capability and a resultant reduction of effective volume.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Supplementary video 1: comparison of dynamics for the two species that were presented in figure 2 (top: E_H = 4, ${k}_{\mathrm{H}}^{\mathrm{d}\mathrm{e}\mathrm{g}}=18.60{k}^{\mathrm{d}\mathrm{e}\mathrm{a}\mathrm{t}\mathrm{h}}$; bottom: E_L = −20, ${k}_{\mathrm{L}}^{\mathrm{d}\mathrm{e}\mathrm{g}}=10.05{k}^{\mathrm{d}\mathrm{e}\mathrm{a}\mathrm{t}\mathrm{h}}$) in separate single-species samples at homeostasis (a subsection of around 3% of the system length is shown for 200 time units, or around 3 generations on average). Dead degrading cells are colored black and denoted by a dashed outline and increased transparency denotes the current softening factor S_i .

Supplementary video 2: dynamics in a competition assay between the same species as in supplementary video 1 for 1000 time units (∼14 generations) at the very start of the simulation and displayed at ×10 speed compared to supplementary video 1.

In such a system with by-products, we see that a species which removes dead matter faster and thus has a higher proportion of active, growing, cells can outcompete a species with a smaller active density even when mechanical pressure would suggest spatial invasion in the reverse direction (figures 2(c)–(f)). This initial evidence raises the prospect of a complementary competition paradigm, which we explore in this study using both theoretical considerations and systematic confirmation in simulations. We start with a simple mean-field continuum theory and show why this level of description does not reveal any particular advantage for species with faster dead matter degradation. However—keeping in mind that a well-mixed state, albeit an adequate description in a coarse-grained continuum model, is inevitably broken at the scale of finite cell granularity—we then provide a theoretical argument for how a fitness advantage could be induced by discrete cellular units, spatially heterogeneous growth rates and a finite mechanical screening length, and confirm that a compatible spatial heterogeneity emerges in our simulations. Finally, we examine systematically in simulations how the active density is correlated to fitness and compare it to the fitness predicted by homeostatic pressure. We show that a homeostatic property, measured in a homogeneous cell sample, can determine fitness in competition and that competition is driven by interfaces between species.

2.1. Theoretical considerations

Before we present a detailed numerical analysis of our system, we first explore the consequences of passive matter production theoretically, in order to narrow down the possible mechanisms that can cause effects such as those seen in figure 2. The most general description of growing matter that adapts its growth/death dynamics to its local mechanical environment is on the level of continuity equations, which we first discuss in the context of a single species:

$$\begin{equation}\frac{\partial {\rho }^{\mathrm{a}}}{\partial t}+\nabla \cdot ({\rho }^{\mathrm{a}}\mathbf{v})=({k}^{\mathrm{d}\mathrm{i}\mathrm{v}}(P)-{k}^{\mathrm{d}\mathrm{e}\mathrm{a}\mathrm{t}\mathrm{h}}){\rho }^{\mathrm{a}}\end{equation} \tag{ 1 }$$

$$\begin{equation}\frac{\partial {\rho }^{\mathrm{p}}}{\partial t}+\nabla \cdot ({\rho }^{\mathrm{p}}\mathbf{v})={k}^{\mathrm{d}\mathrm{e}\mathrm{a}\mathrm{t}\mathrm{h}}{\rho }^{\mathrm{a}}-{k}^{\mathrm{deg}}{\rho }^{\mathrm{p}}.\end{equation} \tag{ 2 }$$

Here, the superscripts a and p indicate active cells (capable of division) and passive (dead) cells, respectively. We assume that all cellular matter has identical mechanical properties, where the isotropic local pressure is determined by an equation of state $P=\mathcal{F}({\rho }^{t})$ with the total density ρ^t = ρ^a + ρ^p (for simplicity, we also assume that $\mathcal{F}$ is an invertible function in the relevant range). In the over-damped limit, the velocity of the otherwise immotile matter is caused solely by intercellular forces according to force balance with a drag force ζρ^t v = −∇P, ζ being the drag coefficient. Therefore, both P and v are identical between equations (1) and (2). Active cells constantly grow but reduce their growth rate in response to pressure, with a monotonically decreasing function for the growth rate k^div(P). However, when cells die with rate k^death, they do not vanish immediately from the system but are first converted to passive matter, which itself degrades with a rate k^deg. These last two rates, that could in principle also depend on pressure, are taken as constant for simplicity. In a system of finite size, active density together with its passive remnants fills the system, growing with a local exponential rate k^div(P) − k^death, until growth is inhibited sufficiently by pressure.

Considering a mean-field approximation, where we assume that ρ^t and P are spatially uniform and thus v = 0, the steady state is determined by local homeostasis, which is defined by equation (1) as k^div(P^h) − k^death = 0, translating to a mean-field homeostatic total density of ${\tilde{\rho }}^{\mathrm{t}}={\mathcal{F}}^{-1}({P}^{\mathrm{h}})$. Meanwhile, equation (2) requires for homeostasis that ${k}^{\text{death}}{\tilde{\rho }}^{\mathrm{a}}={k}^{\mathrm{deg}}{\tilde{\rho }}^{\mathrm{p}}$, such that in steady state

$$\begin{equation}{\tilde{\rho }}^{\mathrm{p}}=\frac{{k}^{\mathrm{d}\mathrm{e}\mathrm{a}\mathrm{t}\mathrm{h}}}{{k}^{\mathrm{deg}}}{\tilde{\rho }}^{\mathrm{a}}\end{equation} \tag{ 3 }$$

$$\begin{equation}{\rho }^{\mathrm{h}}={\tilde{\rho }}^{\mathrm{a}}=\frac{{\tilde{\rho }}^{\mathrm{t}}}{1+\frac{{k}^{\mathrm{d}\mathrm{e}\mathrm{a}\mathrm{t}\mathrm{h}}}{{k}^{\mathrm{deg}}}},\end{equation} \tag{ 4 }$$

where we have defined the homeostatic active density ρ^h as a shorthand for the steady-state of the active density ρ^a, as we will refer to it throughout this study. Note that, even if the mean-field quantities vary in time, ${\tilde{\rho }}^{\mathrm{p}}$ is still entirely determined by the history of the active density and approaches equation (3) on a time scale 1/k^deg if ${\tilde{\rho }}^{\mathrm{a}}$ is constant. Therefore, we can view equation (4) also as a quasi-steady-state relation between ρ^t and ρ^a when ρ^a changes very slowly.

We now turn to a competition scenario with two different species. From the mean-field consideration above, we see that k^deg, the rate at which the passive matter is degraded, determines the relative proportions of active and passive density in homeostasis. In contrast, the homeostatic pressure P^h and the homeostatic total density ${\tilde{\rho }}^{\mathrm{t}}$ are solely determined by the condition k^div(P^h) − k^death = 0, independent of the amount of passive matter at this level of description. If we were to consider two species with different homeostatic pressures (e.g., through different k^div(P)), a situation similar to that of [4, 16, 32] arises, where the species with higher homeostatic pressure has a clear growth advantage due to being able to exhibit net growth at higher pressures. In contrast, we consider here two species with dynamics as in equations (1) and (2) that differ only in their respective passive matter degradation rates ${k}_{\mathrm{H}}^{\mathrm{deg}} > {k}_{\mathrm{L}}^{\mathrm{deg}}$, with ${k}_{\mathrm{H}}^{\mathrm{deg}}$ yielding a high active density ('H') and ${k}_{\mathrm{L}}^{\mathrm{deg}}$ yielding a low active density ('L') in homeostasis according to equation (4) ³ .

Two species only differing in their passive matter degradation rates ${k}_{\mathrm{H}}^{\mathrm{deg}} > {k}_{\mathrm{L}}^{\mathrm{deg}}$ are characterized by two copies of equations (1) and (2) for the densities ${\rho }_{\mathrm{H}}^{\mathrm{a}}$, ${\rho }_{\mathrm{H}}^{\mathrm{p}}$, ${\rho }_{\mathrm{L}}^{\mathrm{a}}$ and ${\rho }_{\mathrm{L}}^{\mathrm{p}}$. Again, v and P are both identical between all four equations. This leads to tight constraints on the possible relative behavior of the two species. For example, the mean-field version of equation (1), ${\dot{\tilde{\rho }}}_{\mathrm{H},\mathrm{L}}^{\mathrm{a}}=({k}^{\mathrm{d}\mathrm{i}\mathrm{v}}(P)-{k}^{\mathrm{d}\mathrm{e}\mathrm{a}\mathrm{t}\mathrm{h}}){\tilde{\rho }}_{\mathrm{H},\mathrm{L}}^{\mathrm{a}}$, is identical for both the H and the L species, implying

$$\begin{equation}\frac{\mathrm{d}}{\mathrm{d}t}\left(\frac{{\tilde{\rho }}_{\mathrm{H}}^{\mathrm{a}}}{{\tilde{\rho }}_{\mathrm{H}}^{\mathrm{a}}+{\tilde{\rho }}_{\mathrm{L}}^{\mathrm{a}}}\right)=0,\end{equation} \tag{ 5 }$$

which excludes any meaningful competition, regardless of the behavior of the passive matter ${\rho }_{\mathrm{H}}^{\mathrm{p}}$, ${\rho }_{\mathrm{L}}^{\mathrm{p}}$. As long as the mean-field assumption of spatially uniform pressure applies, a similar equation holds for the total cell numbers, even if ${\rho }_{\mathrm{H}}^{\mathrm{a}}$, ${\rho }_{\mathrm{H}}^{\mathrm{a}}$ are spatially varying. By integrating equations (1) and (2), we have

$$\begin{equation}{\dot{N}}_{\mathrm{H},\mathrm{L}}^{\mathrm{a}}=({k}^{\mathrm{d}\mathrm{i}\mathrm{v}}(P)-{k}^{\mathrm{d}\mathrm{e}\mathrm{a}\mathrm{t}\mathrm{h}}){N}_{\mathrm{H},\mathrm{L}}^{\mathrm{a}}\end{equation} \tag{ 6 }$$

for both species, where ${N}_{\square }^{\mathrm{a}}=\int {\rho }_{\square }\,\mathrm{d}\mathbf{r}$ (with □ = L or H) indicates particle number corresponding to the relevant densities ${\rho }_{\mathrm{H}}^{\mathrm{a}}$ and ${\rho }_{\mathrm{L}}^{\mathrm{a}}$, again implying the same relative change for both.

When the different species with their active and passive matter components are well-mixed on a microscopic scale, it is therefore not obvious how a species that leaves less passive matter behind could have a particular competitive advantage such as that seen in figure 2. Here, we argue that it is the granularity of the system at the cellular scale and its non-equilibrium nature that leads to a competitive advantage for the H species: due to the growth and removal of discrete cells and finite mobility, the system never reaches a true steady-state, but is actively kept away from mechanical equilibrium (by which we mean the static configuration it would converge to if all growth/degradation processes where halted). It is easy to see that systems such as equation (1) have a well-defined length scale, or screening length, ${l}_{s}=1/\sqrt{-\zeta {\alpha }^{\prime }{\rho }^{\mathrm{h}}}$ for pressure perturbations, where α' = dk^div/dP is the growth response of the system at homeostasis (see appendix C). This screening length is the result of a competition between the diffusive spread of mechanical perturbations (which depends on the drag coefficient ζ) and the recovery of the cellular medium via adaptive growth/removal of matter which seeks to locally bring the system back to the homeostatic pressure P^h. If l_s is comparable to the cell size, a local pressure drop (for example, the removal of a cell) causes only a small number of discrete cells in its vicinity to grow in order to replenish the cellular material. At the interface between the two species, this finite-size selection can potentially lead to small-number effects. Indeed, a lattice-based thought experiment incorporating our most important assumptions—identical growth response for active cells of both species and different equilibrium fractions of passive matter—shows that there is a statistical bias towards finding fewer living cells of species L in this region than necessary to compensate cell removal (see appendix D), leading to a long-term takeover of species H. Therefore, species with a higher proportion of active, growing matter are statistically better at exploiting opportunities to 'conquer' new ground. This bias can also be shown to vanish for increasing screening length. Consequently, a finite screening length in conjunction with the discreteness of the system can yield a growth advantage for species H that is proportional to the difference in active densities ${\Delta}{\rho }^{\mathrm{h}}={\rho }_{\mathrm{H}}^{\mathrm{h}}-{\rho }_{\mathrm{L}}^{\mathrm{h}}$.

We therefore expect that two species only differing in their passive matter degradation rate k^deg should have similar pressure-controlled growth rates during homeostasis when growing separately (because of identical homeostatic pressure P^h). The same should be true for cells of both species when sharing space in the same system given they are far enough away from interfaces (because of a finite mechanical screening length). In contrast, based on the thought experiment mentioned above (see appendix D for details), growth rates near interfaces between the two species should be biased in favor of the species with a higher homeostatic active density ρ^h.

Another prediction based on the above considerations is that, with decreasing drag coefficient ζ, the competitive advantage conferred by a higher active density should diminish, as this increases screening length ${l}_{s}=1/\sqrt{-\zeta {\alpha }^{\prime }{\rho }^{\mathrm{h}}}$. In fact, this is quite intuitive since mechanical relaxation in the system becomes instantaneous in the limit ζ → 0, such that all cells experience the same pressure at all times, leading to identical growth rates, similar to the mean-field scenarios considered earlier.

2.2. Numerical experiments

To test the correlation between homeostatic active density and fitness (as well as the other predictions outlined above), we employ a minimal mathematical model which describes the dynamics of discrete cellular units and incorporates the above-mentioned building blocks (figure 1). In line with our theoretical considerations, we assume identical mechanical properties for all species and only a single differentiating feature: the degradation rate of passive cells. We then run single-species simulations to extract homeostatic properties and link them to the outcome of binary competition simulations.

2.2.1. Agent-based numerical model

First, we briefly describe the cellular model used in the simulations and specifically the explicit cell removal process. The full details of the model can be found in appendix A. We describe a 'dry' cellular system where individual cells are modeled as immotile dumbbells, each composed of two nodes of radius R. Length is henceforth measured in units of node diameter. Nodes of different cells have a steric Hertzian [33] repulsion between them ${f}_{{i}^{\alpha }{j}^{\beta }}\sim \delta {r}_{{i}^{\alpha }{j}^{\beta }}^{3/2}$ as a function of the node overlap $\delta {r}_{{i}^{\alpha }{j}^{\beta }}$ (but only between the nearest nodes of the two cells). Here, Latin indices i, j denote cells, while Greek indices α, β distinguish between nodes of a cell. Nodes of the same cell have a non-linear spring interaction between them (also Hertzian), with rest length ${\ell }_{i}^{\mathrm{r}\mathrm{e}\mathrm{s}\mathrm{t}}$. Cell growth is implemented via elongation of the internal spring's rest length ${\ell }_{i}^{\mathrm{r}\mathrm{e}\mathrm{s}\mathrm{t}}$, expressed equivalently via a growth phase ${g}_{i}^{\mathrm{g}\mathrm{r}\mathrm{o}\mathrm{w}\mathrm{t}\mathrm{h}}={\ell }_{i}^{\mathrm{r}\mathrm{e}\mathrm{s}\mathrm{t}}/2R$ in the range [0, 1]. A cell starts its life with two nodes fully overlapping $({g}_{i}^{\mathrm{g}\mathrm{r}\mathrm{o}\mathrm{w}\mathrm{t}\mathrm{h}}=0)$ and, if it does not die before then, ends it with division when ${g}_{i}^{\mathrm{g}\mathrm{r}\mathrm{o}\mathrm{w}\mathrm{t}\mathrm{h}}$ reaches 1 and the two nodes just touching (zero overlap), as depicted in figure 1.

Each cell has a base growth rate ${\gamma }_{i}=1/{T}_{i}^{\mathrm{d}\mathrm{i}\mathrm{v}}$, where ${T}_{i}^{\mathrm{d}\mathrm{i}\mathrm{v}}$ corresponds to the division time in the absence of external forces. We choose this basal growth to be significantly faster than the death rate k^death (see further below) which cells have to approximately match in homeostasis by decreasing their growth rate. For this purpose, the base growth rate is modulated by a monotonously decreasing sigmoidal function of the cell-sensed pressure $G({P}_{i})=\frac{1}{2}\left(1-\mathrm{e}\mathrm{r}\mathrm{f}\left(\lambda \left(-1+2{P}_{i}/{P}^{\mathrm{s}\mathrm{t}\mathrm{a}\mathrm{l}\mathrm{l}}\right)\right)\right)$ (see appendix A for details), making the overall growth rate of ${\dot{g}}_{i}^{\mathrm{g}\mathrm{r}\mathrm{o}\mathrm{w}\mathrm{t}\mathrm{h}}={\gamma }_{i}\cdot G({P}_{i})$ dependent on external forces and thereby allowing for homeostasis. λ determines how steeply G(P_i ) decays and P^stall is the pressure at which cell growth effectively halts. A visualization of G(P_i ) is provided in appendix A (figure A1). Upon division, each node of a mother cell turns into a new daughter cell (${g}_{{d}_{1}}^{\mathrm{g}\mathrm{r}\mathrm{o}\mathrm{w}\mathrm{t}\mathrm{h}}$ = 0, ${g}_{{d}_{2}}^{\mathrm{g}\mathrm{r}\mathrm{o}\mathrm{w}\mathrm{t}\mathrm{h}}=0$), itself composed of two connected fully overlapping nodes. To avoid synchronization of divisions due to lineage extinctions, the new base growth rates ${\gamma }_{{d}_{1}},{\gamma }_{{d}_{2}}$ are randomly drawn from a normal distribution (as shown e.g. in bacteria in [34, 35]).

Death of a cell is stochastically initiated for each cell at a rate k^death. This fixed rate also sets the simulation time scale. In particular, the growth rate of cells γ_i ⋅ G(P_i ) decreases to approximately match this rate during homeostasis and the cell doubling time in our model is thus ln(2)/k^death. In this work we use the value k^death = 10⁻² that is much lower than the average base growth rate and the cell degradation rate (see below). The corresponding doubling time is ∼69, or around 14 generations (gen) on average per 1000 time units (see appendix E for details on the expected homeostatic pressure underlying this resulting growth rate). The activity of the cells, and any stochasticity, is present only in their growth, and their division and death events. Growth and division, coupled to the inter-cell steric repulsion, result in local forces which can be relaxed to a certain extent through cell motion (see over-damped equations of motion in appendix A) but which, in a closed system, lead to an overall pressure buildup. This in turn feeds back on growth which leads to a dynamic state of homeostatic balance between cell division and death.

To the quite standard model described so far, we now add a new mechanism of mechanical cell degradation and removal (see figure 1). Analogously to the continuum consideration above, the idea is that cells temporarily retain their mechanical integrity when they die, but subsequently fade from the system on an explicit and tunable time scale. In practice, this is achieved by a reduction of the inter-cell steric repulsion via 'node softening'. For this purpose, each cell is assigned a factor S_i which scales the default Hertzian inter-cell forces: ${f}_{{i}^{\alpha }{j}^{\beta }}^{\mathrm{i}\mathrm{n}\mathrm{t}\mathrm{e}\mathrm{r}-\mathrm{c}\mathrm{e}\mathrm{l}\mathrm{l}}{\Rightarrow}{S}_{i}{f}_{{i}^{\alpha }{j}^{\beta }}^{\mathrm{i}\mathrm{n}\mathrm{t}\mathrm{e}\mathrm{r}-\mathrm{c}\mathrm{e}\mathrm{l}\mathrm{l}}{S}_{j}$ (see equation (A.5) in appendix A). For a living growing cell, S_i = 1. However, once death is initiated, growth stops irreversibly so that the cell turns from active to passive matter. The scaling factor S_i then decays according to ${S}_{i}=1-(1-{\mathrm{e}}^{-E\cdot {g}_{i}^{\mathrm{deg}}})/(1-{\mathrm{e}}^{-E})$ with progress tracked by a degradation phase ${g}_{i}^{\mathrm{deg}}=(t-{t}_{i}^{\mathrm{d}\mathrm{e}\mathrm{a}\mathrm{t}\mathrm{h}})/{T}^{\mathrm{r}\mathrm{e}\mathrm{m}\mathrm{o}\mathrm{v}\mathrm{a}\mathrm{l}}$ that grows linearly in time from 0 to 1 over a time period T^removal, starting at the time of death (see visualization in figure 3(a)). The parameter E is used to control the softening profile and is shown below and in figure 3(b) to control the rate of the mechanical degradation of the cell k^deg. When S_i approaches zero, the cell effectively can no longer affect its surroundings (and when it fully decays to zero, the cell is simply removed from the system). In our system, we chose the maximal duration of cell degradation to be T^removal = 10 (the time for ${g}_{i}^{\mathrm{deg}}$ to change from 0 to 1), so that cell degradation happens on the order of our model time units. The degradation time can thus be either very short or can be significant, but is in any case substantially shorter than the typical doubling time of ∼69. The full set of parameters can be found in appendix A.

Zoom In Zoom Out Reset image size

While the above numerical model can be easily generalized to 2D and 3D, here, the long-axis orientation of the cells, their growth and their motion are all constrained to a continuous 1D periodic domain (see figure 1). We use a channel length of size 2560 in node diameter units, which is several orders of magnitude larger than the screening length (cf appendix C), and large enough to eliminate finite size effects. The initial seeding is always done such that cells are not initially in contact (see full details of initial conditions in appendix A).

2.3. Single species homeostatic properties

To measure homeostatic properties such as ρ^h and P^h for a given cell species, we run simulations of the above model, each with a single species. Each such simulation is run for a time ${t}_{\mathrm{h}}^{\mathrm{s}\mathrm{i}\mathrm{m}}$ of 2000 (∼29 gen). Time-averaged measurements, such as for the homeostatic properties ρ^h, P^h, are performed in a time window 1000 < t < 2000, well after the steady-state is reached (see figure 2). In the following, we use the brackets ⟨ ⟩, to denote an average over an ensemble of realizations. As in our theoretical considerations, the species we simulate differ only in the degradation behavior of passive matter, which is controlled by the softening factor E in our agent-based model (see previous section) as shown in figure 3(a). To calculate an effective time scale of cell degradation from this central parameter E, we use the time it takes the force scaling factor S_i to drop by 90% after death, representing the fact that despite the fixed T^removal, cells are effectively removed already when their S_i value decays close enough to zero. This allows us to explicitly calculate the effective degradation rate as

$$\begin{equation}{k}^{\mathrm{deg}}=\frac{1}{t{\vert }_{{S}_{i}=\,0.1}-t{\vert }_{{S}_{i}=1}}=-\frac{E}{{T}^{\,\mathrm{r}\mathrm{e}\mathrm{m}\mathrm{o}\mathrm{v}\mathrm{a}\mathrm{l}}\,\mathrm{ln}\left[1-0.9(1-{\mathrm{e}}^{-E})\right]}.\end{equation} \tag{ 7 }$$

These effective degradation rates k^deg are shown in figure 3(b) for a range of E values. Since k^deg is a monotonically increasing function of E, species H and L with E_H > E_L also have ${k}_{\mathrm{H}}^{\mathrm{deg}} > {k}_{\mathrm{L}}^{\mathrm{deg}}$, which is consistent with the considerations above and figure 2. Note that, besides the overall time scale 1/k^deg on which cells mechanically persist after death, the parameter E also changes the 'collapse' time, i.e. the abruptness with which S_i approaches zero, a process which starts immediately after death for E → ∞ and only shortly before ${g}_{i}^{\mathrm{deg}}$ reaches 1 for E → −∞ (see figure 3(a)). This does not influence the competition mechanism discussed here directly, but it changes the pressure fluctuations in a non-monotonic way, as is explained below and in detail in appendix E.

An intuition for the general behavior of our system can already be gained from the representative examples shown in figures 2(a)–(d) (supplementary video 1) in light of the theoretical considerations presented above. Figures 2(a) and (b) visualize coarse-grained spatio-temporal dynamics for two separate cell populations by tracking cell positions through time (only the first 50 time units out of 2000 are shown). Colored regions represent living cells and any black spaces in between denote either voids or degrading (dead) cells. Figure 2(a) corresponds to a species H with E = 4 and thus a faster degradation rate ${k}_{\mathrm{H}}^{\mathrm{deg}}=18.60{k}^{\mathrm{d}\mathrm{e}\mathrm{a}\mathrm{t}\mathrm{h}}$ compared to E = −20 and ${k}_{\mathrm{L}}^{\mathrm{deg}}=10.05{k}^{\mathrm{d}\mathrm{e}\mathrm{a}\mathrm{t}\mathrm{h}}$ for species L in figure 2(b). The seeded cells fill up the available volume already within 5 time units. The higher density of black lines that is apparent for species L in the steady state represents a greater accumulation of passive matter and supports our expectation from equation (4) that the homeostatic active density for species L is lower than for species H. Indeed, the time series for the active density in figure 2(c), measured as the number of growing cells divided by system length (averaged over 10 realizations, only the first 500 time units are shown), clearly confirms this as does the homeostatic value (denoted in the figure).

It is also useful to consider the pressure, since it is the mechanical property that controls growth. While, a priori, we do not expect an inherent difference in homeostatic pressure (see theoretical considerations above), the time series in figure 2(d) indeed shows a small but significant difference between the two species. The pressure is measured via the Irving–Kirkwood expression [36] which in 1D takes the form: $P(t)=\frac{1}{2}{\sum }_{i > j}\,{f}_{{i}^{\alpha }{j}^{\beta }}{r}_{{i}^{\alpha }{j}^{\beta }}$ (averaged over 10 realizations in the figure) with α and β denoting nearest nodes of the neighbor cells ⁴ . The homeostatic value P^h in this specific case is slightly higher for species L compared to H as denoted in the figure.

We next measure these two homeostatic properties for species with a range of E values and associated degradation rates. One can immediately see in figure 4 that they show systematically different behavior across the parameter space: the homeostatic active density shows a monotonous increase with E and is predicted exceptionally well by the mean-field expression in equation (4). ⁵ The homeostatic pressure, on the other hand, is not constant as the mean-field consideration predicts (although here it increases by only $\sim 2\%$ of its minimal value). The reason for this are local pressure fluctuations generated by degrading dead cells (see appendix E). This is most evident when considering the E → ±∞ limits where a cell vanishes instantaneously and the pressure drops abruptly. Here, the measurements performed for E < −30 and E > 30 correspond to a regime in which the mechanical collapse of the cell is indeed much faster than the recovery response of the surrounding cell medium (S_i approaches a step function; see figure 3(a)). One should note that the fast-cell-collapse limit is not a peculiarity of the current model but in fact common practice in studies of competition and cellular active matter in general which use it implicitly by removing cells instantaneously [4, 16]. At the other extreme of the spectrum, the pressure should be most closely approximated by the mean-field prediction when the local pressure fluctuations are negligible (compared to the mean value). For the numerical model employed here, the global pressure fluctuations are minimal at E ≃ 4 (see appendix E), where indeed also the pressure is minimal: min(P^h) = 1.924, and away from this point increased fluctuations lead to increased homeostatic pressure. The zeroth order mean-field value for the homeostatic pressure is ${\tilde{P}}^{\mathrm{h}}=1.919$ (dashed red line in figure 4(b); see appendix E for derivation) which is even lower than the minimum measured in simulation, compatible with the fact that the local pressure fluctuations are never truly negligible in this numerical model. The observed deviation of P^h from the zeroth order mean-field calculation is in fact useful and will allow us to extract an alternative fitness prediction from the homeostatic pressure [4, 16] that can be compared to that of the active density.

Zoom In Zoom Out Reset image size

2.4. Comparison of competition drivers for cell mixtures

To determine how these homeostatic properties are correlated with species fitness, we perform pairwise competition assays between the species differing in their softening factors E and associated passive matter degradation rates, and observe the relative change in species abundance over time. In each simulation, the channel is initially seeded with a balanced mixture of both species in a random arrangement.

A representative example of such a competition assay, which pits a species L against a species H with a higher degradation rate, is shown in figures 2(e)–(f) (supplementary video 2). From the spatio-temporal dynamics of the cells shown (coarse-grained) in figure 2(e) we see that, besides general coarsening dynamics, patches of both species H and L grow and shrink stochastically. Yet, species H has a clear statistical advantage and increases its proportion in the sample at the expense of species L, as quantified by the species fractions N_H(t)/N(t) and N_L(t)/N(t) in figure 2(f). The observation that a species H with a higher homeostatic active density wins the competition is compatible with the prediction of the lattice-based thought experiment discussed above (and in appendix D).

In the competition assays for other E value combinations, we continue to label cells with a higher E value as species H and those with lower E as species L. We focus on the parameter regime in which we expect the clearest distinction between the two predictions; namely −20 ⩽ E ⩽ 20 (denoted by green shading in figure 4) where both homeostatic pressure and active density change significantly and where the functional form of S_i is more sensitive to changes in E (see figure 3). For reference, we also run simulations where all cells have the same E value and arbitrarily label two equal size subsets of the seeded cells as H and L. Measuring 'fitness' for such fully neutral competition allows us to estimate the natural statistical variability that arises from the model and from the chosen protocol of fitness measurement. We define relative fitness for species H as the fraction of the final proportion of cells ${N}_{\mathrm{H}}^{\mathrm{f}}$ in the entire population N^f, measured after a fixed time ${t}_{\mathrm{c}\mathrm{o}\mathrm{m}\mathrm{p}\mathrm{e}\mathrm{t}\mathrm{e}}^{\mathrm{s}\mathrm{i}\mathrm{m}}=5000$ (∼72 gen). This simulation time was selected to be long enough to show a clear trend for the change of species fraction (cf figure 2(f)), but not long enough to reach the absorbing extinction boundary (as discussed in appendix F), which would have biased the measurement [37]. The relative fitness measure ${N}_{\mathrm{H}}^{\mathrm{f}}/{N}^{\mathrm{f}}$ averaged over an ensemble of simulation realizations is shown in figures 5(a) and (b). As expected, the species fraction in neutral competition (equal E values) is 0.502 (±0.005). At the other extreme, the relative fitness in this numerical experiment grows up to a maximum of 0.75 (±0.02) showing that some species carry a distinct advantage over others.

Zoom In Zoom Out Reset image size

Across the parameter space, the species labeled H exhibits either a similar or a higher fitness than species L, indicated by $\langle {N}_{\mathrm{H}}^{\mathrm{f}}/{N}^{\mathrm{f}}\rangle \,\gtrsim \,0.5$ (figures 5(a) and (b)). It is also worth noting that, in the range −20 ⩽ E ⩽ −5, there is no clear advantage for any of the species, although the homeostatic pressure changes considerably in this range (figure 4(b)). These two observations already indicate a dominant role of the homeostatic active density.

For a more quantitative analysis, we correlated the final cell fraction with the difference in homeostatic active density between the competing species. We do the same for the homeostatic pressure and in both cases normalize by a reference value at E = −60. The results of these two measurements are shown in figures 5(c) and (d). For the active density (figure 5(c)), all the data points fall in the quadrants of positive correlation (yellow shading) within statistical variability. A clear linear trend indicates that the fitness as defined by the species fraction can accurately be predicted from the difference in homeostatic density. From a different perspective, a larger homeostatic active density difference leads to a faster average rate at which the winner species takes over the sample, in line with our theoretical framework. The linear dependence can be fitted as

$$\begin{equation}\langle {N}_{\mathrm{H}}^{\mathrm{f}}/{N}^{\mathrm{f}}\rangle =0.50+2.95({\Delta}\langle {\rho }^{\mathrm{h}}\rangle /\langle {\rho }_{E=-60}^{\mathrm{h}}\rangle ).\end{equation} \tag{ 8 }$$

The root mean square error (rmse) of the fit is 0.02 which matches the scatter of data points for identical species (i.e. at Δρ^h = 0).

For the homeostatic pressure (figure 5(d)), an equally unambiguous correlation is not observed, even taking into account statistical variability (greyed out areas in the figure). For certain mixtures, a species with a lower homeostatic pressure can in fact win the competition (upper left quadrant), as already seen in figure 2. Note that 'winning' here refers to the consistent bias of the species fraction at the chosen simulation time, not necessarily complete take-over, which can be observed on longer time scales for some parameters, but not for others (see appendix F). From this consistent bias, we thus conclude that, in this competition scenario, the homeostatic active density has a higher impact on fitness than does the homeostatic pressure.

We next examine another prediction of our theoretical consideration, namely that preferential growth of a winner species happens only near interfaces between species. For this purpose, we define 'interface' cells as ones that are less than a distance 15 from a cell of another species (either active or passive) and measure their average growth rates (starting from a time of 1000 (∼14 gen.), before which the system has not coarsened enough to provide good statistics for 'bulk' cells). Figures 6(a) and (b) shows that growth rates at interfaces increase with the difference in ρ^h for species H and decrease for species L. In contrast, 'bulk' growth rates do not show any strong trend. A typical time trace of the different growth rate averages for a case of high homeostatic active density difference is shown in figure 6(c), indicating that this advantage, despite considerable stochasticity and the general coarsening dynamics (cf figure 2), is present throughout the simulation. Collectively, this evidence suggests that preferential growth and consequently the competitive advantage are indeed induced by interfaces, supporting our proposed 'opportunistic growth' mechanism.

Zoom In Zoom Out Reset image size

Finally, we examine the effect of friction on the observed competition mechanisms. In our theoretical considerations (as well as appendix C in more detail), we suggested that, for low friction ζ, the screening length ${l}_{s}=1/\sqrt{-\zeta {\alpha }^{\prime }{\rho }^{\mathrm{h}}}$ diverges, taking the system to the mean-field limit of equation (5), in which different passive matter degradation rates cannot induce an advantage via the opportunistic mechanism described in this work. If there are any other competition mechanisms at play, there must therefore be a crossover to a regime where these become more dominant. As shown in figure 4(b), the homeostatic pressure changes with the degradation factor E due to the inherent change in pressure fluctuations (see appendix E), so that homeostatic pressure is one such other competition driver. To explore this, we repeated the homeostasis and competition simulations once with a friction coefficient ζ₀ 10 times larger than in the previously described simulations, and once with a friction coefficient 10 times smaller. The aggregated results are shown in figure 7.

Zoom In Zoom Out Reset image size

For large friction, the dominant competition driver is still the active density and predicts fitness well (panel (a)) compared to homeostatic pressure (panel (b)). For the sake of the following discussion point, we note that the choice of labeling a species with a higher E value as H together with the monotonous relation between E_H − E_L and ${\Delta}{\rho }^{\mathrm{h}}\equiv {\rho }_{\mathrm{H}}^{\mathrm{h}}-{\rho }_{\mathrm{L}}^{\mathrm{h}}$ (see figure 4(a)) means that Δρ^h is always positive and that, on its own, the opportunistic advantage mechanism can only shift $\langle {N}_{\mathrm{H}}^{\mathrm{f}}/{N}^{\mathrm{f}}\rangle$ upwards as evident in figures 5(c) and 7(a).

For lower friction, we see that the opportunistic competitive advantage mechanism weakens and the correlation between active density difference and fitness deteriorates considerably (panel (c)). Simultaneously, the homeostatic pressure difference becomes a better predictor of fitness (panel (d)), in line with our prediction. However, the correlation of Δ⟨P^h⟩ to fitness at low friction is not as high as the correlations to Δ⟨ρ^h⟩ at higher frictions, and additionally the linear fit in panel (d) does not go through the point (0, 0.5). This suggests that, at the chosen low friction value, the opportunistic advantage mechanism still plays a role, providing some small fitness boost to species of higher active density, which can only lead to an upwards shift of $\langle {N}_{\mathrm{H}}^{\mathrm{f}}/{N}^{\mathrm{f}}\rangle$ values (as noted above) compared to a scenario where the homeostatic pressure mechanism acts exclusively.

In this work, we studied collective competition in growing cellular matter; a common phenomenon in living systems that is responsible both for helpful processes in an organism, such as maintaining homeostasis in an organ, as well as in harmful ones, such as tumor development. We addressed a type of competition between cell types that results from purely mechanically regulated growth, independent of the cell type of interaction partners, i.e. without any adversarial response aimed specifically at competing cell types. Instead, we considered a scenario of competitive exclusion where cells in confinement compete for room to grow. We found that the fitness of a given cell type, which indicates how it will fare during competition, can be determined from its emergent homeostatic properties in a homogeneous system where no competition takes place. In our case of pressure-regulated growth and different passive matter degradation rates, the decisive property turned out to be the density of growing cells (active density).

Evidence for homeostatic cell density as a predictor of fitness can be found in both experiments and simulations [1, 16, 38, 39] and indeed homeostatic density has already been suggested as a competition driver, but only in a very limited sense: one that attributes fitness to a species-dependent threshold density above which apoptosis is triggered [40–42]. The quite intuitive reasoning for the homeostatic density driving the competition in that specific case is that the 'loser' species can be above the threshold and diminishes while the 'winner' species is still below its own apoptosis threshold and prospers. Here, on the other hand, we allow for a more general relationship between the death and division rates that still supports homeostasis. Another distinction is that, in [40–42], division and death rates were considered to depend directly on the density of growing cells, while here we consider these to depend on pressure as seen experimentally in many systems including in epithelia [31], tumours [43, 44], morphogenesis [2], yeast [30] and bacteria [45]. When proliferation and/or removal rates depend on pressure, they are affected by the production of all load-bearing matter in the confined system, and not just by the active growing cells. Thus, the continuous production of passive matter such as dead cells or extracellular matrix as well as its removal process, can affect the cells' ability to compete.

Here, specifically, we considered the case of dead cells, by replacing the instantaneous removal employed in many other works with a finite-time degradation process in both theory and an agent-based numeric model. This has the effect of introducing an explicit, tunable time scale for the persistence of passive dead cells that is separate from the time scale at which cells stop growing and die, and which results in a finite density of passive dead matter in the system.

We showed that a mean-field consideration in the continuum model—describing a uniform mixture of active cells and passive dead ones—predicts the active density in homeostasis well and captures its relation to cell removal (degradation) and death rates (figure 4). The continuum model also provided a typical length scale beyond which cells are mechanically screened from the rest of the system. In a competition scenario, this screening length, combined with the fact that the system is never well-mixed at the cellular scale, leads to a specific growth advantage for the species with the higher active density near interfaces between the two species, as can be qualitatively deduced from a discrete lattice-model consideration (see appendix D). Our results suggest that this advantage emerges through an 'opportunistic growth' mechanism that causes a quicker take-over of vacated space.

By having identical pressure dependencies for growth and identical death rates, the species in our system were originally designed such that the active density could be tuned without affecting homeostatic pressure in order to uncover the specific role of dead matter. However, fluctuations in our agent-based model, controlled by a second 'collapse' time scale of the chosen cell removal process (see appendix E), did in fact lead to second order variations in homeostatic pressure across parameters (figure 4), presenting an unexpected opportunity to test its predictions as well, as it was previously suggested to be a competition driver for systems where homeostasis is regulated by pressure [4, 16]. Homeostatic pressure was shown to explain the suppression of tumour growth under mechanical pressure [4, 29, 44, 46] and can reliably predict fitness in systems where an equation of state exists between the pressure and the active density [16, 47]. Such a relation between active density and pressure implies that the predictions of both homeostatic properties—pressure and density—are equivalent. However, we saw that, when passive matter is continuously produced (and degraded) and growing cells are no longer the only mechanical component, these homeostatic properties are decoupled. The added passive matter of course does not rule out an equation of state, but if one exists, it is between pressure and the total density rather then active density. Since competition outcomes were clearly correlated with increased active density, but not with homeostatic pressure (figure 5), our results do not only provide evidence for homeostatic density as an orthogonal competition driver, but also show that other competition drivers such as homeostatic pressure can be, at least to some extent, overpowered by it. The unambiguous correlation of competition outcomes with active density (despite the non-monotonic behavior of the collapse time scale across the explored parameter range, see appendix E) also suggests that the opportunistic competition mechanism is, at least to some extent, insensitive to the details of the removal process and active density is a suitable effective parameter determining competition.

Other recent studies already showed that homeostatic pressure differences can be overpowered, for example if intra-species and cross-species cell adhesion differ, which can lead to additional growth at interfaces through increased surface tension and thus to spatially structured co-existence [48, 49]. While these features are reminiscent of our observations (with long-term co-existence still being an open question; see appendix F), an interesting difference is that the increased growth in these studies relies on mechanically distinct interfaces that arise from model-inherent, explicit disparity in the mechanical interactions within and across species. In contrast, mechanical interactions in our model are identical between all cells, and it is only the species' intrinsic properties and dynamics that underlie competition.

Having established the competitive advantage of fast dead matter elimination, one could hypothesize that fast removal of dead/non-viable cells could be an important survival mechanism when competing against either rival species or mutations, motivating why, under certain circumstances, processes such as cell extrusion in epithelial tissues might be preferred to slow cell structure degradation. However, it should be noted that the balance between different competition drivers is parameter dependent. As we have shown (figures 7(c) and (d)), when friction is reduced competition outcomes show a stronger correlation with differences in homeostatic pressure than with those in active density—the opposite of the behavior seen in figures 5(c) and (d). This is in line with an increased screening length l_s for lower friction (as also suggested in [32, 50]) and a resulting weaker competitive advantage through the opportunistic growth mechanism, as indicated by our theoretical considerations.

The simulations in this study were performed in 1D for simplicity and to distill the relevant mechanism by avoiding emergent collective phenomena arising, e.g., from anisotropic cell shape, but can easily be generalized to higher dimensionality. Though highly idealized, real-world systems resembling this 1D single-file setup do in fact exist, such as bacteria growing in the 'mother machine' [45, 51, 52], Eukaryotic cells in a capillary [53], and Staphylococcus aureus invasion of narrow bone canals (canaliculi) via fission [54]. In order to address the diversity of real biological systems, future work will have to examine the applicability of the uncovered mechanisms to higher dimensions, where interfaces between domains of different species do not have to be sharp [32, 50, 55, 56] and where higher-dimensional collective effects come into play such as nematic order [57–59], cellular jamming [60] and intermittent load bearing force chains [61–63].

On a final note, we wish to point out that all the forces in our model vary continuously with time (except in the extreme cases E → ±∞ that provide the limits of instantaneous cell removal). To the best of our knowledge, this is the first numerical agent-based model of proliferation and removal in dense multicellular systems (even beyond the context of competition) where forces change continuously without any force jumps, throughout the life-cycle of the cell; i.e. both during cell death and removal and during division (see appendix B). Other numerical models (ranging from agent-based models [25, 33, 47] to active vertex models [26], cellular extended Potts models [64–67] and the subcellular elements method [68]), typically employ instantaneous cell division and/or death for computational simplicity leading to discontinuous forces. We could only find two examples where either division or death were implemented without force jumps [69, 70], but not both. Since the rate at which forces change can affect system properties, as was shown in this paper for the fitness and steady-state pressure, we suggest that instantaneous removal (or division) should be used with care. This is especially true when modeling realistic systems, where these processes are never truly instantaneous, but instead can have a wide range of rates as mentioned in the introduction.

We would like to thank Suropriya Saha, Tunrayo Adeleke-Larodo, Aboutaleb Amiri, Benoît Mahault, and Jérémy Vachier for discussions, feedback, and suggestions. We acknowledge support from the Max Planck School Matter to Life and the MaxSynBio Consortium, which are jointly funded by the Federal Ministry of Education and Research (BMBF) of Germany, and the Max Planck Society.

All data that support the findings of this study are included within the article (and any supplementary files).

Here we provide a more detailed account of the numerical model that was described concisely in the main text. A summary of all the parameters used in the model can be found in table A1. The model is described in 1D as it was used in the simulations, but can be easily generalised to higher dimension.

Table A1. Parameter values used in the main simulations.

Parameter	Value	Description
R	0.5	Node radius
kdeath	0.01	Cell death rate
Δt	5 × 10−4	Time step
⟨Γ⟩	10kdeath	Mean of normal distribution for base growth rates γi
σΓ	⟨Γ⟩/2	Std of base growth rates' normal distribution
${t}_{\mathrm{h}}^{\mathrm{s}\mathrm{i}\mathrm{m}}$	2000	Simulation time for single species runs (∼29 generations)
${t}_{\mathrm{c}\mathrm{o}\mathrm{m}\mathrm{p}\mathrm{e}\mathrm{t}\mathrm{e}}^{\mathrm{s}\mathrm{i}\mathrm{m}}$	5000 a	Simulation time for competition runs (∼72 generations)
H	100	Strength of Hertzian interaction
λ	3.5	Threshold at which growth is considered stopped $\frac{1}{2}\mathrm{e}\mathrm{r}\mathrm{f}\mathrm{c}(\lambda )=3.7\times 1{0}^{-7}$
Pstall	3.16	Pressure at which G(Pi ) drops to $\frac{1}{2}\mathrm{e}\mathrm{r}\mathrm{f}\mathrm{c}(\lambda )$ (corresponding to a node overlap of $\delta {r}_{{i}^{\alpha }{j}^{\beta }}=0.1$)
ζ0, ζ∥, ζaxis	1	Drag coefficients
Tremoval	10	Full duration of cell removal process (independent of E)
Lx	2560 a	Channel length

^aSome additional simulations presented in the appendices (see figures E2, F1 and F2) are run using different simulation times and system sizes as denoted explicitly for each.

Notations

Each cell is modeled as a dumbbell composed of two nodes of radius R that have a (time-dependent) preferred separation. The nodes' positions are given by ${r}_{{i}^{\alpha }}$ where Latin indices identify cells and the α, β indices differentiate between the two nodes of a cell. Without loss of generality we designate a cell's left node as negative(−) and its right node as a positive(+). The separation between any two nodes is given by ${r}_{{i}^{\alpha }{j}^{\beta }}\equiv {r}_{{j}^{\beta }}-{r}_{{i}^{\alpha }}$. For forces between any two nodes, we use the convention that the 1st and 2nd indices specify the node on which the force is applied, whereas 3rd and 4th indices specify the node causing the force:

$$\begin{equation}{f}_{{i}^{\alpha }{j}^{\beta }}\equiv {f}_{{i}^{\alpha }{\leftarrow}{j}^{\beta }},\end{equation} \tag{ A.1 }$$

where the indices obey i ≠ j ∪ α ≠ β, since interactions can only happen between nodes (no self-propulsion of a node).

Inter-cell interactions

The interaction between nodes belonging to different cells is, subject to a node screening rule explained below, a Hertzian repulsion that scales with the node overlap $\delta {r}_{{i}^{\alpha }{j}^{\beta }}\equiv 2R-\vert {r}_{{i}^{\alpha }{j}^{\beta }}\vert$:

$$\begin{equation}{f}_{{i}^{\alpha }{j}^{\beta },i\,\ne \,j}^{\mathrm{i}\mathrm{n}\mathrm{t}\mathrm{e}\mathrm{r}-\mathrm{c}\mathrm{e}\mathrm{l}\mathrm{l}}=\begin{cases}-\mathrm{sgn}({r}_{{i}^{\alpha }{j}^{\beta }})H{\left\vert \delta {r}_{{i}^{\alpha }{j}^{\beta }}\right\vert }^{3/2},\enspace \quad \hfill & \delta {r}_{{i}^{\alpha }{j}^{\beta }}\geqslant 0\hfill \\ 0,\enspace \quad \hfill & \mathrm{O}\mathrm{t}\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{w}\mathrm{i}\mathrm{s}\mathrm{e}.\hfill \end{cases}\end{equation} \tag{ A.2 }$$

Furthermore, a node screening rule is applied to enforce that two nodes of different cells cannot interact across other nodes belonging to the same two cells, which means that, without loss of generality, for a cell q located to the left of a second cell n, the only non-zero interactions can be ${f}_{{q}^{+}{n}^{-}}$ and ${f}_{{n}^{-}{q}^{+}}$. This choice makes physical sense, since overlaps of two cells are but a Hertzian approximation for deformations, and should not allow for interaction across physical objects (i.e. nodes). Yet another reason for this choice is force continuity as discussed in appendix B. When a cell dies and starts degrading, these forces are modulated as explained below in equations (A.5) and (A.7).

Intra-cell forces, growth and division

As opposed to the inter-cell forces, the force between two nodes of the same cell can in principle be either repulsive or attractive (in order to keep structural integrity) and is given by a Hertzian 'spring' with a (growing) rest-length ${\ell }_{i}^{\mathrm{r}\mathrm{e}\mathrm{s}\mathrm{t}}$:

$$\begin{equation}{f}_{{i}^{-}{i}^{+}}^{\mathrm{i}\mathrm{n}\mathrm{t}\mathrm{e}\mathrm{r}\mathrm{n}\mathrm{a}\mathrm{l}}=\pm H{\left\vert {\ell }_{i}^{\mathrm{r}\mathrm{e}\mathrm{s}\mathrm{t}}(t)-{r}_{{i}^{-}{i}^{+}}\right\vert }^{3/2},\end{equation} \tag{ A.3 }$$

where H is a constant that provides the scale of the interaction strength, and the sign of the interaction is determined by $-\mathrm{sgn}({\ell }_{i}^{\mathrm{r}\mathrm{e}\mathrm{s}\mathrm{t}}-{r}_{{i}^{-}{i}^{+}})$. The interaction is reciprocal such that we define ${f}_{{i}^{+}{i}^{-}}^{\mathrm{i}\mathrm{n}\mathrm{t}\mathrm{e}\mathrm{r}\mathrm{n}\mathrm{a}\mathrm{l}}=-{f}_{{i}^{-}{i}^{+}}^{\mathrm{i}\mathrm{n}\mathrm{t}\mathrm{e}\mathrm{r}\mathrm{n}\mathrm{a}\mathrm{l}}$. The similarity of this force-law to the inter-cell repulsion of equation (A.2) is important for temporal force continuity as discussed in appendix B. The rest length ${\ell }_{i}^{\mathrm{r}\mathrm{e}\mathrm{s}\mathrm{t}}$ grows linearly with time from 0 (full overlap of the two nodes) immediately after division to a maximum of 2R (zero overlap). We express the growth of the cell equivalently via a growth phase ${g}_{i}^{\mathrm{g}\mathrm{r}\mathrm{o}\mathrm{w}\mathrm{t}\mathrm{h}}={\ell }_{i}^{\mathrm{r}\mathrm{e}\mathrm{s}\mathrm{t}}/2R$ in the range [0, 1]. The growth rate of the cell is comprised of an inherent base growth rate γ_i that is then modulated by pressure

$$\begin{equation}{\dot{g}}_{i}^{\mathrm{g}\mathrm{r}\mathrm{o}\mathrm{w}\mathrm{t}\mathrm{h}}={\gamma }_{i}\cdot G({P}_{i}),\end{equation} \tag{ A.4 }$$

where in the 1D setup the pressure P_i on a cell is simply the average of the external forces on it in the inward direction. The function G(P_i ) regulates the population growth in order to allow for stable homeostasis. The pressure dependence itself is identical for all cells, while the momentary value depends on the current pressure applied on the cell and obeys G(P_i = 0) ≃ 1 and G(P_i = P^stall) ≃ 0. Thus, for an isolated cell with no outside forces acting on it, the life-cycle time of a cell from division to division is 1/γ_i . In this work we chose G(P_i ) to be a monotonously decreasing sigmoidal function as shown in figure A1 (the equation for G(P_i ) that appears in the main text is repeated in the figure caption for convenience). Since G(P_i ) only goes to zero at the limit of infinite pressure, we define an effective stall pressure. We choose for growth to be effectively stopped if a cell is squeezed by its neighbors with an overlap of 10% the node diameter, corresponding to P^stall ≃ 3.16, and the effective threshold for vanishing growth rate is set by λ = 3.5 (see also table A1).

Zoom In Zoom Out Reset image size

Division

If a cell m avoids death for long enough and its growth phase ${g}_{m}=\int {\dot{g}}_{m}^{\mathrm{g}\mathrm{r}\mathrm{o}\mathrm{w}\mathrm{t}\mathrm{h}}\,\mathrm{d}t$ reaches 1 the cell divides, at which point each of the mother cell's two nodes becomes a new cell, itself composed of two overlapping nodes. Each of the new cells d₁, d₂ starts with a growth phase ${g}_{{d}_{1}}^{\mathrm{g}\mathrm{r}\mathrm{o}\mathrm{w}\mathrm{t}\mathrm{h}}={g}_{{d}_{2}}^{\mathrm{g}\mathrm{r}\mathrm{o}\mathrm{w}\mathrm{t}\mathrm{h}}=0$. The base growth rates of the cells ${\gamma }_{{d}_{1}},{\gamma }_{{d}_{2}}$ are drawn from a distribution Γ that we choose here to be normal (compatible e.g. with bacterial phenotypic growth variability [34, 35, 71]). The mean of Γ is set as 10k^death and the standard deviation as half the mean (negative values and non-growing cells are avoided by truncating Γ at 0.1 of the mean). Additionally, when a partial mother–daughter inheritance scheme for the base growth rates was attempted (not shown), in line with e.g. bacterial proliferative self-regulation [35, 72, 73], only small quantitative differences are evident which do not change any of the conclusions presented in this paper.

Degradation process of dying cells

Cell death is taken to be completely stochastic with a fixed rate k^death per cell. Once cell death is initiated, the cell irreversibly stops growing and cannot divide, transforming it from an active part of the system to a passive one that nevertheless still interacts mechanically for a finite amount of time. To model a continuous degradation process each cell is assigned a scaling factor S_i that modulates the forces it can exert on its neighbors as defined by equation (A.2), such that the net force on a given node is:

$$\begin{equation}{f}_{{i}^{\alpha }}^{\mathrm{n}\mathrm{o}\mathrm{d}\mathrm{e}}={f}_{{i}^{\alpha }{i}^{\beta },\beta \ne \alpha }^{\mathrm{i}\mathrm{n}\mathrm{t}\mathrm{e}\mathrm{r}\mathrm{n}\mathrm{a}\mathrm{l}}+{S}_{i}\sum\limits _{j\ne i}{S}_{j}{\,f}_{{i}^{\alpha }{j}^{\beta },\beta \ne \alpha }^{\mathrm{i}\mathrm{n}\mathrm{t}\mathrm{e}\mathrm{r}-\mathrm{c}\mathrm{e}\mathrm{l}\mathrm{l}}.\end{equation} \tag{ A.5 }$$

The force scaling factor S_i is simply 1 as long as the cell i is alive. Immediately after death, the degradation phase begins. It is characterized by a variable ${g}_{i}^{\mathrm{deg}}$ (not to be confused with the growth phase above), which increases linearly from 0 to 1 over a time period T^removal, i.e. ${g}_{i}^{\mathrm{deg}}=(t-{t}_{i}^{\mathrm{d}\mathrm{e}\mathrm{a}\mathrm{t}\mathrm{h}})/{T}^{\,\mathrm{r}\mathrm{e}\mathrm{m}\mathrm{o}\mathrm{v}\mathrm{a}\mathrm{l}}$. The scaling factors S_i decay as a function of ${g}_{i}^{\mathrm{deg}}$ and finally vanish at ${g}_{i}^{\mathrm{deg}}=1$ (after a time T^removal). The functional form of the scaling factor S_i here is exponential decay:

$$\begin{equation}{S}_{i}=A+B\cdot \mathrm{exp}(-E\cdot {g}_{i}^{\mathrm{deg}}),\end{equation} \tag{ A.6 }$$

which together with the boundary conditions ${S}_{i}({g}_{i}^{\mathrm{deg}}=0)=1$, ${S}_{i}({g}_{i}^{\mathrm{deg}}=1)=0$ leads to:

$$\begin{equation}{S}_{i}=\begin{cases}1,\enspace \quad \hfill & \mathrm{A}\mathrm{l}\mathrm{i}\mathrm{v}\mathrm{e}\hfill \\ 1-\frac{1}{1-{\mathrm{e}}^{-E}}(1-{\mathrm{e}}^{-E\cdot {g}_{i}^{\mathrm{deg}}}),\enspace \quad \hfill & \mathrm{D}\mathrm{e}\mathrm{g}\mathrm{r}\mathrm{a}\mathrm{d}\mathrm{i}\mathrm{n}\mathrm{g},\hfill \end{cases}\end{equation} \tag{ A.7 }$$

as shown in figure 3(a). The rate of cell degradation k^deg that is mentioned in the main text and shown in figure 3(b) is taken as the inverse of the time from death initiation (when S_i = 1) until the force scaling factor S_i reaches a value of 0.1 (see also equation (7)).

Equations of motion

The forces on a cell (as opposed to the individual nodes) can be decomposed into a component acting on the center of mass and a component ${F}_{i}^{\mathrm{a}\mathrm{x}\mathrm{i}\mathrm{s}}$ that drives the node separation degree of freedom:

$$\begin{equation}{f}_{i}^{\mathrm{c}.\mathrm{m}.}={f}_{{i}^{+}}^{\mathrm{n}\mathrm{o}\mathrm{d}\mathrm{e}}+{f}_{{i}^{-}}^{\mathrm{n}\mathrm{o}\mathrm{d}\mathrm{e}}\end{equation} \tag{ A.8a }$$

$$\begin{equation}{f}_{i}^{\mathrm{a}\mathrm{x}\mathrm{i}\mathrm{s}}={f}_{{i}^{+}}^{\mathrm{n}\mathrm{o}\mathrm{d}\mathrm{e}}-{f}_{{i}^{-}}^{\mathrm{n}\mathrm{o}\mathrm{d}\mathrm{e}}.\end{equation} \tag{ A.8b }$$

This leads to the following over-damped (zero Reynolds number limit) Newton equations for cell i in terms of its center of mass velocity ${\dot{r}}_{i}^{\mathrm{c}.\mathrm{m}.}$ and the time derivative of the axial node-separation ${\dot{r}}_{{i}^{-}{i}^{+}}$:

$$\begin{equation}{\zeta }^{0}{\dot{r}}_{i}^{\mathrm{c}.\mathrm{m}.}={f}_{i}^{\mathrm{c}.\mathrm{m}.}/{\zeta }^{{\Vert}}\end{equation} \tag{ A.9a }$$

$$\begin{equation}{\zeta }^{0}{\dot{r}}_{{i}^{-}{i}^{+}}={f}_{i}^{\mathrm{a}\mathrm{x}\mathrm{i}\mathrm{s}}/{\zeta }^{\mathrm{a}\mathrm{x}\mathrm{i}\mathrm{s}},\end{equation} \tag{ A.9b }$$

where ζ⁰ is the overall drag factor, while ζ^∥ and ζ^axis define the relative contributions of the translational drag (parallel to cell long axis), and the axial stretch drag, respectively. The nodes of the cells are spatially propagated according to these equations using a forward Euler integration with a fixed time step 4 orders of magnitude smaller than that of cell growth or cell death (see table A1). Note that these equations do not contain self-propulsion nor any explicit noise terms. Thus single-cell diffusion is completely excluded and the only activity and stochasticity are in the growth, division and death which in turn drive the dynamics.

Initial conditions

The simulations start with living seed cells at a number density of 0.5. Cell center of mass positions are first distributed evenly and then perturbed randomly using a uniform distribution in the range (−0.05, 0.05) in node diameter units. This ensures there are no initial overlaps between the seed cells. The initial growth phase ${g}_{i}^{\mathrm{g}\mathrm{r}\mathrm{o}\mathrm{w}\mathrm{t}\mathrm{h}}$ of each cells is uniformly drawn from the range (0, 1). The nodes of each cell are positioned according to the rest length ${\ell }_{i}^{\mathrm{r}\mathrm{e}\mathrm{s}\mathrm{t}}$ so that initially there are no forces between the 2 nodes of the cell. The base growth rates of the seed cells are drawn from the same normal distribution Γ that is used for dividing cells. In the competition assay, cells are then assigned species identities randomly resulting in a well mixed system.

As mentioned in the discussion, the numerical model was constructed carefully in a way that avoids force discontinuities during cell proliferation, death and removal. Aside from the gradual cell degradation that is discussed extensively throughout the paper, two more strategies are employed to avoid force jumps during cell division. First, the 'node screening rule' for the inter-cell forces that is explained below equation (A.2), i.e., the fact that only the two closest nodes (one from each cell) contribute to inter-cell forces, prevents force jumps that would occur due to the instantaneous doubling of nodes upon cell division. Imagine a cell m is about to divide, while its node α has an overlap $\delta {r}_{{m}^{\alpha }{j}^{\beta }}$ with the node of another cell j. When cell m divides, node α suddenly turns to two overlapping nodes of a new cell d₁, both of which have the same overlap as before with the node of cell j and should thus apply double the force previously applied on that node of cell j. This is prevented by the node screening rule.

The second measure for avoiding force jumps upon division, is replacing the growing rigid shape mechanism that is often used for modeling cell growth [33, 74], by a (here non-linear) spring with a growing rest length and then matching the inter-cell force-law of equation (A.2) to the intra-cell force law of equation (A.3) that one gets just prior to division when ${\ell }_{i}^{\mathrm{r}\mathrm{e}\mathrm{s}\mathrm{t}}=2R$. Thus, the functional form of the node–node interactions does not suddenly change when the nodes switch identity (from belonging to the same mother cell, to belonging to two distinct daughter cells) and neither the forces nor their derivatives jump.

These choices together with the explicit degradation process mean that all the forces on a cell change continuously, without any jumps, throughout the cell life-cycle. The exception is of course in the limits E → ±∞. In practice, one also has to make sure that the time step Δt is small enough, since S_i can change quite sharply for very positive and very negative E values (see figure 3(a)). Furthermore, the parameters of the system have to be chosen to prevent excessive crowding, that is, piling of nodes on top of each other, which is not physically meaningful and breaks the node screening rule (the most important requirement being a pressure dependency of growth which stops growth at low enough densities for a given Hertzian repulsion strength). Another benefit of the gradual removal mechanism is the added realism of the system staying closer to confluency, which can be obtained when the collapse time scale is comparable to the system's growth response and no momentary holes are created when cells degrade (e.g. for E = 4, see supplementary video 1).

In the continuum picture, growing active matter follows a continuity equation like equation (1), where we assume that the division and death rates k^div and k^death have suitable dependencies on the local mechanical environment (e.g., pressure or density), such that homeostasis can be reached. Without loss of generality, we assume that the net growth rate α(P) = k^div(P) − k^death(P) depends on P. For simplicity, we also consider a well-mixed system characterized by a single density ρ, such that

$$\begin{equation}{\partial }_{t}\rho +{\partial }_{x}(v\rho )=\alpha (P)\rho .\end{equation} \tag{ C.1 }$$

This equation has a spatially homogeneous steady-state ρ ≡ ρ^h defined by the condition α(P^h) = 0 corresponding to homeostasis, where P^h = P(ρ^h) via a constitute relation for the pressure. In real systems, usually P ≡ 0 below some critical density ρ_c which corresponds to the closely packed (but uncompressed) density of cells and its functional form above ρ_c is determined by the mechanical properties of the cellular material. While the specific form of P(ρ) is not important here, for consistency we have to require that ρ^h > ρ_c, such that in its vicinity

$$\begin{equation}P({\rho }^{\mathrm{h}}+\delta \rho )\approx P({\rho }^{\mathrm{h}})+{\left.\frac{\partial P}{\partial \rho }\right\vert }_{{\rho }^{\mathrm{h}}}\delta \rho ={P}^{\mathrm{h}}+\frac{1}{\kappa {\rho }^{\mathrm{h}}}\delta \rho ,\end{equation} \tag{ C.2 }$$

where κ is the compressibility of the cellular material at homeostatic density.

To achieve stability, we also require that

$$\begin{equation}{\alpha }^{\prime }={\left.\frac{\partial \alpha }{\partial P}\right\vert }_{{P}^{\mathrm{h}}}< 0,\end{equation} \tag{ C.3 }$$

as this ensures that, again near ρ^h, the net growth rate

$$\begin{equation}\alpha (P({\rho }^{\mathrm{h}}+\delta \rho ))=\alpha ({P}^{\mathrm{h}})+{\left.\frac{\mathrm{d}\alpha }{\mathrm{d}\rho }\right\vert }_{{\rho }^{\mathrm{h}}}\delta \rho ={\left.\frac{\partial \alpha }{\partial P}\right\vert }_{{P}^{\mathrm{h}}}{\left.\frac{\partial P}{\partial \rho }\right\vert }_{{\rho }^{\mathrm{h}}}\delta \rho =\frac{{\alpha }^{\prime }}{\kappa {\rho }^{\mathrm{h}}}\delta \rho \end{equation} \tag{ C.4 }$$

opposes δρ. Local force balance in this over-damped system reads

$$\begin{equation}\zeta v\rho =-{\partial }_{x}P,\end{equation} \tag{ C.5 }$$

where ζ is the friction coefficient with the environment or substrate. Plugging relations (C.2), (C.4) and (C.5) into equation (C.1), we arrive at the following equation for a perturbation ρ = ρ^h + δρ of the homogeneous, homeostatic state:

$$\begin{equation}{\partial }_{t}\delta \rho -\frac{1}{\zeta \kappa {\rho }^{\mathrm{h}}}{\partial }_{x}^{2}\delta \rho =\frac{{\alpha }^{\prime }}{\kappa {\rho }^{\mathrm{h}}}\delta \rho ({\rho }^{\mathrm{h}}+\delta \rho ).\end{equation} \tag{ C.6 }$$

Neglecting higher order terms, we have

$$\begin{equation}\kappa {\partial }_{t}\delta \rho -\frac{1}{\zeta {\rho }^{\mathrm{h}}}{\partial }_{x}^{2}\delta \rho ={\alpha }^{\prime }\delta \rho .\end{equation} \tag{ C.7 }$$

Perturbations therefore both spread diffusively by mechanical relaxation with a diffusion constant D = 1/(κζρ^h) and decay actively (through growth/degradation of cellular matter) with a rate γ = α'/κ. This defines a typical length scale ${l}_{s}=\sqrt{-D/\gamma }=1/\sqrt{-\zeta {\alpha }^{\prime }{\rho }^{\mathrm{h}}}$ which perturbations are able to affect. Note that the independence of l_s from κ is a consequence of α being defined as a function of P (and not, for example, ρ directly), such that both mechanical relaxation and growth adaptation scale equally with κ. For the value of the drag coefficient in our simulations (see table A1), a total homeostatic density of order 0.8 (see, e.g., the right end of figure 4(a)) and α' ≈ ln(2)⟨Γ⟩G'(P^h) ≈ −0.044 with ⟨Γ⟩ = 0.1 (table A1) and P^h ≈ 1.95 (figure 4(b)), this length scale evaluates to l_s ≈ 5.3 in our simulations.

l_s can be interpreted as a screening length: consider a continuous external perturbation (for example, a density sink) which fixes δρ(x = 0) = δρ₀ at the boundary of ${\mathbb{R}}_{\geqslant 0}$. Then, the steady state of equation (C.7) is given by δρ(x) = δρ₀ exp(−x/l_s ), indicating exponential screening of the external perturbation. An increased friction ζ with the substrate (lower mobility) or a stronger change in the growth-degradation balance for a local deviation from the homeostatic pressure (α') will therefore lead to a stronger attenuation away from the perturbation and consequently a shorter l_s . Conversely, in the limit of large mobilities (small friction coefficient ζ) or slow growth responses α', the screening length diverges and the system approaches the mean-field description in which pressure/density equilibration is the fastest process.

We would like to examine the possible influence of the finite granularity of the system on competition dynamics. By granularity we mean the fact that the system is comprised of discrete units, which can be either passive (dead) or active (have the potential to divide). In the main text, we showed that an entirely mean-field, continuum consideration with well-mixed active and passive matter predicts no particular advantage for a cell species which produces more or less passive matter. However, we also showed that, in reality, there exists a finite screening length in the system, so that growth in response to a local perturbation in pressure happens only in a finite region around a perturbation (see appendix C). If this length scale is on the order of the cell size, consequently, only a finite number of discrete cells will respond to a perturbation (for example, the removal of a cell), possibly triggering small-number effects. To investigate the possible consequences of this granularity, we consider the limit of very small screening lengths, in which only neighbors are able to contribute significantly to growth in response to a pressure drop.

To this end, we construct a stochastic lattice model which is characterized by only two different rates: k^death and k^deg, the rate at which an active cell turns into a passive (non-dividing) cell, and the rate at which passive cells are removed from the system, respectively. The role of the growth-limiting mechanical environment (e.g., through pressure) in this model is played by the finite availability of lattice sites, where growth can only occur if a lattice site is freed up by cell removal. This also conserves the property of both the particle-based and the continuum model, where the homeostatic density/pressure do not depend on composition, as the growth behavior of all active cells is identical, irrespective of whether they belong to a species that leaves more or less dead matter behind.

In a single species, lattice sites will be occupied by a mixture of active cells that become passive with rate k^death and passive cells that are removed with rate k^deg. Once a cell is removed, one of the nearby active cells immediately divides and fills the gap (we will clarify below what exactly we mean by 'nearby'). In steady-state, when the total number of cells is constant, the fraction of active cells is therefore

$$\begin{equation}a=\frac{{k}^{\mathrm{deg}}}{{k}^{\text{death}}+{k}^{\text{deg}}},\end{equation} \tag{ D.1 }$$

which is the discrete equivalent of equation (4). At the same time, this is the fraction of time the average cell spends in its active state, or the probability of finding a randomly picked cell to be in an active state.

We consider a one-dimensional configuration where each half-space is filled by a different species. Similar to the particle-based model that we consider later in the main text, they only differ in the rate k^deg with which passive cells are removed. Without loss of generality, we assume that the species on the left has a high degradation rate ${k}_{\mathrm{H}}^{\mathrm{deg}}$ compared to the species on the right with ${k}_{\mathrm{L}}^{\mathrm{deg}}$, which implies a higher equilibrium active fraction a_H > a_L as well.

What then happens at the interface, where a cell of the left species directly neighbors a cell from the right species? First, we note that in the event of the removal of one of the two central cells, a symmetrical distribution of species arises: the gap is lined by H cells to its left and L cells to its right. Since we consider two species with identical growth dynamics (for the active cells), we also assume that it is only the number and position of the active cells of each species that determine which species will fill the gap. Let ${\left\{{r}_{i}^{\mathrm{H}}\right\}}_{i=1,\dots ,\infty }$ be the distances of the active cells of species H, numbered in increasing distance from the interface, and ${r}_{i}^{\mathrm{L}}$ the corresponding distances for species L. Consistent with our screening length consideration, we assume a weighting function w(r) that determines the growth contribution of an active cell at distance r from the interface, with w(r) → 0 as r → ∞. Given a specific configuration characterized by $\left\{{r}_{i}^{\mathrm{H}}\right\}$ and $\left\{{r}_{i}^{\mathrm{L}}\right\}$, the probability that species H will fill the gap is then equal to its total relative growth contribution, i.e. ${\sum }_{i}w({r}_{i}^{\mathrm{H}})/\left({\sum }_{i}w({r}_{i}^{\mathrm{H}})+{\sum }_{i}w({r}_{i}^{\mathrm{L}})\right)$. With this probability for a single given configuration, the probability of species H filling the gap in a large ensemble of such events with arbitrary configurations is

$$\begin{equation}{p}_{\mathrm{H}}^{\text{fill}}=\sum\limits _{\left\{{r}_{i}^{\mathrm{H}}\right\},\left\{{r}_{i}^{\mathrm{L}}\right\}}\frac{{\sum }_{i}w({r}_{i}^{\mathrm{H}})}{{\sum }_{i}w({r}_{i}^{\mathrm{H}})+{\sum }_{i}\,w({r}_{i}^{\mathrm{L}})}\mathrm{P}\mathrm{r}(\left\{{r}_{i}^{\mathrm{H}}\right\},\left\{{r}_{i}^{\mathrm{L}}\right\}),\end{equation} \tag{ D.2 }$$

where the sum goes over all possible configurations that have a non-zero total contribution and $\mathrm{Pr}(\left\{{r}_{i}^{\mathrm{H}}\right\},\left\{{r}_{i}^{\mathrm{L}}\right\})$ is the probability of the configuration $\left\{{r}_{i}^{\mathrm{H}}\right\},\left\{{r}_{i}^{\mathrm{L}}\right\}$ among these. For simplicity, we now assume a sharp cutoff where a maximum of d cells on either side of the interface can contribute to filling the gap (essentially a discrete equivalent of the screening length), i.e.

$$\begin{equation}w(r)=\begin{cases}1\quad r\leqslant d\hfill \\ 0\quad r > d\hfill \end{cases} .\end{equation} \tag{ D.3 }$$

In this case, the sums over the weights in equation (D.2) simply count the active cells in each half-space within distance d from the interface, and the probability of having a specific number of active cells, e.g. k on the left and l on the right, is binomially distributed, such that

$$\begin{equation}{p}_{\mathrm{H}}^{\text{fill}}(d)=\frac{1}{1-{(1-{a}_{\mathrm{H}})}^{d}{(1-{a}_{\mathrm{L}})}^{d}}\sum\limits _{k+l\geqslant 1}\,\frac{k}{k+l}\cdot \left(\genfrac{}{}{0pt}{}{d}{k}\right){a}_{\mathrm{H}}^{k}{(1-{a}_{\mathrm{H}})}^{d-k}\left(\genfrac{}{}{0pt}{}{d}{l}\right){a}_{\mathrm{L}}^{l}{(1-{a}_{\mathrm{L}})}^{d-l}.\end{equation} \tag{ D.4 }$$

To estimate the maximum effect, we explicitly calculate ${p}_{\mathrm{H}}^{\text{fill}}$ for d = 1 from equation (D.4), such that only the nearest neighbors determine growth:

$$\begin{align}\hfill {p}_{\mathrm{H}}^{\text{fill}}(1)& =\frac{1}{1-(1-{a}_{\mathrm{H}})(1-{a}_{\mathrm{L}})}\left(\frac{1}{2}\cdot {a}_{\mathrm{H}}{a}_{\mathrm{L}}+1\cdot {a}_{\mathrm{H}}(1-{a}_{\mathrm{L}})+0\cdot (1-{a}_{\mathrm{H}}){a}_{\mathrm{L}}\right)\hfill \\ \hfill & =\frac{{a}_{\mathrm{H}}(2-{a}_{\mathrm{L}})}{2({a}_{\mathrm{L}}+{a}_{\mathrm{H}}-{a}_{\mathrm{H}}{a}_{\mathrm{L}})}.\hfill \end{align} \tag{ D.5 }$$

For large d, the relative numbers of active cells in each half space become increasingly more narrowly distributed around their expectation values a_H d and a_L d, so we can approximate equation (D.4) as

$$\begin{equation}{p}_{\mathrm{H}}^{\text{fill}}(d)\approx \frac{{a}_{\mathrm{H}}}{{a}_{\mathrm{H}}+{a}_{\mathrm{L}}}+\frac{{a}_{\mathrm{H}}{a}_{\mathrm{L}}}{{({a}_{\mathrm{H}}+{a}_{\mathrm{L}})}^{3}}\frac{{a}_{\mathrm{H}}-{a}_{\mathrm{L}}}{d}.\end{equation} \tag{ D.6 }$$

The second term corresponds to a finite-size correction with respect to the limit ${\mathrm{lim}}_{d\,\to \,\infty }\,{p}_{\mathrm{H}}^{\text{fill}}(d)=\frac{{a}_{\mathrm{H}}}{{a}_{\mathrm{H}}+{a}_{\mathrm{L}}}$, where exactly a fraction of a_H cells in the left half-space and a fraction of a_L cells in the right half space are active. This shows that, in any symmetrical finite sample of cells around an interface, the growth contribution from active cells is statistically biased in favor of the species with the higher active fraction a_H. Note that the limit for infinite d should hold for any reasonable weighting function w(r) where d represents the scale on which cells contribute significantly to growth (not just equation (D.3)).

Having determined the probability with which each side will fill the gap, we need to calculate the frequency with which gaps are created in the first place. Overall, species H will only have an advantage if, on average, gaps created by the removal of L cells are filled by H cells with a higher rate than gaps created by the removal of H cells are filled by L cells. Given a cell that is passive (which is the case with probability 1 − a), it is removed with rate k^deg, leading to a total gap creation rate for a single site of

$$\begin{equation}{k}^{\text{gap}}=(1-a){k}^{\text{deg}}=a{k}^{\text{death}}=\frac{{k}^{\text{death}}{k}^{\text{deg}}}{{k}^{\text{death}}+{k}^{\text{deg}}}.\end{equation} \tag{ D.7 }$$

Consistently, for very large k^deg (i.e. no passive cells), we simply have k^gap = k^death, whereas for small k^deg, removal of passive cells is the bottleneck. Note that ${k}_{\mathrm{H}}^{\text{gap}} > {k}_{\mathrm{L}}^{\text{gap}}$, i.e. gaps are created more frequently on the side with a higher percentage of active cells. However, we have to consider both gap creation and cell replacement to calculate the net bias. Consequently, the total rate k_L→H with which cells of type L are replaced with cells of type H directly at the interface is

$$\begin{equation}{k}_{\mathrm{L}\to \mathrm{H}}(d)={k}_{\mathrm{L}}^{\text{gap}}{p}_{\mathrm{H}}^{\text{fill}}(d)-{k}_{\mathrm{H}}^{\text{gap}}{p}_{\mathrm{L}}^{\text{fill}}(d).\end{equation} \tag{ D.8 }$$

For the d = 1 case, this evaluates to

$$\begin{equation}{k}_{\mathrm{L}\to \mathrm{H}}(1)={k}^{\text{death}}\frac{{a}_{\mathrm{H}}{a}_{\mathrm{L}}({a}_{\mathrm{H}}-{a}_{\mathrm{L}})}{2({a}_{\mathrm{L}}+{a}_{\mathrm{H}}-{a}_{\mathrm{H}}{a}_{\mathrm{L}})}\approx \frac{1}{2}{k}^{\text{death}}({a}_{\mathrm{H}}-{a}_{\mathrm{L}}),\end{equation} \tag{ D.9 }$$

where the last step is valid for a_H, a_L close to 1. Therefore, for d = 1 (and, in fact, for any finite d), the lower ${k}_{\mathrm{L}}^{\text{gap}}$ with which species L opens a gap is more than outweighed by species H's statistically increased probability to fill it. This creates a bias towards species H which has a higher active fraction a_H. Note that by using equation (D.6) we can see that for large d,

$$\begin{equation}{k}_{\mathrm{L}\,\to \,\mathrm{H}}(d)\approx \frac{1}{4d}{k}^{\text{death}}({a}_{\mathrm{H}}-{a}_{\mathrm{L}}),\end{equation} \tag{ D.10 }$$

and consequently,

$$\begin{equation}\underset{d\,\to \,\infty }{\mathrm{lim}}\,{k}_{\mathrm{L}\to \mathrm{H}}(d)=0,\end{equation} \tag{ D.11 }$$

which underscores the nature of this bias as a small-number effect due to the granularity of the system in conjunction with the screening length.

It should be noted that the correspondence between the true, agent-based model (or a continuum model) and this lattice model should be regarded as conceptual rather than formal. While d above does play the role of a screening length in the sense that it allows only cells within a certain distance of the interface to replenish the system after the removal of an interface cell, there are many differences which forbid the interpretation of these lattice model results as a quantitative prediction. Most important among these are the instantaneous (rather than continuous) growth and the absence of explicit mechanical relaxation (which is only incorporated in an effective sense through the screening length). For example, with a more continuous growth process in mind, the probability p_H of filling the gap with a cell of type H should rather be viewed as the fraction of the gap which is filled with cellular material of type H. Therefore, this analysis should be interpreted as 'supporting evidence' for our hypothesis that a finite screening length in conjunction with the discreteness of cells can lead to a biased selection of active cells from the more active species, rather than a formal proof.

As already pointed out in the main text, the homeostatic pressure measured in the simulations (figure 4(b)) shows a deviation from the constant value predicted by the zeroth-order mean-field approximation ${\tilde{P}}^{\mathrm{h}}={\mathcal{J}}^{-1}({k}^{\mathrm{d}\mathrm{e}\mathrm{a}\mathrm{t}\mathrm{h}})$, where $\mathcal{J}$ defines the pressure dependency of the division rate, ${k}^{\mathrm{d}\mathrm{i}\mathrm{v}}\equiv \mathcal{J}(P)$ (and a tilde denotes a mean-field quantity). The reason for this is the degradation of cells, which entails a constant turnover of matter. Due to the screening length, degrading cells cause local deviations from the mean-field uniform pressure derived in the continuum description in the main text.

First, let us derive this zeroth-order prediction of the homeostatic pressure (i.e. where division and death exactly balance) in the context of the agent-based numerical model: in homeostasis, the probabilities for death and for survival until division have to be equal. Cells die with a constant rate k^death, leading to an exponentially distributed time of death, such that the condition for the survival probability becomes p^survival = exp(−k^death τ^div) = 1/2, where τ^div is the resulting (pressure modulated) division time in homeostasis. This leads to ${k}^{\mathrm{d}\mathrm{e}\mathrm{a}\mathrm{t}\mathrm{h}}=\mathrm{l}\mathrm{n}(2)/{\tau }^{\mathrm{d}\mathrm{i}\mathrm{v}}\simeq \mathrm{l}\mathrm{n}(2){\bar{\dot{g}}}^{\mathrm{g}\mathrm{r}\mathrm{o}\mathrm{w}\mathrm{t}\mathrm{h}}\simeq \mathrm{l}\mathrm{n}(2)\bar{\gamma }G(P)$, where $\bar{\gamma }$ is the average base growth rate (which we approximate below with the mean of the distribution of growth rates for new cells Γ). Inverting G one gets for the homeostatic pressure

$$\begin{equation}{\tilde{P}}^{\mathrm{h}}\simeq {G}^{-1}\left({k}^{\mathrm{d}\mathrm{e}\mathrm{a}\mathrm{t}\mathrm{h}}/(\mathrm{ln}(2)\langle {\Gamma}\rangle )\right)=1.919.\end{equation} \tag{ E.1 }$$

Second, we examine how the time-averaged balance between the division rate and the death rate is modified by local pressure fluctuations P = P^h + δP about the unknown homeostatic pressure P^h, given the nonlinear dependence of division on pressure through $\mathcal{J}$. For this purpose, we expand the time-averaged division rate

$$\begin{equation}{\bar{k}}^{\mathrm{d}\mathrm{i}\mathrm{v}}=\mathcal{J}({P}^{\mathrm{h}})+\frac{\mathrm{d}\mathcal{J}}{\mathrm{d}P}\bar{\delta P}+\frac{1}{2}\frac{{\mathrm{d}}^{2}\mathcal{J}}{{\mathrm{d}P}^{2}}\bar{{(\delta P)}^{2}}+\mathcal{O}(\bar{{(\delta P)}^{3}}).\end{equation} \tag{ E.2 }$$

The term linear in δP drops out in the time averaging and extracting the homeostatic pressure we get

$$\begin{equation}{P}^{\mathrm{h}}\simeq {\mathcal{J}}^{-1}\left({k}^{\text{death}}-\frac{1}{2}\frac{{\mathrm{d}}^{2}\mathcal{J}}{{\mathrm{d}P}^{2}}\bar{{(\delta P)}^{2}}\right).\end{equation} \tag{ E.3 }$$

Since $\mathcal{J}$ is a decreasing function of pressure and $\frac{{\mathrm{d}}^{2}\mathcal{J}}{{\mathrm{d}P}^{2}}$ is positive (see figure A1 with 0.61 ⩽ P^h/P^stall ⩽ 0.62 in simulation), we see that, for larger local fluctuations, the pressure increases (if one can neglect the change of ${\left.\frac{{\mathrm{d}}^{2}\mathcal{J}}{{\mathrm{d}P}^{2}}\right\vert }_{{P}^{\mathrm{h}}}$ itself).

With this expectation in mind, it is worth re-examining figure 4(b). If local pressure fluctuations are small, we expect a homeostatic pressure equal to the zeroth-order prediction (E.1) indicated by the dashed orange line. The fact that the homeostatic pressure is always above this value (the minimum value measured being 1.924) indicates that the fluctuations are always relevant for the chosen cell removal mechanism and parameters used. Furthermore, we expect that the fluctuations would be smallest at around E = 4 where the measured homeostatic pressure is minimal, and for the fluctuations to increase away from this value towards both higher and lower E values ⁶ .

Indeed, looking at the global pressure fluctuations in figure E1(a), we see that they are minimal at around the same E value. Note that a global pressure fluctuation represents, at any point in time, the sum of contributions of individual local pressure fluctuations with only weak spatial correlations beyond the screening length. Thus, global fluctuations should decay $\propto \!1/\sqrt{{L}_{x}}$ for system sizes L_x much larger than the screening length. This is indeed the case as shown in figure E2(b).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Now that we know that (up to said $1/\sqrt{{L}_{x}}$ scaling) global pressure fluctuations are a good proxy for how pressure fluctuates on the scale of the screening length, we would like to show how changes in E, i.e. the dynamics of the degradation mechanism, affect the latter local pressure fluctuations. For this purpose, we consider the change in length of a single degrading cell with time. As an approximation, we start with the rest length of a cell and modulate it by the decaying force scaling factor S_i while keeping the effects of the surrounding cell medium fixed (such as the pressure P^ref applied by it on the cell):

$$\begin{align}\hfill {\tilde{\ell }}_{i}& =2R+{r}_{{i}^{\alpha }{i}^{\beta }}-\delta {r}_{{i}^{\alpha }{(i+1)}^{\beta }}-\delta {r}_{{i}^{\alpha }{(i-1)}^{\beta }}=2R+{g}_{i}^{\mathrm{i}\mathrm{n}\mathrm{i}\mathrm{t}}-\delta {r}_{{i}^{\alpha }{i}^{\beta }}-\delta {r}_{{i}^{\alpha }{(i+1)}^{\beta }}-\delta {r}_{{i}^{\alpha }{(i-1)}^{\beta }}\hfill \\ \hfill & =2R+{g}_{i}^{\mathrm{i}\mathrm{n}\mathrm{i}\mathrm{t}}-{({P}^{\mathrm{r}\mathrm{e}\mathrm{f}}/H)}^{2/3}(1+{({S}_{i}(E,{g}_{i}^{\mathrm{deg}}))}^{-2/3}),\hfill \end{align} \tag{ E.4 }$$

where the indices i ± 1 denote immediate (living) neighbors and where for the third equality the overlaps are calculated using equations (A.2), and (A.5) which give ${P}^{\mathrm{r}\mathrm{e}\mathrm{f}}\simeq {P}_{i}={S}_{i}{S}_{i\pm 1}H{(\delta {r}_{{i}^{\alpha }{(i\pm 1)}^{\beta }})}^{3/2}$ and ${P}^{\mathrm{r}\mathrm{e}\mathrm{f}}\simeq {P}_{i}=H{(\delta {r}_{{i}^{\alpha }{i}^{\beta }})}^{3/2}$. Here ${g}_{i}^{\mathrm{i}\mathrm{n}\mathrm{i}\mathrm{t}}$ is the growth phase at which the cell died and started degrading. Examples of the decay of this length with time are shown in figure E1(b) for several E values (using representative values for ${g}_{i}^{\mathrm{i}\mathrm{n}\mathrm{i}\mathrm{t}},{P}^{\mathrm{r}\mathrm{e}\mathrm{f}}$). For very negative E values, the cell length decays only slowly for most of the degradation time before finally dropping sharply. We call this sharp drop in length the 'collapse' stage and associate to it a large local pressure drop. Such collapse could describe for example the loss of structural integrity due to tearing of the cell membrane in eukaryotes or breaking of the cell wall in bacterial lysis [24]. For very high positive E values, the cell starts collapsing as soon as it dies. As a measure of the collapse time scale τ^coll we consider the time it takes the length of a degrading cell to decay from 0.9 to 0.1 of its initial value (when S_i = 1):

$$\begin{equation}{\tau }^{\mathrm{c}\mathrm{o}\mathrm{l}\mathrm{l}}(E)={t}_{\,\tilde{\ell }=0.1{\tilde{\ell }}^{\mathrm{i}\mathrm{n}\mathrm{i}\mathrm{t}}}-t{\vert }_{\tilde{\ell }=0.9{\tilde{\ell }}^{\mathrm{i}\mathrm{n}\mathrm{i}\mathrm{t}}}.\end{equation} \tag{ E.5 }$$

Though this measure for the collapse time scale depends on the growth phase of a cell when it died ${g}_{i}^{\mathrm{i}\mathrm{n}\mathrm{i}\mathrm{t}}$, we verified that, for our model, it only does so weakly. A visualization of the calculated collapse time scale is shown in figure E1(c). It is reasonable to associate faster cell collapse (small τ^coll) with larger local pressure fluctuations. Indeed, a simple linear transformation of τ^coll matches the measured global pressure fluctuations as shown in figure E1(a), which we showed above to be a proxy for the local fluctuations. This suggests that pressure fluctuations in our system are indeed determined by the characteristics of the degradation process. It is thus apparent that the chosen removal strategy of mechanical softening in fact introduces two time scales that are relevant for competition; (1) the inverse of the degradation rate k^deg which controls the opportunistic growth mechanism that is exposed in this work, and (2) the collapse time τ^coll, which via the pressure fluctuations controls the homeostatic pressure advantage. The fact that these two time scales change inter-dependently with E and that the former is monotonous in E (figure E1(d)) while the latter changes non-monotonously (figure E1(c)), allows us to explore also the cases where these two competition drivers are opposed as discussed in the main text; something that would not be possible with, for example, a simple linear time-decay of S_i .

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We remind the reader that the original simulation time of 5000 (∼72 gen) for the competition assay was chosen to be short enough to avoid measurement bias of the fitness caused by the absorbing boundary (where an entire sample is taken over by a single species). At longer times, at least for some of the species pairs, the competition does not saturate at some finite species fraction and the less fit species is indeed eradicated from the system at the long-time limit. Specifically we see this when the difference in active density is high enough and when the difference in homeostatic pressure also favors the same species as shown in figures F1(a) and (b) for simulations that are 8 times as long as the ones in the main text (t^sim = 40 000, ∼576 gen). Conversely, some other species pair combinations do not result in extinction of the loser despite a clear initial disadvantage, as shown for example in figures F2(a) and (b). One possible explanation for this is that the effectiveness of the opportunistic competition mechanism, acting inherently only at interfaces as shown above, diminishes as the system coarsens and the number of interfaces decays. It is thus conceivable that, at some point, this mechanism and the one of the homeostatic pressure mechanism become equally strong, and, when acting oppositely, can lead to co-existence. However, it cannot be easily determined whether a stable co-existence state exists since the species ratio at the long-time limit differs significantly for different realizations of the same simulation and, in fact, does not seem to converge to any steady-state value. It is thus also possible that, for the relevant species pairs, the time scale for extinction is simply longer than the longest simulation time attempted here. Indeed, for the values of E_H, E_L where the two competition mechanisms are opposed, we expect a slower sample take-over. One can see this when comparing figures 4(a) and (b); under the constraint of opposed competition mechanisms (P^h for species H being lower than for species L), the largest difference of ρ^h is around half the maximal possible difference, which means that, when the two competition drivers are opposed, we are also exploring a subset of the parameter range where the fitness advantage of the opportunistic growth mechanism alone is in any case reduced.

Finally, is also useful to look at the species fraction dynamics at start of the competition assay as shown in figure F1(c). One can see that the fraction of species H can go below 0.5 due to the stochastic nature of the dynamics, giving way to species L, before increasing again and eventually winning. This shows that the dynamics is robust to the choice of initial conditions and that an initial larger fraction of L cells will not change the identity of the winner.