In our experiments, we examined competition between wild-type Madin–Darby Canine kidney epithelial cells (winners, MDCK^WT) and cells depleted for the polarity protein scribble (losers, MDCK^Scrib) (Norman et al., 2012). To allow for simple image analysis, each cell type stably expressed a histone marker fused to a different fluorophore (MDCK^WT:GFP and MDCK^Scrib:mRFP). Cells were seeded in various ratios of loser:winner cells (10:90, 50:50, 90:10) as well as colonies and imaged for up to 96 hr at 4 min intervals (Figure 1A). Cell segmentation and tracking allowed to determine population measurements (such as the evolution of cell count, the number of mitoses, and the number of apoptoses) as well as cellular-scale measurements (such as local cell density, number of neighbours, identity of neighbours, and cell state) for each cell type (Bove et al., 2017). These data provided the metrics to compare simulations to experiments (Figure 1B, C).

Figure 1 with 5 supplements see all

Download asset Open asset

To understand the emergence of cell competition, we implemented a multi-scale agent-based model that simulates mechanical interactions between cells and with their underlying substrate (Figure 1B, Figure 1—figure supplement 1, grey shaded area) and implements cell-autonomous decisions for growth, mitosis, and apoptosis (Figure 1B, Figure 1—figure supplement 1, pink shaded area). In contrast to existing computational approaches (Rejniak and Anderson, 2011; Zhang et al., 2009), our model includes the coupling between cellular mechanics and decision-making.

In our simulation, epithelial cells are modelled using a cellular Potts model (CPM) (Graner and Glazier, 1992), which enables physical interactions at the cell-cell and cell-substrate interfaces to be simulated (Figure 1—figure supplements 1 and 2, Materials and methods). This implementation was preferred to the less computationally costly vertex model (Fletcher et al., 2014) because we compare our model to our in vitro competition experiments (Bove et al., 2017) that start from a sub-confluent state.

In the CPM, each cell ${\sigma }^k$ possesses a cell type $\tau$ (winner or loser) and is represented by a set of pixels (i,j). The free energy of a group of N cells ${\sigma }^k$ is given by the Hamiltonian $H$:

The energy function represents energetic contributions due to intercellular adhesion, cell adhesion to the extracellular matrix, cell elasticity, and active cell movement. The first term in the Hamiltonian describes the adhesive interactions between cells at their shared interface, where δ represents the Kronecker delta function and Θ is the Heaviside theta function (see Materials and methods). When a cell interacts with other cells, it engages in either homotypic adhesion if they are of the same cell type or heterotypic adhesion if they are not. The respective adhesion energies are given by J_homotypic and J_heterotypic. From a biological perspective, $J$ is the difference between the surface tension and intercellular adhesion. Therefore, higher $J$ implies lower intercellular adhesion. Cells can also interact with the substrate at their periphery via integrin binding to the extracellular matrix, parametrised by an adhesion energy J_{cell-substrate}. The second term describes an elastic energy with an elastic modulus $\lambda$ arising from the cytoskeleton. This energy scales with the difference between a cell’s actual area $A\left({\sigma }^k\right)$ and its target area $A_T\left({\sigma }^k\right)$, which it would occupy in the absence of crowding due to other cells. The actual area is determined by the height and volume of a cell. The third term reflects energy due to cell motion and is parametrised by a kinetic energy ${\lambda }_m$ and a unit polarity vector $\widehat{m}\left(\sigma ,t\right)$ that defines the direction of cell motion and undergoes rotational diffusion. An empirical conversion between computational time and experimental time was obtained by comparing the mean square displacement of isolated cells in experiments and simulations.

In addition to the CPM that determines cell shapes based on mechanical equilibrium, a second computational layer based on cell automata rules regulates changes in cell size due to growth and implements changes in cell fate (division and apoptosis) (Figure 1B, Figure 1—figure supplement 1, pink shaded area). It is in this layer that cellular decision-making is implemented at each time point based on a set of probabilistic rules that we determine from our experimental data (Materials and methods). We now briefly describe these rules and the calibration of their associated parameters (Figure 1—figure supplement 2A).

Experimental work has shown that MDCK cells maintain cell size homeostasis by following an ‘adder’ mechanism, in which each cell cycle adds a set volume to the cell (Cadart et al., 2018; Figure 1—figure supplement 1B). Following the adder model, we increase each cell’s target area A_T(t) at each time point t and cells divide when a threshold area $∆A_{tot}$ has been added since the start of their cell cycle (Figure 1—figure supplement 3). The target area of cells at birth A_T(0), the threshold area added at each cell cycle $∆A_{tot}$, and the maximum growth rate G were all calibrated from movies of isolated cells to reflect the cell cycle time and size distribution measured in experiments (Figure 1—figure supplement 2A). Above a certain cell density, proliferation ceases – a phenomenon known as contact inhibition of proliferation (Abercrombie, 1979; Lieberman and Glaser, 1981). Arrest in proliferation is accompanied by a decrease in protein synthesis due to a drop in ribosome assembly and downregulation of the synthesis of cyclins (Azar et al., 2010). We incorporated contact inhibition of proliferation (Figure 1—figure supplement 1E) by making the target area growth rate $d{A}_T\left(t\right)/dt$ dependent on the difference between the actual cell area A(t) and the target area A_T(t) as $d{A}_T\left(t\right)/dt=G{e}^{-k{\left(A\right(t)-A_T(t\left)\right)}^2}$, where $G$ is the growth rate in the absence of crowding and $k$ is a heuristic parameter that quantifies the sensitivity to contact inhibition. As a result, cellular growth rate slows down exponentially as crowding increases, leading to an increase in cell cycle time.

We implemented two separate rules for cell apoptosis in mechanical and biochemical competition (Figure 1—figure supplement 1C). First, under crowded conditions where mechanical competition is dominant, the probability of apoptosis p_apo increases with local cell density ρ following a sigmoid curve (Figure 1—figure supplement 4A). Our experimental data suggested that winner cell apoptosis followed the same law as losers except shifted towards higher local densities (Bove et al., 2017). Second, in biochemical competition, experimental data indicates that the probability of apoptosis of loser cells depends on the percentage of their perimeter in contact with the winner cells (heterotypic contact, Figure 1—figure supplement 5A; Levayer et al., 2015). This was implemented as a Hill function as a function of percentage of perimeter occupied by heterotypic contact (Materials and methods). In the absence of experimental measurements, we chose a maximum probability ${\displaystyle p_{apomax}}$ of death per frame of a similar magnitude to that measured in mechanical competition. This is justified by the fact that mechanical and biochemical competition takes place over comparable durations in MDCK cells ~2–4 days (Hogan et al., 2009; Kajita et al., 2010).

In experiments, cell elimination can also occur through live cell extrusions when the cell apical area decreases substantially compared to the population average (Eisenhoffer et al., 2012; Kocgozlu et al., 2016). In our simulations, cells delaminated when their actual area A(t) became smaller than ${\displaystyle \frac{<A_i>}{2}}$, with <A_i> the average area of all cells in the simulation at time t (Figure 1—figure supplement 1D).

Taken together, the combination of cellular mechanics and decision-making strategies provides a multi-scale agent-based model to investigate how the interplay between short-range and long-range competitive interactions determines tissue composition. Many parameters are used to describe each cell type’s behaviour, some can be measured directly from experiments while others must be empirically determined based on comparison of the output of the simulations and the experimental data (Figure 1—figure supplement 2A, Supplementary files 1 and 2). We sought to fix as many parameters as possible to restrict the parameter space explored.

To validate the predictive power of our model, we first simulated homeostasis in pure cell populations of winner and loser cells undergoing proliferation and apoptosis. To calibrate our model parameters (Figure 1—figure supplement 2A), we compared simulations to experiments on the basis of the temporal evolution of population metrics, such as cell count and average cell density (Figure 2C, D, F, G). In addition, as experimental and theoretical work has shown that cell organisation in monolayers can be described by the distribution of number of neighbours each cell possesses (Gibson et al., 2006) and their area relative to the population mean (Farhadifar et al., 2007), we also used these cell-scale metrics for confluent epithelia (Figure 2E, H, Figure 2—figure supplement 1).

Figure 2 with 2 supplements see all

Download asset Open asset

Video 1

Download asset

At the start of our experiments, many isolated cells can be observed (t=0h, Figure 2B, Figure 2—figure supplement 2A–C). In pure populations, MDCK^Scrib cells spread markedly more than MDCK^WT (Figure 2—figure supplement 2B, D, compare MDCK^WT in Video 2 to MDCK^Scrib in Video 3), consistent with Wagstaff et al., 2016. Interestingly, in mixed populations dominated by WT cells, MDCK^Scrib cell area diminished compared to pure populations even prior to confluence (Figure 2—figure supplement 2C–E). We used these measurements to set the distribution in cell areas at birth A_T(0) and the area added at each cell cycle $∆A_{tot}$ for each cell type in pure and competitive conditions. Transcriptomic data comparing both cell types indicates that MDCK^Scrib do not express more integrins than MDCK^WT (Wagstaff et al., 2016). Therefore, we assigned the same value of J_{cell-substrate} to both cell types.

Video 2

Download asset

Video 3

Download asset

The maximum growth rate G was parameterised based on the measured distribution of cell cycle durations prior to confluence (Bove et al., 2017), which showed that losers grew significantly slower than winners (mean cell cycle times: MDCK^Scrib ~21.6 hr vs. MDCK^WT ~18 hr). Therefore, we assigned a smaller G to losers than to winners.

After confluence, cell shape is controlled by the interplay between intercellular adhesion energy J_homotypic and the stiffness modulus λ. The value of J_homotypic was adjusted such that, at confluence, the distributions in apical area, number of neighbours, and area relative to the population mean matched experiments (Figure 2—figure supplement 1A, B), as done by others (Farhadifar et al., 2007). An accurate replication of sidedness of cells is particularly important for simulating biochemical competition because the probability of apoptosis of loser cells is linked to the fraction of their perimeter contacting winner cells (Levayer et al., 2015). Previous experimental work showed that, in competition experiments, the area of loser cells became significantly smaller than that of winners when monolayers reached high densities post-confluence (Wagstaff et al., 2016; Bove et al., 2017). Therefore, we assigned losers a smaller value of λ than winners.

The strength of contact inhibition, $k$ (Materials and methods), was adjusted empirically such that the average local cell density at long time scales in simulations reached a plateau that matched the one observed in experiments (Figure 2D, G).

In mechanical competition experiments, loser cells undergo apoptosis when they are in crowded environments. Therefore, we implemented a relationship between probability of apoptosis per cell per unit time (p_apo) and the local cell density ρ parametrised by fitting our experimental data with a sigmoid function (Materials and methods, Figure 1—figure supplement 4A; Bove et al., 2017). The local cell density ρ was defined as the inverse of the sum of the area of the cell of interest and its first neighbours (Materials and methods, Figure 1—figure supplement 4B). While p_apo saturates at high densities for loser cells, experimental data was not available for the highest densities for winner cells (data points, Figure 1—figure supplement 4A). Therefore, we assumed that the maximum probability of apoptosis p_{apo max} was the same for winner and loser cells (Figure 1—figure supplement 4A).

When the simulations were initialised with the calibrated parameters and the same initial cell number as in experiments, cell count and density in the simulations qualitatively reproduced our experimental observations (Figure 2A, B, Videos 1 and 2). MDCK^WT cell count increased for ~70 hr before reaching a plateau at a normalised cell count of 5.5, indicative of homeostasis (Figure 2C). The temporal evolution of the average local cell density and the distribution of number of neighbours at confluence were also faithfully replicated by our simulations (Figure 2D, E). Similarly, our parametrisation of MDCK^Scrib accurately replicated the temporal evolution of cell count and density, as well as the distribution of the number of cell neighbours (Figure 2F–H, Video 3). In particular, the loser cell count and density stayed fairly constant throughout the whole simulation and experiment.

The maintenance of an intact barrier between the internal and the external environment is a key function of epithelia. This necessitates exact balancing of the number of cell deaths and divisions. Failure to do so results in hyperplasia, an early marker of cancer development. Previous work has revealed that epithelia possess a preferred density to which they return following perturbation, signifying that they seek to maintain a homeostatic density (Eisenhoffer et al., 2012; Marinari et al., 2012; Gudipaty et al., 2017). In the experiments performed in Eisenhoffer et al., 2012, cells were grown to confluence on stretchable substrates and subjected to a step deformation in one axis. When deformation increased cellular apical area, the frequency of cell division increased (Gudipaty et al., 2017), while a decrease in apical area resulted in increased live cell extrusion and apoptosis (Eisenhoffer et al., 2012). Therefore, the existence of a homeostatic density is an essential property of epithelia that relates to their sensitivity to crowding – a key factor in mechanical competition. However, current models of epithelia do not implement this.

In our model, we implemented two mechanisms shown experimentally to decrease cell density: cell extrusions and density-dependent apoptoses (Figure 1—figure supplement 1C, D, Materials and methods). We simulated the response of a confluent epithelium to a sudden 30% increase in homeostatic density. In experiments on confluent MDCK^WT epithelia (Eisenhoffer et al., 2012), a sudden increase in crowding was followed by a gradual decrease in cell density resulting from a combination of apoptoses and live cell extrusion, before returning to the initial homeostatic density after ∼6 hr (green data points, Figure 3A).

Figure 3

Download asset Open asset

In our simulations, we allowed MDCK^WT cells to reach their homeostatic density before suddenly increasing the cell density by a percentage similar to experiments. In response to this, cell density decreased gradually over a period of 6 hr with dynamics similar to those determined experimentally (solid black line, Figure 3B). Thereafter, cell density remained at the homeostatic density for long time periods (Figure 3B). Thus, our model implementation and parametrisation replicate the return to homeostasis of pure populations of MDCK^WT cells. As the relationship between p_apo and ρ was fitted using the data from our competition experiments, this also raises the possibility that cell apoptosis in response to crowding is a cell-autonomous process rather than specific to cell competition.

Our experimental work indicated that, in competitions between MDCK^WT and MDCK^Scrib, two processes might be at play, a density-dependent apoptosis of MDCK^Scrib and an upregulation of division of MDCK^WT in MDCK^Scrib-dominated neighbourhoods (Bove et al., 2017). The former is central to mechanical competition, while the latter implies that contact between cell types may control division rate. To determine which process was dominant, we tested whether cell-type differences in density-mediated apoptosis alone were sufficient to explain competition.

We used our model of winner and loser cells with different sensitivities to crowding parametrised from experiments on pure cell populations (Figure 2, Supplementary file 1). We initialised our simulations by seeding a 90:10 winner-to-loser cell ratio, as in experiments. Our simulations were able to quantitatively reproduce the experimental data for competition dynamics, with no further adjustment in parameters. As in the experiments, simulated winner cells (green) rapidly proliferated while loser cell numbers (red) increased weakly until ~50 hr before diminishing (Figure 4A, B, Videos 4 and 5). Furthermore, the evolution of cell count was quantitatively replicated over the entire duration of the experiment for both winner and loser cells (Figure 4C). Cumulative divisions and apoptoses in simulations closely matched those observed in experiments (Bove et al., 2017; Figure 4—figure supplement 1B, C). One of the most striking features of experimental data is that the local density of loser cells in competition increases dramatically compared to pure populations (approximately fivefold increase, Figure 4—figure supplement 1D for comparison), while the local density of winner cells in competition follows the same trend as in pure populations (Figure 4—figure supplement 1E for comparison) (Bove et al., 2017). The sharp increase in local density of loser cells is replicated in our simulations (red curve, Figure 4D) and likely arises from their lower stiffness modulus $\lambda$. In addition, the probability of apoptosis and division as a function of density computed from simulation data matched the experimentally measured ones for both cell types (Figure 4E, F). While the former is an input to our simulation, the latter is an output. Finally, when we compared the probability of division of winner cells in contact with at least one loser cell to that of winner cells in contact with only winner cells, we found that the probability of division of winner cells increased when they were in contact with loser cells (Figure 4G). This was consistent with our experimental observations and occurred despite the fact that we did not implement any aspect of biochemical signalling in our simulations, signifying it represents an emergent property.

Figure 4 with 1 supplement see all

Download asset Open asset

Video 4

Download asset

Video 5

Download asset

Overall, differences in density-dependent apoptosis alone are sufficient to replicate the evolution of cell count and density observed in competition between MDCK^WT and MDCK^Scrib as well as the upregulation of MDCK^WT division in MDCK^Scrib-dominated neighbourhoods (Bove et al., 2017), suggesting that mechanical competition represents the dominant mechanism of population change in these experiments.

To understand the mechanistic origin of density-mediated cell competition, we varied the growth rate G, the intercellular adhesion J_heterotypic, the stiffness λ, and the contact inhibition parameter k for the individual cell types starting from an initial set of values that gave rise to mechanical competition (Supplementary file 2). We reasoned that, in a competition setting, the values of each parameter in one cell type relative to the other were likely more important than their absolute values. Therefore, we varied each parameter in only one of the two cell types.

We first varied the growth rate G of loser cells (Supplementary file 2). G controlled the time required for elimination and the peak loser cell count but did not affect the outcome of competition with losers being eliminated for all growth rates examined (Figure 5—figure supplement 1A, B), consistent with simulations of population dynamics based on ordinary differential equations (Basan et al., 2009).

Second, we varied the heterotypic adhesion energy, J_heterotypic, between the winners and the losers between ±50% of the value of the homotypic adhesion J_homotypic (Supplementary file 2). When J_heterotypic is larger than J_homotypic, cells preferentially adhere to cells of their own type. In our simulations, J_heterotypic did not appear to change the kinetics or the outcome of competition (Figure 5—figure supplement 1C). Note that varying J_homotypic in one of the cell types only would have similar effects to a variation in J_heterotypic.

In our simulations, sensitivity to contact inhibition k was chosen to be the same for both cell types. This parameter constrains how far cells can deviate from their target area A_T before their growth rate G(t) approaches 0 and they stop growing (Materials and methods, Supplementary file 2). In pure winner cell populations, the average local density reached a plateau after confluence defining a homeostatic density ($HD$), which decreased with increasing contact inhibition k (Figure 5A, blue line). However, this effect was not observed in pure loser populations (Figure 5A, red line) because their probability of apoptosis is high even for densities below the homeostatic density dictated by $k$ (Figure 1—figure supplement 4A). Indeed, under normal growth conditions, we predict that loser cells never reach densities where contact inhibition parametrised by k becomes active. In all cases, the homeostatic density of winner cells was higher than in loser cells but the difference in homeostatic density, $∆HD$, decreased with increasing $k$ (Figure 5A). Thus, in winner cells, homeostatic density is controlled by a decrease in growth controlled by the contact inhibition parameter $k$, while in loser cells it is controlled by density-dependent apoptosis (Figure 1—figure supplement 4A).

Figure 5 with 1 supplement see all

Download asset Open asset

In competitions, when we varied the homeostatic density of the winner cells (by changing $k$, Supplementary file 2), we found that, after 80 hr, loser cells were completely eliminated for high values of $∆HD$ but they survived when $∆HD$ was lower (Figure 5B, C, E). In addition, the time required for elimination of 50% of loser cells increased with decreasing $∆HD$ and increasing k (Figure 5A, F). Interestingly, loser cell count appeared to converge towards a non-zero plateau for values of k larger than 0.1 at long time scales (Figure 5—figure supplement 1D). Therefore, the difference in homeostatic density, $∆HD$, between the winner and loser cells governs the kinetics and the outcome of mechanical competition for the durations examined in this study (Figure 5A, C ,F, Figure 5—figure supplement 1D).

As in our initial parameterisation the winner cells have a higher stiffness $\lambda$, the loser cells are compressed by the winners during competition. As a result, the average local density of loser cells is larger than that of winners, which, combined with their greater sensitivity to crowding, leads to increased apoptosis. To determine the impact of $\lambda$ on competition, we varied the loser cell stiffness while maintaining winner cell stiffness constant. When the relative stiffness parameter $\Lambda ={\lambda }_{loser}/{\lambda }_{winner}$ was smaller than 1, loser cells were eliminated (Figure 5D, Figure 5—figure supplement 1E, F). By contrast, when $\Lambda$ was equal to or higher than 1, loser cells survived (Figure 5D, Figure 5—figure supplement 1E, G). Akin to $∆HD$, changes in the ratio of winner-to-loser cell stiffness altered the kinetics of competition and its outcome over durations considered in this study (Figure 5E, Figure 5—figure supplement 1E). These results are consistent with experiments showing that competition is decreased in the presence of an inhibitor of Rho-kinase (Y27632) that reduces cell contractility in both populations (Wagstaff et al., 2016).

When $∆HD$ was low or $\Lambda$ was larger than 1, the change in competition outcome occurred because of a decrease in the local density of loser cells in mixed populations, which in turn led to decreased apoptosis. However, winners have an extra competitive edge because when free space becomes available due to cell death or cell area compressibility, they take advantage of the free space due to faster growth using a squeeze and take or a kill and take tactic (Gradeci et al., 2020). At very long times, this effect alone may be sufficient for them to dominate in mixed populations.

Our simulation can also be used to gain mechanistic insights into biochemical competition. Recent work has shown that, during biochemical competition, apoptosis in loser cells is governed by the extent of their contact with winner cells (Levayer et al., 2016) and that perturbations that increase mixing between cell types increase competition (Levayer et al., 2015).

To study biochemical competition in isolation from any mechanical effect, we assumed that both cell types have identical stiffnesses λ, equal sensitivities to contact inhibition k, and high but equal homeostatic densities (Supplementary file 3). In both cell types, we modelled the dependency of apoptosis on the proportion of cell perimeter $p$ engaged in heterotypic contact by using a Hill function parameterised by a steepness $S$ and an amplitude ${\displaystyle p_{apomax}}$:

${\displaystyle p_{apo}\left(p\right)=\frac{p_{apomax}{p}^n}{\left(S^n+p^n\right)}}$, where $n$ is the Hill coefficient (Figure 1—figure supplement 5A).

When $S$ decreases, the probability of apoptosis increases rapidly with the extent of heterotypic contact (Figure 1—figure supplement 5A). For winner cells, we chose a low ${\displaystyle p_{apomax}}$ and high $S$ because we do not expect their apoptosis to show sensitivity to contact with loser cells. In contrast, for loser cells, we chose ${\displaystyle p_{apomax}}$ to be 10-fold higher than in winners, giving an amplitude similar to the maximal probability of apoptosis observed in losers during mechanical competition and similar kinetics of elimination, as observed in experiments (Hogan et al., 2009; Figure 1—figure supplement 4A, Figure 1—figure supplement 5A, B). To investigate biochemical competition, we varied parameters modulating contact between cells (the heterotypic adhesion J_heterotypic), apoptosis of losers (the Hill function parameter S of loser cells), as well as tissue organisation.

First, we assumed a homogenous seeding of each cell type with a 50:50 ratio between winners and losers. Competition depends on the relative probability of apoptosis p_apo in winners and losers as a function of the fraction of their perimeter p in heterotypic contact (Figure 1—figure supplement 1C, right). Therefore, varying ${\displaystyle p_{apomax}}$ or S has a qualitatively similar overall effect on competition. We examined the dependency of competition outcome on $S_{Loser}$ with S_Winner, ${\displaystyle p_{apomax,winner}}$, and ${\displaystyle p_{apomax,loser}}$ fixed (Supplementary file 3). For all values of $S_{Loser}$, loser cells first increased in number until overall confluence at ~60 hr, before decreasing after that (Figure 1—figure supplement 5B). When $S_{Loser}$ was low, losers were eliminated because p_apo for losers was higher than for winners for all heterotypic contact extents (Figure 1—figure supplement 5B, C, Video 6), whereas when $S_{Loser}$ was high, winners and losers had comparable p_apo when in heterotypic contact, leading to coexistence because no competition took place (Figure 1—figure supplement 5B, C).

Video 6

Download asset

Next, we investigated the dependency of cell competition on initial seeding conditions for values of S_Loser that gave rise to competition (S_Loser = 0.3, S_Winner = 0.5, Supplementary file 3). We examined three different initial seeding conditions: fully mixed (Figure 6A, middle column), partially sorted with loser cells grouped into a few colonies (Figure 6B, middle column, Videos 7–9), and fully sorted with loser cells and winner cells occupying opposite sides of the field of view (Figure 6C, middle column, Video 10). Strikingly, in mechanical competition, the normalised count of loser cells reached a maximum around 60 hr before continuously decreasing thereafter, consistent with our experimental observations in fully mixed (Figure 4B) and partially sorted conditions (Figure 6—figure supplement 2A–C, Video 9). By 160 hr, loser cells had been eliminated for all configurations (Figure 6A–C, left-hand column, E). In contrast, in biochemical competition, the outcome of competition appeared strongly dependent on initial seeding conditions with large differences in normalised count of losers after 100 hr (Figure 6A–C, right-hand column, D). Indeed, loser cell normalised count was close to 0 for fully mixed seeding but remained larger than 1 for partially and fully sorted seedings. By 200 hr, loser cell count had dropped to 0 in partially sorted seeding but only decreased gradually for fully sorted seeding (Video 10). This suggests that the kinetics of biochemical competition is sensitive to tissue organisation. To quantitatively compare tissue organisations, we computed the evolution of mixing entropy, a measure of local tissue organisation, in each competition (Materials and methods). We found that, when cells reached confluence, entropy of mixing was highest in the fully mixed seeding and lowest in the fully sorted seeding (Figure 6—figure supplement 1A). In the fully mixed and partially sorted seedings, mixing entropy decreased after overall confluence as the competition progressed, whereas for fully sorted seeding, mixing entropy stayed constant because the interface between winners and losers maintained its shape over time even though the number of loser cells gradually decreased (Figure 6C, Figure 6—figure supplement 1A, Video 10). Loser colony size and geometry may therefore determine the kinetics and outcome of biochemical competition. To gain further insight, we systematically varied the size of the loser colony and determined the time to elimination (Figure 6F). This revealed that the time to elimination monotonously increased with cell number in the colony but that no true steady-state coexistence was reached.

Figure 6 with 2 supplements see all

Download asset Open asset

Video 7

Download asset

Video 8

Download asset

Video 9

Download asset

Video 10

Download asset

As we found that the kinetics of biochemical competition was controlled by tissue organisation and the intermixing of cells, we examined how the relative magnitude of heterotypic versus homotypic intercellular adhesion energy affects competition. In our simulations, we varied the heterotypic adhesion, $J_{heterotypic}$, while keeping J_homotypic constant. In the Potts model, J represents a surface energy which is the difference between the surface tension and adhesion. Therefore, when J is high, intercellular adhesion is low. We found that the normalised loser cell count after 100 hr decreased with decreasing $J_{heterotypic}$ in fully mixed and partially sorted tissue organisation (Figure 6—figure supplement 1C, D). This is because for low $J_{heterotypic}$, cell intermixing is favoured and winner cells can invade colonies of loser cells, consistent with experimental observations (Levayer et al., 2015). However, in the fully sorted configuration, where the interface between the two cell types is minimal (Figure 6C), changes in $J_{heterotypic}$ only have a weak effect (Figure 6—figure supplement 1B). Thus, the kinetics of biochemical competition is sensitive to changes in parameters that lead to tissue reorganisation, such as the relative magnitude of homotypic and heterotypic adhesion energy.

In summary, our simple implementation of biochemical competition was sufficient to qualitatively reproduce current experimental observations, although the precise experimental curves relating probability of apoptosis to extent of heterotypic contact remain to be accurately determined experimentally.

The CPM is implemented in Compucell3D (Swat et al., 2012). We chose a 2D lattice-based model, where cells are composed of a collection of lattice sites (pixels). Cells interact at their interfaces through predefined adhesion energies, and several different cell types can be implemented. Each cell is then given attributes characterising their mechanical and adhesive properties. For example, each cell is assigned a cell type $\tau$, which in turn has some value of surface contact energy $J_{cell-cell}$ with other cell types and adhesion energy $J_{cell-substrate}$ with the substrate. Cells are also assigned a target area $A_T\left(t\right)$ (that represents the area a cell would occupy at time t if it were isolated) and an area expansion modulus $\lambda$ (that represents the energetic cost of increasing cell area and originates from the mechanical properties of the cytoskeleton). In the computational cell decision-making in our simulation, $A_T$ and $\lambda$ play important roles in the implementation of growth and division dynamics. In addition, we also incorporate active cell motility. The free energy of the system is given by the Hamiltonian $H$:

where the first term describes the interaction of lattice sites due to the adhesion energy between the cell types or between cells and the substrate. The coefficient $J$ is the surface energy between cell type $\tau$ of the target lattice site ${\sigma }^k$ and the cell type $\tau \left({\sigma }_{i^{\text{'}}{j}^{\text{'}}}\right)$ of its nearest-neighbour lattice points. By convention, for free space τ = 0. From a biological perspective, $J$ is the difference between the surface tension and intercellular adhesion. Therefore, higher $J$ implies lower intercellular adhesion. The multiplicative term ${\displaystyle \left(1-\delta \left({\sigma }_{ij}^k,{\sigma }_{i^{\prime }{j}^{\prime }}\right)\right)}$ prevents cells from interacting energetically with themselves, where

The second term in the Hamiltonian describes an additional energy cost due to deviation of the actual area $A\left(\sigma \right)$ of a cell from its target area $A_T\left(\sigma \right)$, specific to each cell at time t. In the second computational layer of the simulation, A_T is varied at each time point to reflect cell growth. The coefficient $\lambda$ represents the area expansion modulus in 2D, which is related to planar cell stiffness or the ratio between apical and lateral contractility that control cell height. We introduce the term $\Theta \left(\tau \right)$ to treat the free space pixels differently from pixels belonging to cells. In contrast to cells, the free space does not have a target area, and hence no associated mechanical energy.

The final term in the Hamiltonian assigns active motility to the cells along a random unit vector $\widehat{m}$ (Li and Lowengrub, 2014). Here, $\widehat{s}$ is the spin flip direction between the lattice site in question and one of its neighbouring lattice sites.

To describe epithelial cell dynamics using the Potts model, we parametrised it using our experimental data (Bove et al., 2017). For simplicity, we chose the same length scale for pixels in our simulation as in our experimental images of competition experiments. The lattice size and cell sizes are chosen to match the experimental data. The lattice is chosen to be 1200 × 1600 pixels, where each pixel is 0.33 × 0.33 μm².

The target area and stiffness of each cell type were determined based on the average cell areas measured from cells isolated from one another in brightfield images (Figure 2—figure supplement 2).

In the simulation, one Monte Carlo timestep (MCS) is defined by each lattice point being given the possibility of changing identity. A conversion between experimental time and computational time was derived empirically by comparing the mean squared displacements of isolated cells in experiments to those in the simulations. We found that 10 MCS represented one frame of a time-lapse movie in our experiments (4 min).

The agent-based part of the model requires the introduction of cellular behaviour in the form of probabilistic rules for cell growth, division, extrusion, and apoptosis. In our simulations, cells grow linearly by increasing their target areas A_T(t) at a rate G, which was chosen to replicate the average cell doubling time measured in experiments pre-confluence (Bove et al., 2017). In line with recent experimental work (Cadart et al., 2018), we assume that MDCK cells follow an ‘adder’ mechanism for cell size control, such that cells divide along their major axis once a threshold volume ΔA_tot has been added since birth (Figure 1—figure supplement 1, Figure 1—figure supplement 3). In our simulations, the added cell volume at each time point was a random value distributed around the mean experimental value G, so as to capture cell-to-cell variability. When the simulation was initialised, cell areas had a homogenous distribution to mimic a uniform probability for cells of being in any given stage of their cell cycle at the start of experiments.

In our simulations, cells possess a target area, $A_T\left(t\right)$, which they would occupy at that time if they had no neighbours, and an actual area, $A\left(t\right)$, which they currently occupy. As A_T increases at each timestep due to cell growth, the difference between their target and actual area A increases. If this difference becomes too large, the second term of the Hamiltonian dominates, leading to energetically unfavourable swaps and a collapse of the network. To mimic reduced protein synthesis reported due to contact inhibition of proliferation (Azar et al., 2010), we assume that the effective growth rate depends on the difference between A and A_T:

where $G$ is the growth rate for cells with no neighbours, $A_T$ the target cell area, $A$ the actual area, and $k$ quantifies the sensitivity to contact inhibition. k parametrises how much deviation can be tolerated between the target area and current cell area before growth stalls. Note that this condition is applied iteratively at every frame for each cell, such that, when free space becomes available, growth can immediately resume nearby.

In crowded conditions such as those present in mechanical competition, the probability of apoptosis increases with local cell density (Wagstaff et al., 2016; Bove et al., 2017; Eisenhoffer et al., 2012). To implement this, each cell was assigned a probability of apoptosis p_apo at each timestep that depended on its local cell density ρ. In our simulation, ρ was defined as the sum of inverse of areas of the cell of interest σ^k and its first neighbours σⁱ: $\rho \left({\sigma }^k\right)=\frac{1}{A\left({\sigma }^k\right)}+{\sum }_{i=1}^n\frac{1}{A\left({\sigma }^i\right)}$ (Figure 1—figure supplement 4B). Based on our experimental data, we decided to describe the relationship between p_apo and ρ as: $p_{apo}\left(\rho \right)=\frac{p_{apo,max}}{\left(1+e^{-\alpha \left(\rho -{\rho }_{1/2}\right)}\right)}$ (Figure 1—figure supplement 4A). p_apo,max was fixed to the same value for both populations, and ρ_1/2 was determined for each population separately based on experimental data (Bove et al., 2017).

In biochemical competition, apoptosis occurs when loser cells are in direct contact with winner cells. Recent work has shown that in Drosophila the probability of apoptosis of loser cells depends on the percentage of the perimeter in contact with the winner cells (Levayer et al., 2015). Following this, we chose to implement the probability of apoptosis as a sigmoid function (Hill function) following the relationship $p_{apo}\left(p\right)=p_{apo,max}{p}^n/(S^n+p^n)$, where $p$ is the percentage of perimeter in heterotypic contact, $n$ is the Hill coefficient, $p_{apo,max}$ is the maximum probability, and $S$ is the steepness. We chose a maximum probability $p_{apo,max}$ of death per frame similar to that encountered in mechanical competition and a Hill coefficient $n=3$ (Figure 1—figure supplement 5A). This is justified by the fact that mechanical and biochemical competition take place over comparable durations in MDCK cells ~2–4 days (Hogan et al., 2009; Kajita et al., 2010). For both mechanical and biochemical competition, the execution of apoptosis was implemented by setting the target area A_T of the cell to 0 and area expansion modulus λ to 2. This allows for a quick but not instantaneous decrease of the cell area until the cell is completely removed.

Under conditions where the local cell density increases rapidly, live cells can be extruded from monolayers, likely because they have insufficient adhesion with the substrate to remain in the tissue (Eisenhoffer et al., 2012; Kocgozlu et al., 2016; Figure 1—figure supplement 1D). We assumed that cells underwent live extrusion when their area A_i was $Ai\le ⟨A\left(t\right)⟩/2$ with <A(t)> the average area of all cells in the simulation. Once this occurs, the cell is eliminated immediately from the tissue. Unlike apoptosis, cell elimination via extrusion is implemented as an instantaneous removal of the qualifying cell from the lattice to reflect the faster rate of live extrusions compared to programmed cell death.

All simulated data for mechanical competition were compared quantitatively to experiments acquired in Bove et al., 2017 or performed specifically for this publication. Methods for cell culture, image acquisition, segmentation, and analysis are described in detail in Bove et al., 2017. Briefly, MDCK wild-type cells (MDCK^WT) were winners in these competitions and their nuclei were labelled with H2B-GFP, while MDCK scribble knock-down cells (MDCK^Scrib, described in Norman et al., 2012) were the losers and labelled with H2B-RFP (Bove et al., 2017). MDCK^Scrib cells conditionally expressed an shRNA targeting scribble that could be induced by addition of doxycycline to the culture medium. All cell lines were regularly tested for mycoplasma infection and were found to be negative (MycoAlert Plus Detection Kit, Lonza, LT07-710).

MDCK^WT cells were grown in DMEM (Thermo Fisher) supplemented with 10% fetal bovine serum (Sigma-Aldrich), HEPES buffer (Sigma-Aldrich), and 1% penicillin/streptomycin in a humidified incubator at 37°C with 5% CO₂. MDCK^Scrib cells were cultured as MDCK^WT, except that we included tetracycline-free bovine serum (Clontech, 631106) to supplement the culture medium. To induce expression of scribble shRNA, doxycycline (Sigma-Aldrich, D9891) was added to the medium at a final concentration of 1 µg/ml.

For competition assays, winner and loser cells were seeded in the chosen proportion to reach an overall density of 0.07 cells per 100 μm² and left to adhere for 2 hr. Cells were then imaged every 4 min for 4 days using a custom-built incubator microscope and the appropriate wavelengths (Bove et al., 2017). Movies were then automatically analysed to track the position, state, and lineage of the cells using deep-learning-based image classification and single-cell tracking as detailed in Bove et al., 2017. Single-cell tracking was performed using bTrack (github.com/quantumjot/BayesianTracker; Lowe, 2021). Cell neighbours are determined using a Voronoi tessellation dual to the Delaunay triangulation of nuclei.

In most experiments, knock-down of scribble was induced by addition of doxycycline 48 hr before the beginning of the experiment to ensure complete depletion. However, for experiments examining cell competition in partially sorted conditions, we first seeded the same number of MDCK^Scrib cells as in 90:10 competitions and cultured them for 48 hr without doxycycline until they formed well-defined colonies. Then, we added the number of MDCK^WT cells that would be expected after 48 hr competition, left them to adhere for 4 hr, and added doxycycline. We then started imaging the following day, signifying that scribble depletion was incomplete at the start of the experiment.

For the analysis of experiments and simulations, fate information for each cell is dynamically recorded to a file and analysed using a custom software written in Matlab (Bove et al., 2017).

The entropy of mixing was calculated as the Shannon entropy of a two-state system, where the states considered are the cell types (winner/loser). The entropy was then calculated as $s=-P_1\ln P_1-P_2\ln P_2$ for each cell at each frame, where ${\displaystyle P_1=\frac{\mathrm{\#}winnerneighbors}{total\mathrm{\#}neighbors}}$ and ${\displaystyle P_2=\frac{\mathrm{\#}loserneighbors}{total\mathrm{\#}neighbors}}$. The entropy of the whole system was then calculated as $S=<s>/\sum cells$.

To calculate the probability of apoptosis and division for each cell type, cells were binned appropriately (by density or by time). Then, we determined the number of events n (apoptosis or division) and the total number of cells of each type N (the observations) in each bin. The probability p was then computed as $p=\frac{n}{N}$. Because we are examining rare events, we calculated the coefficient of variation cv that measures the relative precision of our estimator of probability as $cv=\sqrt{\frac{(1-p)}{pn}}$.

Our model has been deposited in Github (https://github.com/DGradeci/cell_competition_paper_models; copy archived at swh:1:rev:55f8b189c6f5d998cc5b2819f672ad80b547c956; Gradeci, 2021). The data used for model calibration will be deposited in doi: 10.5522/04/12287465.

The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.

This work was supported by Engineering and Physical Sciences Research Council (EPSRC) PhD studentships to DG and AB. SB acknowledges support from the Royal Society (URF/R1/180187). GV was supported by BBSRC grant BB/S009329/1 to AL and GC. The authors also wish to acknowledge the reviewers whose insightful comments greatly improved the study.

1. Received: July 13, 2020
2. Accepted: April 23, 2021
3. Version of Record published: May 20, 2021 (version 1)

© 2021, Gradeci et al.

This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.