Mechanical basis and topological routes to cell elimination

Version of Record

Accepted for publication after peer review and revision.

Download
Cite
Share
CommentOpen annotations (there are currently 0 annotations on this page).

Version of Record published: April 18, 2023 (This version)
Accepted: March 22, 2023
Received: August 3, 2022
Preprint posted: August 18, 2021 (Go to version)

1. Of interest
Morphogenesis: The enigma of cell intercalation

Raphaël Clément

Insight May 3, 2024
Further reading

Abstract
Editor's evaluation
Introduction
Results and discussion
Materials and methods
Appendix 1
Data availability
References
Article and author information
Metrics

Abstract

Cell layers eliminate unwanted cells through the extrusion process, which underlines healthy versus flawed tissue behaviors. Although several biochemical pathways have been identified, the underlying mechanical basis including the forces involved in cellular extrusion remains largely unexplored. Utilizing a phase-field model of a three-dimensional cell layer, we study the interplay of cell extrusion with cell–cell and cell–substrate interactions in a flat monolayer. Independent tuning of cell–cell versus cell–substrate adhesion forces reveals that extrusion events can be distinctly linked to defects in nematic and hexatic orders associated with cellular arrangements. Specifically, we show that by increasing relative cell–cell adhesion forces the cell monolayer can switch between the collective tendency towards fivefold, hexatic, disclinations relative to half-integer, nematic, defects for extruding a cell. We unify our findings by accessing three-dimensional mechanical stress fields to show that an extrusion event acts as a mechanism to relieve localized stress concentration.

Editor's evaluation

In this work, Monfared et al. construct a valuable three-dimensional phase-field model for cell monolayers and use this to investigate the relationship between single-cell extrusion events and topological defects in cellular arrangement. The extension of existing 2D phase field models to three dimensions is an important contribution of this paper, which will be of general interest to the theoretical modelling of epithelial monolayers. Here the model is used to study the importance of cell-cell and cell-substrate interaction in extrusion from cell monolayers, which will be of practical interest to biologists and physicists working on this process. This paper presents convincing evidence that extrusion events are distinctly linked to defects in nematic and hexatic orders in the cell monolayer.

https://doi.org/10.7554/eLife.82435.sa0

Introduction

The ability of cells to self-organize and collectively migrate drives numerous physiological processes including morphogenesis (Chiou and Collins, 2018; Vafa and Mahadevan, 2022), epithelial–mesenchymal transition (Barriga et al., 2018), wound healing (Brugués et al., 2014), and tumor progression (De Pascalis and Etienne-Manneville, 2017). Advanced experimental techniques have linked this ability to mechanical interactions between cells (Maskarinec et al., 2009; Ladoux, 2009; Ladoux and Mège, 2017). Specifically, cells actively coordinate their movements through mechanosensitive adhesion complexes at the cell–substrate interface and cell–cell junctions. Moreover, cell–cell and cell–substrate adhesions seem to be coupled (Balasubramaniam et al., 2021), further complicating the interplay of mechanics with biochemistry.

While advances in experimental techniques are followed by more nuanced theoretical and computational developments, a majority of current approaches to simulate multicellular layers are limited to two-dimensional systems, hindering in-depth exploration of intrinsically three-dimensional nature of the distinct forces that govern cell–cell and cell–substrate interactions. Furthermore, some of the most fundamental processes in cell biology such as cell extrusion – responsible for tissue integrity – are inherently three-dimensional. Thus, studying the underlying mechanisms necessitates access to both in-plane and out-of-plane forces in the cell layers.

Cell extrusion refers to the process of removal of excess cells to prevent accumulation of unnecessary or pathological cells (Rosenblatt et al., 2001). This process can get initiated through apoptotic signaling (Rosenblatt et al., 2001), oncogenic transformation (Hogan et al., 2009), and overcrowding of cells (Marinari et al., 2012; Eisenhoffer et al., 2012; Levayer et al., 2016) or induced by replication stress (Dwivedi et al., 2021). Most importantly, cell extrusion plays an important role in developmental (Toyama et al., 2008), homeostatic (Eisenhoffer et al., 2012; Le et al., 2021), and pathological processes (Slattum and Rosenblatt, 2014), including cancer metastasis. However, the underlying mechanisms that facilitate cell extrusion are still unclear.

The similarities between cellular systems and liquid crystals, studied both theoretically and experimentally, featuring both nematic order (Saw et al., 2017; Kawaguchi et al., 2017; Duclos et al., 2018; Blanch-Mercader et al., 2018; Tan et al., 2020; Zhang et al., 2021) and hexatic order (Classen et al., 2005; Sugimura and Ishihara, 2013; Pasupalak et al., 2020; Maitra et al., 2020; Hoffmann et al., 2022) with the two phases potentially coexisting (Armengol-Collado et al., 2022) and interacting provide a fresh perspective for understanding cellular processes. The fivefold disclinations in hexatic arrangement of cells are numerically shown to favor overlaps between the cells in two-dimensions (Loewe et al., 2020), potentially contributing to the cell extrusion in three-dimensions. In this vein, it is shown that a net positive charge associated with hexatic disclinations can be associated with the maximum curvature of dome-like structures in model organoids and in epithelial cell layers (Rozman et al., 2020; Rozman et al., 2021; Hoffmann et al., 2022). Moreover, in cellular monolayers, comet- and trefoil-shaped half-integer topological defects, corresponding to +1/2 and -1/2 charges, respectively, are prevalent (Doostmohammadi et al., 2015; Doostmohammadi et al., 2016). These are singular points in cellular alignment that mark the breakdown of orientational order (de Genne and Prost, 1998). Recent experiments on epithelial monolayers found a strong correlation between extrusion events and the position of a subset of +1/2 defects in addition to a relatively weaker correlation with -1/2 defects (Saw et al., 2017). These recently introduced purely mechanical routes to cell extrusion have opened the door to new questions on the nature of forces that are involved in eliminating cells from the monolayer and challenge the purely biological consensus that an extruding cell sends a signal to its neighbor that activates its elimination process (Rosenblatt et al., 2001). Nevertheless, it is not clear whether these different mechanisms are related, and whether, depending on the mechanical features of the cells, the cell layers actively switch between different routes to eliminate the unwanted cells. Since all the existing studies so far have only focused on effective two-dimensional models of the cell layers, fundamental questions about the three-dimensional phenomenon of cell extrusion and its connection to the interplay between cell-generated forces at the interface between cells and the substrate, with multicellular force transmission across the cell layer, remain unanswered.

In this article, we explore three-dimensional collective cell migration in cellular monolayers. Based on large-scale simulations, we examine (i) the underlying mechanisms responsible for cell extrusion, including any correlations with ±1/2 topological defects and fivefold disclinations, and (ii) the interplay of cell–cell and cell–substrate adhesion with extrusion events in cellular systems. Moreover, by mapping the full three-dimensional mechanical stress field across the entire monolayer, we identify localized stress concentration as the unifying factor that governs distinct topological routes to cell extrusion.

Results and discussion

Topological routes to cell extrusion: Nematic and hexatic disclinations

In the absence of self-propulsion forces, the initial configuration tends to equilibrate into a hexagonal lattice (see Appendix 1—figure 11 in Appendix 1 for an example). As we introduce self-propulsion forces associated with front-rear cell polarity (see ‘Materials and methods’ for polarization dynamics), the system is pushed away from its equilibrium hexagonal configuration, resulting in defects manifested as fivefold and sevenfold disclinations, as shown in Figure 2b. Figure 1a shows a simulation snapshot with two extrusion events taking place. An extrusion event is detected if the vertical displacement of a cell, relative to other cells in the monolayer, exceeds $R_{0} / 2$ , where R₀ is the initial cell radius. Figure 1b displays the out-of-plane normalized velocity profile, ${\vec{\tilde{v}}}_{z} = (\vec{v} (\vec{x}) \cdot {\vec{e}}_{z}) / v_{z}^{max}$ where $v_{z}^{max}$ is the maximum value of the velocity component in ${\vec{e}}_{z}$ direction in the displayed cross-section of the monolayer, clearly marking the extruding cells as they get expelled from the monolayer and lose contact with the substrate.

Figure 1

Download asset Open asset

Cell extrusion in a 3D representation of a confluent cell layer.

(a) A representative simulation snapshot (cell–substrate adhesion $ω_{cw} = 0.0025$ and relative cell–cell adhesion $Ω = ω_{cc} / ω_{cw} = 0.4$ ) of a three-dimensional cell monolayer. Two cells are visibly extruding. (b) A cross-section (dotted yellow line $ᵉ - ᵉ^{^{'}}$ ) of the cell monolayer highlighting the two extruding cells via the normalized out-of-plane velocity ( ${\tilde{v}}_{z} = (\vec{v} \cdot {\vec{e}}_{z}) / v_{z}^{max}$ ), where $v_{z}^{max}$ is the maximum value of the v_z component of the velocity field $\vec{v}$ in the shown cross-section.

In order to probe the possible mechanical routes to cell extrusion, we begin by characterizing topological defects in cell orientation field and disclinations in cellular arrangements. To this end, we first map the orientation field of the cells from the 2D projected cell shape profile on $x y -$ plane ( $z = 0$ , i.e., the basal side) and identify topological defects as the singularities in the orientation field. The results (example snapshot in Figure 2a) show the continuous emergence of half-integer (±1/2), nematic, topological defects that spontaneously nucleate in pairs and follow chaotic trajectories before annihilation (see Appendix 1—figure 9 in Appendix 1 for energy spectra characterization). It is noteworthy that unlike previous studies of active nematic behavior in 2D cell layers (Mueller et al., 2019; Wenzel and Voigt, 2021), the nematic defects here emerge in the absence of any active dipolar stress or subcellular fields as the only active driving in these simulations is the polar force that the cells generate. Therefore, although the cells are endowed with polarity in terms of their self-propulsion, the emergent symmetry in terms of their orientational alignment is nematic, which is in line with experimental observations in cell monolayers (Saw et al., 2017; Blanch-Mercader et al., 2018), discrete models of self-propelled rods (Bär et al., 2020; Meacock et al., 2021), and recently proposed continuum model of polar active matter (Amiri et al., 2022).

Figure 2

Download asset Open asset

Nematic and hexatic disclinations govern cell extrusion.

A representative analysis corresponding to the configuration shown in Figure 1a and projected into $x y -$ plane ( $z = 0$ , i.e., the basal side). (a) A coarse-grained director field with coarse-graining length of one cell size $ℓ_{dir.} = R_{0}$ and +1/2 (filled circles with the line indicating orientation) and -1/2 (three connected lines with threefold symmetry) nematic defects. (b) Number of neighbors $z$ for each cell, including fivefold and sevenfold disclinations mapped into the monolayer. The symbol + denotes the center of mass for two extruding cells. (**c, d**) Probability densities of the normalized minimum distance between extruding cells and the nearest ±1/2 defect, ${\tilde{d}}_{min} = d_{min} / R_{0}$ , for varying cell–cell to cell–substrate adhesion ratios $Ω$ for (c) -1/2 and (d) +1/2 topological defects (inset: distribution mean $m = ⟨ {\tilde{d}}_{min}^{+ 1 / 2} ⟩$ vs. $Ω$ ). (e) The probability density of average coordination number $\bar{z}$ for an extruding cell during $\tilde{t} = (t / τ_{0}) \in [{\tilde{t}}_{e} - 2.5, {\tilde{t}}_{e} + 0.3125]$ , where ${\tilde{t}}_{e}$ denotes extrusion time, $τ_{0} = ξ R_{0} / α$ and for varying cell–substrate to cell–cell adhesion ratios $Ω$ . The data in (c–e) corresponds to four different realizations.

Remarkably, in accordance with experimental observations (Saw et al., 2017), we find that the extrusion events can be correlated with the position of both +1/2 comet-shaped and -1/2 trefoil-shaped topological defects. To quantify this, Figure 2c and d display the probability density of the normalized minimum distance ${\tilde{d}}_{min}^{\pm 1 / 2} = d_{min}^{\pm 1 / 2} / R_{0}$ between an extruding cell and ±1/2 topological defects in the interval $\tilde{t} \in [{\tilde{t}}_{e} - 5.625, {\tilde{t}}_{e} + 0.625]$ , where $\tilde{t} = t / τ_{0}$ is the normalized time, $τ_{0} = ξ R_{0} / α$ , $ξ$ corresponds to cell–substrate friction, $α$ denotes the strength of polarity force, and ${\tilde{t}}_{e}$ is the (normalized) extrusion time. This temporal window is chosen based on the first moment of a defect’s lifetime distribution (see Appendix 1—figure 5 in Appendix 1). The data in Figure 2c and d is based on four distinct realizations and for varying cell–substrate to cell–cell adhesion ratios, $Ω = ω_{cc} / ω_{cw}$ . For both defect types, the probability density peaks in the vicinity of the eliminated cell ( $\approx 1.5 R_{0}$ ), at a much smaller distance relative to a typical distance between two defects (see Appendix 1—figure 7 in Appendix 1), and falls off to nearly zero for $d_{min}^{\pm 1 / 2} 5 R_{0} (= 40)$ . Furthermore, laser ablation experiments have established that an induced extrusion event does not favor the nucleation of a pair of nematic defects (Saw et al., 2017).

In a hypothesis-testing approach, we check whether these peaks in the minimum distance represent a correlation between extrusion events and nematic defects. To this end, we set out to falsify the hypothesis that the extrusion events are uncorrelated with the nematic defects. We utilize a Poisson point process to randomly generate positions for extrusion events and quantify the minimum distance between each event and the nearest half-integer nematic defect. For each simulation, we generate five different realizations for the extrusion events using a Poisson point process with the intensity set equal to the number of extrusions in that particular simulation. The extrusion time is also a random variable described by a uniform distribution, $t_{e} \sim U (1, n_{sim})$ , where $n_{sim} = 29, 000$ is total number of time steps. As an example, Figure 3a shows probability density function for ${\tilde{d}}_{min}^{+ 1 / 2}$ for $Ω = 0.4$ using simulation data as well as data randomly generated with Poisson point process. Finally, a Kolmogorov–Smirnov (KS) test is used to measure if the two samples, one based on our simulations and one based on randomly generated extrusion events, belong to the same distribution. The results of the KS test reject this (see Appendix 1—table 1 and Appendix 1—figure 1 in Appendix 1) and thus falsify the hypothesis that simulation-based extrusion events are uncorrelated with the half-integer nematic defects.

Figure 3

Download asset Open asset

Topological, rather than geometrical, route to cell extrusion.

(a) Probability density functions for normalized minimum distance between an extrusion event and a +1/2 defect, ${\tilde{d}}_{min}^{+ 1 / 2}$ , based on simulation results and randomly generated through a Poisson point process and for $Ω = 0.4$ . (b) Comparison of Lewis’s linear and quadratic relations with our simulations. ${\bar{A}}_{z}$ is the average area for cells with $z$ neighbors and $\bar{A}$ is the average area of all cells. The color bar indicates simulation time step, and the data correspond to the case $ω_{cw} = 0.0025$ and $Ω = 0.4$ .

Next, we explore the other possible mechanical route to cell extrusion based on the disclinations in cellular arrangement. To this end, we compute the coordination number of each cell based on their phase-field interactions and identify the fivefold and sevenfold disclinations (see Figure 2b). To quantify the relation between extrusion events and the disclinations, the probability density of the coordination number of an extruding cell ( ${\tilde{d}}_{min} = 0$ ) averaged over the time interval, $\tilde{t} \in [{\tilde{t}}_{e} - 5.625, {\tilde{t}}_{e} + 0.625]$ , $\bar{z}$ , for all the realizations is shown in Figure 2e, clearly exhibiting a sharp peak near $\bar{z} = 5$ . The coordination number is determined based on the interactions of cells (see Appendix 1) and this property is independent of apical or basal considerations (Kaliman et al., 2021), unlike geometrical structures called scutoids that have been identified in curved epithelial tubes (Gómez-Gálvez et al., 2018). In our setup, the asymmetric interactions of cells with apical and basal sides are captured by varying the strength of cell–substrate adhesion. In our simulations, increasing cell–substrate adhesion leads to lower extrusion events (see Appendix 1—figure 8 in Appendix 1).

Thus far, our results suggest topological rather than geometrical routes to cell extrusion. To probe the role of geometrical constraints further, we investigate the existence of any correlation between cell area and its number of neighbors. The best known such a correlation – for cellular matter with no gaps between them, that is, confluent state – is a linear one and it is due to Lewis, 1928 with other types of relations, for example, quadratic, proposed since his work (Kokic et al., 2019). We compare our simulation results against both the linear ( ${\bar{A}}_{z} / \bar{A} = (z - 2) / 4$ where ${\bar{A}}_{z}$ is the average area of cell with $z$ neighbors and $\bar{A}$ is the average area of all cells) and quadratic relations ( ${\bar{A}}_{z} / \bar{A} = {(z / 6)}^{2}$ ) and find the agreement poor, as shown for the case of $ω_{cw} = 0.0025$ and $Ω = 0.4$ in Figure 3b (see also Appendix 1—figures 3 and 4 in Appendix 1). While in our simulations the cell monolayers are not always confluent due to the extrusion events, other studies with confluent cellular layers have also found such relations to not be valid (Kim et al., 2014; Wenzel et al., 2019). In our simulations, the projected area of an extruding cell decreases prior to extrusion, but the number of interacting neighbors generally does not change in that time frame (see Appendix 1—figure 4 in Appendix 1). Together, these results suggest mechanical rather than geometrical routes to cell extrusion. Specifically, in our approach cell extrusion emerges as a consequence of cells pushing and pulling on their neighbors due to their intrinsic activity. This contrasts with inherently threshold-based vertex models (see, e.g., Okuda and Fujimoto, 2020) for both cellular rearrangements (T1 transitions) and extrusions (T2 transitions).

Mechanical stress localization unifies distinct topological routes to cell extrusion

The correlation between disclinations and extrusion events is also related to the mechanical stress localization at the fivefold disclinations: The occurrence of disclinations in a flat surface produces local stress concentration (Irvine et al., 2010). Generally, it is energetically favorable to bend a flat surface, rather than to compress or to stretch it (Landau et al., 1986). Thus, the local stress concentration can lead to a fivefold (positive Gaussian curvature) or a sevenfold (negative Gaussian curvature) disclination (Seung and Nelson, 1988; Guitter and Kardar, 1990). In our set-up and given that we consider a rigid substrate, fivefold disclinations are much more likely to provide relief for the high local stress concentration. This can change if the rigidity of substrate is relaxed or extrusion in three-dimensional spheroids are considered. Since we conjecture that both topological defect- and disclination-mediated extrusion mechanisms are closely linked with stress localization, we characterize the in-plane and out-of-plane stresses associated with the simulated monolayer. We compute a coarse-grained stress field (Christoffersen et al., 1981; Li et al., 2022), $σ_{i j} = (1 / (2 V_{cg})) \sum_{m \in V_{cg}} ({\vec{T}}_{i} ({\vec{x}}_{m}) \otimes {\vec{e}}_{j}^{n} + {\vec{T}}_{j} ({\vec{x}}_{m}) \otimes {\vec{e}}_{i}^{n})$ , where ${\vec{x}}_{0}$ represents the center of the coarse-grained volume, $V_{cg} = ℓ_{stress}^{3}$ , corresponding to coarse-grained length $ℓ_{stress}$ and unit vector ${\vec{e}}_{i}^{n} = ({\vec{x}}_{0} - {\vec{x}}_{m}) / | {\vec{x}}_{0} - {\vec{x}}_{m} |$ . Herein, the stress fields are computed using $ℓ_{stress} = R_{0} / 4$ .

For the example simulation snapshot displayed in Figure 1a, at the onset of two extrusion events, we visualize normalized isotropic stress ${\tilde{σ}}^{iso} (\vec{x}) = σ^{iso} (\vec{x}) / σ_{max}^{iso}$ and out-of-plane shear ${\tilde{σ}}_{xz} (\vec{x}) = σ_{xz} (\vec{x}) / σ_{xz}^{max}$ , where $σ_{max}^{iso}$ and $σ_{x z}^{max}$ are the maximum values in their corresponding fields (see Figure 4a and b). We observe a high, out-of-plane, shear stress concentration (Figure 4b) as well as tensile and compressive stress pathways (Figure 4a) reminiscent of force chains in granular systems (Majmudar and Behringer, 2005).

Figure 4

Download asset Open asset

Temporal build-up of mechanical stress before extrusion events.

A representative analysis corresponding to the configuration shown in Figure 1a and projected into $x y -$ plane ( $z = 0$ , i.e., the basal side). (a) Normalized isotropic stress field ${\tilde{σ}}^{iso} (\vec{x}) = σ^{iso} (\vec{x}) / σ_{max}^{iso}$ , where $σ_{max}^{iso}$ is the maximum value of the isotropic stress field, and (b) normalized shear stress field, ${\tilde{σ}}_{xz} (\vec{x}) = σ_{xz} (\vec{x}) / σ_{xz}^{max}$ , where $σ_{xz}^{max}$ is the maximum value of $σ_{xz} (\vec{x})$ field. The symbol + denotes the center of mass for two extruding cells. (c) Cell (spatially) averaged normalized isotropic stress ${\bar{\tilde{σ}}}_{i}^{iso} (\tilde{t}) = {⟨ σ^{iso} (\vec{x}, \tilde{t}) ⟩}_{\vec{x} \in ℛ_{i}} / {⟨ σ^{iso} (\vec{x}, \tilde{t}) ⟩}_{\vec{x} \in ℛ}$ and (d) shear stress ${\bar{\tilde{σ}}}_{x z}^{i} (\tilde{t}) = {⟨ σ_{x z} (\vec{x}, \tilde{t}) ⟩}_{\vec{x} \in ℛ_{i}} / {⟨ σ_{x z} (\vec{x}, \tilde{t}) ⟩}_{\vec{x} \in ℛ}$ for an extruding cell $i$ during $\tilde{t} = (t / τ_{0}) \in [{\tilde{t}}_{e} - 2.5, {\tilde{t}}_{e} + 0.3125]$ , where ${\tilde{t}}_{e}$ denotes extrusion time and $τ_{0} = ξ R_{0} / α$ . The data shown in (**c, d**) correspond to all the considered parameters for cell–substrate ( $ω_{cw}$ ) and relative cell–cell adhesions ( $Ω$ ) and for four distinct realizations. Each gray line in the background represents an extruding cell, and the red line shows the mean and the standard deviation of the normalized stresses.

Figure 4c shows the evolution of spatially averaged normalized isotropic stress for extruding cell $i$ , ${\bar{\tilde{σ}}}_{i}^{iso} (\tilde{t}) = {⟨ σ^{iso} (\vec{x}, \tilde{t}) ⟩}_{\vec{x} \in ℛ_{i}} / {⟨ σ^{iso} (\vec{x}, \tilde{t}) ⟩}_{\vec{x} \in ℛ}$ , demonstrating a clear stress build-up, followed by a drop near $\tilde{t} - {\tilde{t}}_{e} = 0$ , as a cell detaches the substrate and loses contact with other cells, where ${\tilde{t}}_{e}$ is detected by our stress-independent criterion, $ℛ = ⋃_{i = 1}^{N} ℛ_{i}$ and $ℛ_{i}$ is the domain associated with cell i, $ℛ_{i} := {\vec{x} | ϕ_{i} (\vec{x}) \geq 0.5}$ .

Similarly, Figure 4d displays the spatially averaged normalized out-of-plane shear stress, ${\bar{\tilde{σ}}}_{x z}^{i} (\tilde{t}) = {⟨ σ_{x z} (\vec{x}, \tilde{t}) ⟩}_{\vec{x} \in ℛ_{i}} / {⟨ σ_{x z} (\vec{x}, \tilde{t}) ⟩}_{\vec{x} \in ℛ}$ , prior to a cell extrusion and for all extrusion events in our simulations, that is, nine cases and four realizations for each case. The shear stress prior to extrusion exhibits oscillations with large magnitudes relative to the mean field, a stark departure from their non-extruding counterparts (see Appendix 1—figure 4 in Appendix 1). This may indicate a hindrance to cell movement as we explore further next.

Interestingly, the association of cell extrusion events with regions of high out-of-plane shear stress has parallels with the phenomenon of plithotaxis, where it was shown that cells collectively migrate along the orientation of the minimal in-plane intercellular shear stress (Tambe et al., 2011). In this context, based on the association of cell extrusion events with regions of high out-of-plane shear stress, we conjecture that high shear stress concentration hinders collective cell migration with cell extrusion providing a mechanism to re-establish the status quo. This is also consistent with the observation we made earlier about large oscillations in ${\bar{\tilde{σ}}}_{x z}^{i}$ prior to an extrusion event for all extruding cells (Figure 4d).

Shifting tendencies towards extrusion at nematic and hexatic disclinations

The results so far clearly demonstrate the existence of mechanical routes for cell removal that are associated with nematic and hexatic disclinations and are governed by the in-plane and out-of-plane mechanical stress patterns in the cell assembly. The relative strength of cell–cell to cell–substrate adhesion, $Ω$ , further alters the likelihood of an extrusion event being associated with a +1/2 defect or a fivefold disclination. This is clearly observed in Figure 2d, which shows the first moment of the distribution for ${\tilde{d}}_{min}^{+ 1 / 2}$ , $m = ⟨ {\tilde{d}}_{min}^{+ 1 / 2} ⟩$ , increases with $Ω = ω_{cc} / ω_{cw}$ (see inset) while the peak of the probability density decreases with increasing $Ω$ . At the same time, the probability of an extrusion occurring at a fivefold disclination increases with increasing $Ω$ , as displayed in Figure 2e. However, nematic and hexatic order parameters do not show any clear trends with $Ω$ (see Appendix 1—figure 10 in Appendix 1). To better understand this tendency, we characterize the average isotropic stress fields around a +1/2 defect. This involves tracking each +1/2 defect starting from its nucleation and mapping the isotropic stress field, for each time step during the defect’s lifetime, in a square domain of size $L \approx 1.5 R_{0}$ centered on the defect location and accounting for its orientation, where $L$ is chosen based on the peak in ${\tilde{d}}_{min}^{+ 1 / 2}$ (see Figure 2d). An example for the normalized average isotropic stress field corresponding to $ω_{cw} = 0.0025$ and $Ω = 0.4$ is shown in Figure 5c, where ${\tilde{\bar{σ}}}^{iso} (\vec{x}) = {\bar{σ}}^{iso} (\vec{x}) / {\bar{σ}}_{max}^{iso}$ , with ${\bar{σ}}^{iso}$ representing the average field during defect lifetime, for all nucleated defects, and ${\bar{σ}}_{max}^{iso}$ is the maximum of the average field. This is in agreement with experimental measurements on epithelial monolayers (Saw et al., 2017; Balasubramaniam et al., 2021), with a compressive stress region near the head of the defect and a tensile region near the tail. Interestingly, there is an asymmetry in the intensity of stress in the compressive region at the head of the comet as opposed to the tensile region at the tail (≈5× higher). To expand on this observation, we focus on the probability density function for ${\tilde{\bar{σ}}}^{iso}$ and various $Ω$ . Remarkably, as shown in Figure 5a, with increasing $Ω$ , the peak of the probability density function decreases and the shoulders become wider, that is, the stress localization becomes more spread. At the same time, it is worth noting that the probability for the occurrence of a fivefold disclination increases as $Ω$ is increased, as shown in Figure 5b, while no clear trend is observed for the density of half-integer defects (see Appendix 1—figure 6 in Appendix 1). Therefore, the more spread localized stress is more likely to only clear the lower energetic barrier associated with buckling of a fivefold disclination (Seung and Nelson, 1988; Guitter and Kardar, 1990) – forming a positive Gaussian curvature – as opposed to cells with six neighbors. Furthermore, for a single disclination, this energy is higher for a sevenfold disclination (Deem and Nelson, 1996) and in our case the rigid substrate defies any attempts by a sevenfold disclination to buckle and form a negative Gaussian curvature. Together, these results provide a potential explanation for why as $Ω$ is increased, cells collectively have a tendency towards leveraging fivefold disclinations instead of +1/2 defects for extruding an unwanted cell.

Figure 5

Download asset Open asset

Spatial localization of mechanical stress leading to extrusion events.

(a) Probability density function for the normalized, ensemble average of isotropic stress field, ${\tilde{\bar{σ}}}^{iso}$ , projected into $x y$ -plane with $z = 0$ , that is, the basal side, around +1/2 defects for various cell–cell to cell–substrate adhesion ratios $Ω$ (colors correspond to legend in [a]). (b) Probability density function for the number of neighbors $z$ and various $Ω$ , for all cells and simulation time steps. (c) Normalized, ensemble average of isotropic stress field, ${\tilde{\bar{σ}}}^{iso} (\vec{x}) = {\bar{σ}}^{iso} (\vec{x}) / {\bar{σ}}_{max}^{iso}$ , with ${\bar{σ}}^{iso}$ representing the average field during defect life time, for nucleated defects, and ${\bar{σ}}_{max}^{iso}$ is the maximum of the average field, for the case $ω_{cw} = 0.0025$ and $Ω = 0.4$ .

Furthermore, one may naively think that only the distance between a half-integer nematic defect and an extrusion site is of importance. Such a view implicitly assumes the statistics we have presented (e.g., ${\tilde{d}}_{min}^{\pm 1 / 2}$ , ${\tilde{\bar{σ}}}^{iso}$ ) correspond to independent events, disregarding the highly heterogeneous nature of such a complex, active system. This heterogeneous nature manifests in stress fields, as shown in Figure 5a and c for the normalized ensemble average around a +1/2 defect. Therefore, the distance between a defect and an extrusion site, the intensity and the extent of the stress fields around that defect all play a role and are embedded in the statistics that we present in this work. In the future, it can be illuminating to study the effect of heterogeneity in the apical–basal mechanical response due to different mechanical properties and/or the nature of activity.

Conclusions

Our study presents a three-dimensional model of the collective migration-mediated cell elimination. Importantly, this framework allows for cell–substrate and cell–cell adhesion forces to be tuned independently. Our findings indeed suggest that varying the relative strength of cell–cell and cell–substrate adhesion can allow cells to switch between distinct mechanical pathways – leveraging defects in nematic and hexatic phases – to eliminate unwanted cells through: (i) cell extrusion at ±1/2 topological defects in the cell orientation field, consistent with experimental observations (Saw et al., 2017); and (ii) cell extrusion at fivefold disclinations in cell arrangement, where our results show a direct role of these disclinations in extruding the cells. Focusing on the extruded cells, the results demonstrate that increasing relative cell–cell adhesion increases the probability of an extruded cell being a fivefold disclination while weakening the correlation with +1/2 topological defects. This seems to emerge with a confluence of factors at play: (i) higher likelihood for a cell to be a fivefold disclination as $Ω = ω_{cc} / ω_{cw}$ is increased, (ii) more spread stress concentration around a +1/2 defect with increasing $Ω$ , and (iii) a higher likelihood for such a diffused local stress field to only reach the lower energy barrier associated with buckling a fivefold disclination (forming a positive Gaussian curvature) as opposed to cells with six neighbors as well as sevenfold disclinations. In the latter case, in addition to higher energy barrier, the rigid substrate denies a sevenfold disclination to create any negative Gaussian curvatures.

Additionally, the presented framework provides access to the local stress field, including the out-of-plane shear components. Access to this information led us to conjecture that high shear stress concentration frustrates collective cell migration with cell extrusion providing a pathway to re-establish the status quo. We expect these results to trigger further experimental studies of the mechanical routes to live cell elimination and probing the impact of tuning cell–cell and cell–substrate interactions, for example, by molecular perturbations of E-cadherin adhesion complexes between the cells and/or focal adhesion between cells and substrate, as performed recently in the context of topological defect motion in cell monolayers (Balasubramaniam et al., 2021). In this study, we intentionally narrowed our focus to the interplay of cell–cell and cell–substrate adhesion, without accounting for cell proliferation. In its absence, simulations with high extrusion events may lose confluency. However, the identified mechanical routes to extrusion prevail in cases with both high and low number of extrusions, where confluency is maintained.

Finally, we anticipate that this modeling framework opens the door to several interesting and unresolved problems in studying three-dimensional features of cell layers. In particular, the mechanics can be coupled with biochemistry to study a wider range of mechanisms that affect live cell elimination. Additionally, using our framework the substrate rigidity can be relaxed in the future studies to further disentangle the impacts of cell–substrate adhesion from substrate deformation due to cell generated forces. Similarly, three-dimensional geometries, such as spheroids or cysts, can be examined. The links between collective cell migration and granular physics, in terms of force chains and liquid-to-solid transition, as well as probing the impact of three-dimensionality and out-of-plane deformations on these processes, are exciting directions for future studies. Lastly, the coexistence of nematic and hexatic phases, their potential interactions, and their interplay with curved surfaces are promising avenues for extending the work presented here.

Materials and methods

We consider a cellular monolayer consisting of $N = 400$ cells on a substrate with its surface normal ${\vec{e}}_{n} (= {\vec{e}}_{z}) = {\vec{e}}_{x} \times {\vec{e}}_{y}$ and periodic boundaries in both ${\vec{e}}_{x}$ and ${\vec{e}}_{y}$ , where $({\vec{e}}_{x}, {\vec{e}}_{y}, {\vec{e}}_{z})$ constitute the global orthonormal basis (Figure 1). Cells are initiated on a two-dimensional simple cubic lattice and inside a cuboid of size $L_{x} = L_{y} = 320$ , $L_{z} = 64$ , grid size $a_{0} = 1$ and with radius $R_{0} = 8$ . The cell–cell and cell–substrate interactions have contributions from both adhesion and repulsion, in addition to self-propulsion forces associated with cell polarity. To this end, each cell $i$ is modeled as an active deformable droplet in three-dimensions using a phase-field, $ϕ_{i} = ϕ_{i} (\vec{x})$ . The interior and exterior of cell i corresponds to $ϕ_{i} = 1$ and $ϕ_{i} = 0$ , respectively, with a diffuse interface of length $λ$ connecting the two regions and the midpoint, $ϕ_{i} = 0.5$ , delineating the cell boundary. A three-dimensional extension of the 2D free energy functional (Palmieri et al., 2015; Aranson, 2016; Camley and Rappel, 2017; Mueller et al., 2019) is considered with additional contributions to account for cell–cell and cell–substrate adhesions:

\begin{array}{ll} F = \sum_{i}^{N} \frac{γ}{λ} \int d \vec{x} {4 ϕ_{i}^{2} {(1 - ϕ_{i})}^{2} + λ^{2} {(\vec{\nabla} ϕ_{i})}^{2}} + \\ \sum_{i}^{N} μ {(1 - \frac{1}{V_{0}} \int d \vec{x} ϕ_{i}^{2})}^{2} + \sum_{i}^{N} \sum_{j \neq i} \frac{κ_{c c}}{λ} \int d \vec{x} ϕ_{i}^{2} ϕ_{j}^{2} + \\ \sum_{i}^{N} \sum_{j \neq i} \frac{ω_{cc}}{λ^{2}} \int d \vec{x} \vec{\nabla} ϕ_{i} \cdot \vec{\nabla} ϕ_{j} + \sum_{i}^{N} \frac{κ_{c w}}{λ} \int d \vec{x} ϕ_{i}^{2} ϕ_{w}^{2} + \\ \sum_{i}^{N} \frac{ω_{cw}}{λ^{2}} \int d \vec{x} \vec{\nabla} ϕ_{i} \cdot \vec{\nabla} ϕ_{w}, \end{array}

where $ℱ$ contains a contribution due to the Cahn–Hilliard free energy (Cahn and Hilliard, 1958), which stabilizes the cell interface, followed by a soft constraint for cell volume around $V_{0} (= (4 / 3) π R_{0}^{3})$ , such that cells – each initiated with radius R₀ – are compressible. Additionally, $κ$ and $ω$ capture repulsion and adhesion between cell–cell (subscript $c c$ ) and cell–substrate (subscript $c w$ ), respectively. Moreover, $γ$ sets the cell stiffness and µ captures cell compressibility and $ϕ_{w}$ denotes a static phase-field representing the substrate. This approach resolves the cellular interfaces and provides access to intercellular forces. The dynamics for field $ϕ_{i}$ can be defined as:

\partial_{t} ϕ_{i} + {\vec{v}}_{i} \cdot \vec{\nabla} ϕ_{i} = - \frac{δ ℱ}{δ ϕ_{i}}, i = 1, \dots, N,

where $ℱ$ is defined in Equation (1), and ${\vec{v}}_{i}$ is the total velocity of cell $i$ . To resolve the forces generated at the cellular interfaces, we consider the following over-damped dynamics for cells:

{\vec{t}}_{i} = ξ {\vec{v}}_{i} - {\vec{F}}_{i}^{sp} = \int d \vec{x} ϕ_{i} \vec{\nabla} \cdot Π^{int},

where ${\vec{t}}_{i}$ denotes traction as defined for Bayesian Inversion Stress Microscopy in Saw et al., 2017, $ξ$ is substrate friction, and ${\vec{F}}_{i}^{sp} = α {\vec{p}}_{i}$ represents self-propulsion forces due to polarity, constantly pushing the system out-of-equilibrium. In this vein, $α$ characterizes the strength of polarity force and $Π^{int} = (\sum_{i} - (δ F / δ ϕ_{i})) 1$ . While only passive interactions are considered here, active nematic interactions can be readily incorporated in this framework (Mueller et al., 2019; Balasubramaniam et al., 2021). To complete the model, the dynamics of front-rear cell polarity is introduced based on contact inhibition of locomotion (CIL) (Abercrombie and Heaysman, 1954; Abercrombie, 1979) by aligning the polarity of the cell to the direction of the total interaction force acting on the cell (Smeets et al., 2016; Peyret et al., 2019). As such, the polarization dynamics is given by

\partial_{t} θ_{i} = - J | {\vec{t}}_{i} | Δ θ_{i} + D_{r} η,

where $θ_{i} \in [- π, π]$ is the angle associated with polarity vector, ${\vec{p}}_{i} = (\cos θ_{i}, \sin θ_{i}, 0)$ , and $η$ is the Gaussian white noise with zero mean, unit variance, $D_{r}$ is rotational diffusivity, $Δ θ_{i}$ is the angle between ${\vec{p}}_{i}$ and ${\vec{t}}_{i}$ , and positive constant $J$ sets the alignment time scale. It is worth noting that the self-propulsion forces, ${\vec{F}}_{i}^{sp}$ , associated with cell polarity, ${\vec{p}}_{i}$ , act in-plane but can induce out-of-plane components in force and velocity fields as a cell described by $ϕ_{i} (\vec{x})$ deforms in three-dimensions (see Equation 3).

We perform large-scale simulations with a focus on the interplay of cell–cell and cell–substrate adhesion strengths and its impact on cell expulsion from the monolayer. To this end, we set the cell–substrate adhesion strength $ω_{cw} \in {0.0015, 0.002, 0.0025}$ and vary the cell–substrate to cell–cell adhesion ratio in the range $Ω = ω_{cc} / ω_{cw} \in {0.2, 0.4, 0.6}$ . For each case in this study (total of nine), we simulate four distinct realizations with a total of $n_{sim} = 29, 000$ time steps. All results are reported in dimensionless units, introduced throughout the text, and the simulation parameters are chosen within the range that was previously shown to reproduce defect flow fields in epithelial layers (Balasubramaniam et al., 2021; see Appendix 1).

Appendix 1

Kolmogorov–Smirnov test for correlation between extrusion events and nematic defects

We use a hypothesis-testing approach to explore the existence of a correlation between the extrusion events in our simulations and nucleation of nematic topological defects. Specifically, we hypothesize that the extrusion events are uncorrelated with the topological defects. To falsify this, we use a Poisson point process to randomly generate extrusion events and quantify the minimum distance, ${\bar{d}}_{min} = d_{min} / R_{0}$ between the extrusion location and the nearest half-integer topological defects. To this end, for each simulation, we generate five realization of extrusion events using a Poisson point process with intensity set equal to the number of extrusions in that particular simulation. Furthermore, we assign an extrusion time, t_e, using a uniform distribution $t_{e} \sim U (1, n_{sim})$ . Then, we quantify the minimum distance ${\bar{d}}_{min}$ as we have done for our simulations. The results are shown in Appendix 1—figure 1. To falsify the stated hypothesis, we use a Kolmogorov–Smirnov (KS) test to measure if the two samples, one based in our simulations and one based on randomly generated extrusion events belong to the same distribution. As shown in Appendix 1—table 1, the results of the KS test falsify this hypothesis, that is, simulation based extrusion events are uncorrelated with the topological defects.

Appendix 1—figure 1

Download asset Open asset

Establishing the correlation between extrusion events and half-integer nematic defects via hypothesis testing.

Probability density of the normalized minimum distance, ${\bar{d}}_{min}^{+ 1 / 2}$ between an extrusion event and the nearest +1/2 topological defect (**a–c**) for $Ω = 0.2$ (a), $Ω = 0.4$ (b), and $Ω = 0.6$ (c). Probability density of the normalized minimum distance, ${\bar{d}}_{min}^{- 1 / 2}$ between an extrusion event and the nearest -1/2 topological defect (**d–f**) for $Ω = 0.2$ (d), $Ω = 0.4$ (e), and $Ω = 0.6$ (f).

Appendix 1—table 1

Results of the Kolmogorov–Smirnov (KS) tests for various cell–cell to cell–substrate adhesion ratio ( $Ω = ω_{cc} / ω_{cw}$ ) and for half-integer topological defects.

Both statistics and p-value are KS test results and $n$ corresponds to the number of samples in each distribution, for both simulations and the extrusion events generated through a Poisson point process.

Probability density	statistics	p-Value	$n$ (simulations)	$n$ (randomly generated)
$Ω = 0.2$ (+1/2)	0.8221	$1.95 \times 10^{- 154}$	426	402
$Ω = 0.4$ (+1/2)	0.228	$2.22 \times 10^{- 15}$	648	605
$Ω = 0.6$ (+1/2)	0.827	$5.12 \times 10^{- 212}$	551	570
$Ω = 0.2$ (–1/2)	0.802	$1.22 \times 10^{- 15}$	426	423
$Ω = 0.4$ (–1/2)	0.840	$8.17 \times 10^{- 260}$	648	660
$Ω = 0.6$ (–1/2)	0.824	$2.84 \times 10^{- 198}$	551	507

Simulation parameters

We perform large-scale simulations with a focus on the interplay of cell–cell and cell–substrate adhesion strengths on collective cell migration and its impact on cell expulsion from the monolayer. Following Mueller et al., 2019, the space and time discretization in our simulations are based on the average radius of MDCK cells, $\sim 5 μ m$ , velocity $\sim 20 μ m / h$ , and average pressure of ∼100 Pa, measured experimentally in MDCK monolayers (Saw et al., 2017), corresponding to $Δ x \sim 0.5 μ m$ , $Δ t \sim 0.1 s$ , and $Δ F \sim 1.5$ nN for force. In this study, we set the cell–substrate adhesion strength $ω_{cw} \in {0.0015, 0.002, 0.0025}$ and vary the cell–substrate to cell–cell adhesion ratio in the range $Ω = ω_{cc} / ω_{cw} \in {0.2, 0.4, 0.6}$ . Based on previous experimental and theoretical studies (Mueller et al., 2019; Peyret et al., 2019; Balasubramaniam et al., 2021), the other simulation parameters are $κ_{cc} = 0.5$ , $κ_{cw} = 0.15$ , $ξ = 1$ , $α = 0.05$ , $λ = 3$ , $μ = 45$ , $D_{r} = 0.01$ , and $J = 0.005$ , unless stated otherwise.

Lewis’s empirical relation

The empirical relationship proposed by Lewis, 1928 is generally valid for cellular matter that fill in the space without gaps. In our simulations, we can lose confluency due to cellular extrusions. Furthermore, even in the case of a confluent cellular layer, that is, no gaps between cells, Lewis’s law fails to capture the correlation between normalized area and number of neighbors (see Kim et al., 2014; Wenzel et al., 2019). However, we still investigated the existence of such correlation in our simulations. To this end, we used the following linear and quadratic relationships (Kokic et al., 2019):

\frac{{\bar{A}}_{z}}{\bar{A}} = \frac{z - 2}{4}

\frac{{\bar{A}}_{z}}{\bar{A}} = {(\frac{z}{6})}^{2}

where ${\bar{A}}_{z}$ is the average area of cells with $z$ neighbors and $\bar{A}$ is the average of area of the cells in the monolayer. As shown in Appendix 1—figure 2, the agreement is poor. Furthermore, while the projected area of an extruding cell decreases prior to extrusion, the number of neighbors its interacting with remains generally unchanged. This is shown in Appendix 1—figure 2.

Appendix 1—figure 2

Download asset Open asset

Lewis's relationship for extruding cells.

(**a–i**) Lewis’s relationship compared with our simulations for $ω_{cw} = 0.0015$ and $Ω = 0.2$ (a), $Ω = 0.4$ (b), and $Ω = 0.6$ (c); $ω_{cw} = 0.002$ and $Ω = 0.2$ (d), $Ω = 0.4$ (e), and $Ω = 0.6$ (f); and $ω_{cw} = 0.0025$ and $Ω = 0.2$ (g), $Ω = 0.4$ (h), and $Ω = 0.6$ (i).

Appendix 1—figure 3

Download asset Open asset

Coordination and area statistics of extruding cells.

(a) The average and standard deviation of coordination number – number of interacting neighbors – for an extruding cell. The data corresponds to all simulations and the four distinct realizations considered in the article. (b) Percent of extrusion events with a given $z$ – characterizing the change in number of interacting neighbors at extrusion time ${\tilde{t}}_{e}$ and ${\tilde{t}}_{e} - 2.5$ . The data corresponds to all simulations and the four distinct realizations considered in the article. (c) The temporal evolution of area for extruding cells normalized with the area at ${\tilde{t}}_{e} - 2.5$ , for one of the realizations.

Coordination number computation

To compute the coordination number, we use interaction between the cells instead of Voronoi tessellation. This is because when confluency is lost and there is a heterogeneous density of cells on the substrate, Voronoi tessellation would overestimate a cell’s number of neighbors. To this end, we consider two cells, $i$ and $j$ , as interacting cells if the following is satisfied:

{ϕ_{i} | ϕ_{i} > 0.25} \cap {ϕ_{j} | ϕ_{j} > 0.25} \neq \emptyset

Out-of-plane shear, $σ_{x z}$ , for non-extruding cells

The fluctuations in out-of-plane shear, ${\bar{\tilde{σ}}}_{x z}^{i} (\tilde{t}) = {⟨ σ_{x z} (\vec{x}, \tilde{t}) ⟩}_{\vec{x} \in ℛ_{i}} / {⟨ σ_{x z} (\vec{x}, \tilde{t}) ⟩}_{\vec{x} \in ℛ}$ for an extruding cell $i$ normalized by the maximum out-of-plane shear for all non-extruding cells in the same temporal window, as shown in Appendix 1—figure 4.

Appendix 1—figure 4

Download asset Open asset

Oscillations of out-of-plane shear stress for extruding cells relative to non-extruding cells.

(a) Temporal evolution of the out-of-plane shear stress for an extruding cell normalized by the maximum of the out-of-plane shear for all non-extruding cells within the same temporal window.

Half-integer defect statistics

We have computed the lifetime for a half-integer defect by tracking that defect from its nucleation to annihilation in our simulations. The probability density for defect lifetimes are shown in Appendix 1—figure 5. Furthermore, we computed the defect density, which we define as the number of defects, either +1/2 or -1/2, detected at each simulation (time) frame divided by the domain of the simulation, $a = L_{x} \times L_{y}$ . This is shown in Appendix 1—figure 6. We also computed the distances between half-integer defects for each simulation (time) frame, as shown in Appendix 1—figure 7. The peak in these probability densities are much larger than the peak of the minimum distance to ±1/2 defects, as shown in Figure 2c and d.

Appendix 1—figure 5

Download asset Open asset

Defect statistics.

(a) The probability density of +1∕2 and (b) -1/2 defect lifetimes.

Appendix 1—figure 6

Download asset Open asset

Probability density of (a) +1∕2 and (b) −1∕2 defect density distributions.

Appendix 1—figure 7

Download asset Open asset

Probability density of pairwise distance between (a) +1∕2 and +1∕2 defects, (b) −1∕2 and −1∕2 defects, and (c) +1∕2 and −1∕2 defects.

Phase diagram for extrusion intensity

To further explore the impact of asymmetric interaction of the cells with apical and basal sides, we have performed additional simulations varying the strength of the cell–substrate interactions. Appendix 1—figure 8 shows how changing cell–substrate adhesion (basal), $ω_{c w}$ , affects the extrusion rate. The results show that increasing cell–substrate adhesion leads to less extrusion events, while the ratio $Ω = ω_{c c} / ω_{c w}$ does not seem to play a significant role on the likelihood of an extrusion event occurring.

Appendix 1—figure 8

Download asset Open asset

Extrusion phase diagram.

(a) Mean and (b) standard deviation of the number of extrusions for a range of values for $Ω = ω_{c c} / ω_{c w}$ and cell–substrate adhesion, $ω_{c w}$ , over the range of 29,000 times steps and for four realizations.

Energy spectra

We calculated energy spectra for different cell-cell adhesion strengths, which suggests different power-law regimes, as shown in Appendix 1—figure 9. The kinetic energy spectrum, ${\hat{E}}_{v} = \frac{1}{2} ⟨ {\hat{v}}_{i} (k) {\hat{v}}_{i} (k) ⟩$ , where ${\hat{v}}_{i}$ is the Fourier transforms of the velocity field, and ${\tilde{E}}_{v} = {\hat{E}}_{v} (k) / {\hat{E}}_{v}^{max} (k)$ . Furthermore, $\tilde{k} = k / (2 π / R_{0})$ .

Appendix 1—figure 9

Download asset Open asset

Energy spectra for three sample simulations with (a) $Ω = 0.2$ , (b) $Ω = 0.4$ , and (c) $Ω = 0.6$ .

$p -$ atic order

We computed the $p -$ atic order parameter, associated with a liquid crystal exhibiting $p -$ fold rotational symmetry (Nelson and Halperin, 1979), $ψ_{p}^{i} = \frac{1}{N_{n}^{i}} \sum_{j}^{N_{n}^{i}} \exp (p i θ_{i j})$ , where $N_{n}^{i}$ is the number of neighbors for cell $i$ , $θ_{i j}$ is the angle between vector ${\vec{r}}_{i j}$ , connecting cell $i$ and neighboring cell $j$ , and ${\vec{e}}_{x}$ . Lastly, $p = 2$ for nematic and $p = 6$ for hexatic phases. The mean of the absolute value, $\bar{| ψ_{2} |}$ and $\bar{| ψ_{6} |}$ for various $Ω = ω_{c c} / ω_{c w}$ is shown in Appendix 1—figure 10.

Appendix 1—figure 10

Download asset Open asset

Temporal evolution of nematic order parameter (a) and hexatic order parameter (b) for various relative cell–cell adhesions, $Ω$ .

Example for equilibrated monolayer configuration

In the absence of activity, cells tend to equilibrate into a hexagonal lattice. An example is shown in Appendix 1—figure 11a along with the temporal evolution of the mean hexatic order parameter, $\bar{| ψ_{6} |}$ displayed in Appendix 1—figure 11b, plateauing at $\bar{| ψ_{6} |} = 1$ indicative of perfect hexatic order.

Appendix 1—figure 11

Download asset Open asset

In the absence of active forces, cells tend to equilibrate into a hexagonal lattice.

An example configuration (a) and evolution of the hexatic order parameter as the systems equilibrates (b).

Data availability

The current manuscript is a computational study, so no data have been generated for this manuscript. Modelling code is uploaded on the first author's GitHub page (https://github.com/siavashmonfared/siavashmonfared.github.io copy archived at Monfared, 2023).

References

1. Abercrombie M
2. Heaysman JEM
(1954) Observations on the social behaviour of cells in tissue culture
Experimental Cell Research 6:293–306.
https://doi.org/10.1016/0014-4827(54)90176-7
- PubMed
- Google Scholar
1. Abercrombie M
(1979) Contact inhibition and malignancy
Nature 281:259–262.
https://doi.org/10.1038/281259a0
- PubMed
- Google Scholar
(2022) Unifying polar and nematic active matter: Emergence and co-existence of half-integer and full-integer topological defects
Journal of Physics A 55:094002.
https://doi.org/10.1088/1751-8121/ac4abe
- Google Scholar
Book
1. Aranson IS
(2016) Physical Models of Cell Motility. Biological and Medical Physics, Biomedical Engineering
Cham: Springer International Publishing.
https://doi.org/10.1007/978-3-319-24448-8
- Google Scholar
Preprint
(2022) Epithelia are multiscale active liquid crystals
bioRxiv.
https://doi.org/10.1101/2022.02.01.478692
- Google Scholar
1. Balasubramaniam L
2. Doostmohammadi A
3. Saw TB
4. Narayana GHNS
5. Mueller R
6. Dang T
7. Thomas M
8. Gupta S
9. Sonam S
10. Yap AS
11. Toyama Y
12. Mège RM
13. Yeomans JM
14. Ladoux B
(2021) Author correction: Investigating the nature of active forces in tissues reveals how contractile cells can form extensile monolayers
Nature Materials 20:1167.
https://doi.org/10.1038/s41563-021-00974-9
- PubMed
- Google Scholar
(2020) Self-propelled rods: Insights and perspectives for active matter
Annual Review of Condensed Matter Physics 11:441–466.
https://doi.org/10.1146/annurev-conmatphys-031119-050611
- Google Scholar
1. Barriga EH
2. Franze K
3. Charras G
4. Mayor R
(2018) Tissue stiffening coordinates morphogenesis by triggering collective cell migration in vivo
Nature 554:523–527.
https://doi.org/10.1038/nature25742
- PubMed
- Google Scholar
(2018) Turbulent dynamics of epithelial cell cultures
Physical Review Letters 120:208101.
https://doi.org/10.1103/PhysRevLett.120.208101
- PubMed
- Google Scholar
1. Brugués A
2. Anon E
3. Conte V
4. Veldhuis JH
5. Gupta M
6. Colombelli J
7. Muñoz JJ
8. Brodland GW
9. Ladoux B
10. Trepat X
(2014) Forces driving epithelial wound healing
Nature Physics 10:683–690.
https://doi.org/10.1038/nphys3040
- PubMed
- Google Scholar
1. Cahn JW
2. Hilliard JE
(1958) Free energy of a nonuniform system
I. Interfacial Free Energy. Journal of Chemical Physics 28:258–267.
https://doi.org/10.1063/1.1744102
- Google Scholar
1. Camley BA
2. Rappel WJ
(2017) Physical models of collective cell motility: From cell to tissue
Journal of Physics D 50:11.
https://doi.org/10.1088/1361-6463/aa56fe
- PubMed
- Google Scholar
1. Chiou K
2. Collins EMS
(2018) Why we need mechanics to understand animal regeneration
Developmental Biology 433:155–165.
https://doi.org/10.1016/j.ydbio.2017.09.021
- PubMed
- Google Scholar
(1981) A micromechanical description of granular material behavior
Journal of Applied Mechanics 48:339–344.
https://doi.org/10.1115/1.3157619
- Google Scholar
(2005) Hexagonal packing of Drosophila wing epithelial cells by the planar cell polarity pathway
Developmental Cell 9:805–817.
https://doi.org/10.1016/j.devcel.2005.10.016
- PubMed
- Google Scholar
1. Deem MW
2. Nelson DR
(1996) Free energies of isolated five- and sevenfold disclinations in hexatic membranes
Physical Review. E, Statistical Physics, Plasmas, Fluids, and Related Interdisciplinary Topics 53:2551–2559.
https://doi.org/10.1103/physreve.53.2551
- PubMed
- Google Scholar
Book
1. de Genne PG
2. Prost J
(1998)
The Physics of Liquid Crystals. Number 83 in The International Series of Monographs on Physics

Clarendon Press.
- Google Scholar
1. De Pascalis C
2. Etienne-Manneville S
(2017) Single and collective cell migration: The mechanics of adhesions
Molecular Biology of the Cell 28:1833–1846.
https://doi.org/10.1091/mbc.E17-03-0134
- PubMed
- Google Scholar
1. Doostmohammadi A
2. Thampi SP
3. Saw TB
4. Lim CT
5. Ladoux B
6. Yeomans JM
(2015) Celebrating soft matter’s 10th anniversary: cell division: a source of active stress in cellular monolayers
Soft Matter 11:7328–7336.
https://doi.org/10.1039/c5sm01382h
- PubMed
- Google Scholar
(2016) Defect-Mediated morphologies in growing cell colonies
Physical Review Letters 117:048102.
https://doi.org/10.1103/PhysRevLett.117.048102
- PubMed
- Google Scholar
(2018) Spontaneous shear flow in confined cellular nematics
Nature Physics 14:728–732.
https://doi.org/10.1038/s41567-018-0099-7
- PubMed
- Google Scholar
1. Dwivedi VK
2. Pardo-Pastor C
3. Droste R
4. Kong JN
5. Tucker N
6. Denning DP
7. Rosenblatt J
8. Horvitz HR
(2021) Replication stress promotes cell elimination by extrusion
Nature 593:591–596.
https://doi.org/10.1038/s41586-021-03526-y
- PubMed
- Google Scholar
(2012) Crowding induces live cell extrusion to maintain homeostatic cell numbers in epithelia
Nature 484:546–549.
https://doi.org/10.1038/nature10999
- PubMed
- Google Scholar
(2018) Author correction: Scutoids are a geometrical solution to three-dimensional packing of epithelia
Nature Communications 9:4210.
https://doi.org/10.1038/s41467-018-06671-7
- PubMed
- Google Scholar
1. Guitter E
2. Kardar M
(1990) Tethering, crumpling, and melting transitions in hexatic membranes
Europhysics Letters 13:441–446.
https://doi.org/10.1209/0295-5075/13/5/011
- Google Scholar
(2022) Theory of defect-mediated morphogenesis
Science Advances 8:15.
https://doi.org/10.1126/sciadv.abk2712
- Google Scholar
1. Hogan C
2. Dupré-Crochet S
3. Norman M
4. Kajita M
5. Zimmermann C
6. Pelling AE
7. Piddini E
8. Baena-López LA
9. Vincent J-P
10. Itoh Y
11. Hosoya H
12. Pichaud F
13. Fujita Y
(2009) Characterization of the interface between normal and transformed epithelial cells
Nature Cell Biology 11:460–467.
https://doi.org/10.1038/ncb1853
- PubMed
- Google Scholar
(2010) Pleats in crystals on curved surfaces
Nature 468:947–951.
https://doi.org/10.1038/nature09620
- PubMed
- Google Scholar
1. Kaliman S
2. Hubert M
3. Wollnik C
4. Nuić L
5. Vurnek D
6. Gehrer S
7. Lovrić J
8. Dudziak D
9. Rehfeldt F
10. Smith AS
(2021) Mechanical regulation of epithelial tissue homeostasis
Physical Review X 11:031029.
https://doi.org/10.1103/PhysRevX.11.031029
- Google Scholar
(2017) Topological defects control collective dynamics in neural progenitor cell cultures
Nature 545:327–331.
https://doi.org/10.1038/nature22321
- PubMed
- Google Scholar
(2014) Lewis’ law revisited: The role of anisotropy in size-topology correlations
New Journal of Physics 16:015024.
https://doi.org/10.1088/1367-2630/16/1/015024
- PubMed
- Google Scholar
Preprint
(2019) Minimisation of Surface Energy Drives Apical Epithelial Organisation and Gives Rise to Lewis’ Law
bioRxiv.
https://doi.org/10.1101/590729
- Google Scholar
1. Ladoux B
(2009) Cells guided on their journey
Nature Physics 5:377–378.
https://doi.org/10.1038/nphys1281
- Google Scholar
1. Ladoux B
2. Mège RM
(2017) Mechanobiology of collective cell behaviours
Nature Reviews. Molecular Cell Biology 18:743–757.
https://doi.org/10.1038/nrm.2017.98
- PubMed
- Google Scholar
Book
(1986)
Theory of elasticity

In: Landau LD, Lifshitz EM, editors. Number Vol. 7 in Course of Theoretical Physics. Pergamon Press. pp. 7–44.
- Google Scholar
1. Le AP
2. Rupprecht JF
3. Mège RM
4. Toyama Y
5. Lim CT
6. Ladoux B
(2021) Adhesion-mediated heterogeneous actin organization governs apoptotic cell extrusion
Nature Communications 12:397.
https://doi.org/10.1038/s41467-020-20563-9
- PubMed
- Google Scholar
(2016) Tissue crowding induces caspase-dependent competition for space
Current Biology 26:670–677.
https://doi.org/10.1016/j.cub.2015.12.072
- PubMed
- Google Scholar
1. Lewis FT
(1928) The correlation between cell division and the shapes and sizes of prismatic cells in the epidermis of Cucumis
The Anatomical Record 38:341–376.
https://doi.org/10.1002/ar.1090380305
- Google Scholar
1. Li H
2. Matsunaga D
3. Matsui TS
4. Aosaki H
5. Kinoshita G
6. Inoue K
7. Doostmohammadi A
8. Deguchi S
(2022) Wrinkle force microscopy: A machine learning based approach to predict cell mechanics from images
Communications Biology 5:361.
https://doi.org/10.1038/s42003-022-03288-x
- PubMed
- Google Scholar
(2020) Solid-liquid transition of deformable and overlapping active particles
Physical Review Letters 125:038003.
https://doi.org/10.1103/PhysRevLett.125.038003
- PubMed
- Google Scholar
(2020) Chiral active hexatics: Giant number fluctuations, waves, and destruction of order
Physical Review Letters 125:238005.
https://doi.org/10.1103/PhysRevLett.125.238005
- PubMed
- Google Scholar
1. Majmudar TS
2. Behringer RP
(2005) Contact force measurements and stress-induced anisotropy in granular materials
Nature 435:1079–1082.
https://doi.org/10.1038/nature03805
- PubMed
- Google Scholar
1. Marinari E
2. Mehonic A
3. Curran S
4. Gale J
5. Duke T
6. Baum B
(2012) Live-Cell delamination counterbalances epithelial growth to limit tissue overcrowding
Nature 484:542–545.
https://doi.org/10.1038/nature10984
- PubMed
- Google Scholar
(2009) Quantifying cellular traction forces in three dimensions
PNAS 106:22108–22113.
https://doi.org/10.1073/pnas.0904565106
- PubMed
- Google Scholar
(2021) Bacteria solve the problem of crowding by moving slowly
Nature Physics 17:205–210.
https://doi.org/10.1038/s41567-020-01070-6
- Google Scholar
Software
1. Monfared S
(2023) Siavashmonfared.github.io, version swh:1:rev:544e446b755bb5949e3e2623fb18f38a98f22c7e
Software Heritage.

https://archive.softwareheritage.org/swh:1:dir:b0ac6611c1cebf2886d77d440f40c25deebaa83e;origin=https://github.com/siavashmonfared/siavashmonfared.github.io;visit=swh:1:snp:7efe820c2215b9bc9d5211ecc67e77f0f564dca2;anchor=swh:1:rev:544e446b755bb5949e3e2623fb18f38a98f22c7e
(2019) Emergence of active nematic behavior in monolayers of isotropic cells
Physical Review Letters 122:048004.
https://doi.org/10.1103/PhysRevLett.122.048004
- PubMed
- Google Scholar
1. Nelson DR
2. Halperin BI
(1979) Dislocation-mediated melting in two dimensions
Physical Review B 19:2457–2484.
https://doi.org/10.1103/PhysRevB.19.2457
- Google Scholar
1. Okuda S
2. Fujimoto K
(2020) A mechanical instability in planar epithelial monolayers leads to cell extrusion
Biophysical Journal 118:2549–2560.
https://doi.org/10.1016/j.bpj.2020.03.028
- PubMed
- Google Scholar
1. Palmieri B
2. Bresler Y
3. Wirtz D
4. Grant M
(2015) Multiple scale model for cell migration in monolayers: Elastic mismatch between cells enhances motility
Scientific Reports 5:11745.
https://doi.org/10.1038/srep11745
- PubMed
- Google Scholar
(2020) Hexatic phase in a model of active biological tissues
Soft Matter 16:3914–3920.
https://doi.org/10.1039/d0sm00109k
- PubMed
- Google Scholar
(2019) Sustained oscillations of epithelial cell sheets
Biophysical Journal 117:464–478.
https://doi.org/10.1016/j.bpj.2019.06.013
- PubMed
- Google Scholar
(2001) An epithelial cell destined for apoptosis signals its neighbors to extrude it by an actin- and myosin-dependent mechanism
Current Biology 11:1847–1857.
https://doi.org/10.1016/s0960-9822(01)00587-5
- PubMed
- Google Scholar
(2020) Collective cell mechanics of epithelial shells with organoid-like morphologies
Nature Communications 11:3805.
https://doi.org/10.1038/s41467-020-17535-4
- PubMed
- Google Scholar
(2021) Morphologies of compressed active epithelial monolayers
The European Physical Journal. E, Soft Matter 44:99.
https://doi.org/10.1140/epje/s10189-021-00094-x
- PubMed
- Google Scholar
1. Saw TB
2. Doostmohammadi A
3. Nier V
4. Kocgozlu L
5. Thampi S
6. Toyama Y
7. Marcq P
8. Lim CT
9. Yeomans JM
10. Ladoux B
(2017) Topological defects in epithelia govern cell death and extrusion
Nature 544:212–216.
https://doi.org/10.1038/nature21718
- PubMed
- Google Scholar
1. Seung HS
2. Nelson DR
(1988) Defects in flexible membranes with crystalline order
Physical Review. A, General Physics 38:1005–1018.
https://doi.org/10.1103/physreva.38.1005
- PubMed
- Google Scholar
1. Slattum GM
2. Rosenblatt J
(2014) Tumour cell invasion: An emerging role for basal epithelial cell extrusion
Nature Reviews. Cancer 14:495–501.
https://doi.org/10.1038/nrc3767
- PubMed
- Google Scholar
1. Smeets B
2. Alert R
3. Pešek J
4. Pagonabarraga I
5. Ramon H
6. Vincent R
(2016) Emergent structures and dynamics of cell colonies by contact inhibition of locomotion
PNAS 113:14621–14626.
https://doi.org/10.1073/pnas.1521151113
- PubMed
- Google Scholar
1. Sugimura K
2. Ishihara S
(2013) The mechanical anisotropy in a tissue promotes ordering in hexagonal cell packing
Development 140:4091–4101.
https://doi.org/10.1242/dev.094060
- PubMed
- Google Scholar
1. Tambe DT
2. Hardin CC
3. Angelini TE
4. Rajendran K
5. Park CY
6. Serra-Picamal X
7. Zhou EH
8. Zaman MH
9. Butler JP
10. Weitz DA
11. Fredberg JJ
12. Trepat X
(2011) Collective cell guidance by cooperative intercellular forces
Nature Materials 10:469–475.
https://doi.org/10.1038/nmat3025
- PubMed
- Google Scholar
1. Tan TH
2. Liu J
3. Miller PW
4. Tekant M
5. Dunkel J
6. Fakhri N
(2020) Topological turbulence in the membrane of a living cell
Nature Physics 16:657–662.
https://doi.org/10.1038/s41567-020-0841-9
- Google Scholar
(2008) Apoptotic force and tissue dynamics during Drosophila embryogenesis
Science 321:1683–1686.
https://doi.org/10.1126/science.1157052
- PubMed
- Google Scholar
1. Vafa F
2. Mahadevan L
(2022) Active nematic defects and epithelial morphogenesis
Physical Review Letters 129:098102.
https://doi.org/10.1103/PhysRevLett.129.098102
- PubMed
- Google Scholar
(2019) Topological and geometrical quantities in active cellular structures
The Journal of Chemical Physics 150:164108.
https://doi.org/10.1063/1.5085766
- PubMed
- Google Scholar
1. Wenzel D
2. Voigt A
(2021) Multiphase field models for collective cell migration
Physical Review. E 104:054410.
https://doi.org/10.1103/PhysRevE.104.054410
- PubMed
- Google Scholar
1. Zhang J
2. Yang N
3. Kreeger PK
4. Notbohm J
(2021) Topological defects in the mesothelium suppress ovarian cancer cell clearance
APL Bioengineering 5:036103.
https://doi.org/10.1063/5.0047523
- PubMed
- Google Scholar

Article and author information

Author details

Siavash Monfared

Division of Engineering and Applied Science, California Institute of Technology, Pasadena, United States

Contribution
Conceptualization, Software, Formal analysis, Investigation, Visualization, Methodology, Writing – original draft, Writing – review and editing

For correspondence
monfared@caltech.edu

Competing interests
No competing interests declared

"This ORCID iD identifies the author of this article:" 0000-0002-7629-7977
Guruswami Ravichandran

Division of Engineering and Applied Science, California Institute of Technology, Pasadena, United States

Contribution
Conceptualization, Resources, Supervision, Funding acquisition, Investigation, Writing – review and editing

Competing interests
No competing interests declared
José Andrade

Division of Engineering and Applied Science, California Institute of Technology, Pasadena, United States

Contribution
Conceptualization, Resources, Supervision, Funding acquisition, Investigation, Writing – review and editing

Competing interests
No competing interests declared
Amin Doostmohammadi

Niels Bohr Institute, University of Copenhagen, Copenhagen, Denmark

Contribution
Conceptualization, Resources, Supervision, Funding acquisition, Investigation, Methodology, Writing – review and editing

For correspondence
doostmohammadi@nbi.ku.dk

Competing interests
No competing interests declared

Funding

Multidisciplinary University Research Initiative (W911NF-19-1-0245)

Guruswami Ravichandran
José Andrade

Novo Nordisk Fonden (NNF18SA0035142)

Amin Doostmohammadi

Villum Fonden (29476)

Amin Doostmohammadi

Novo Nordisk Foundation (NNF21OC0068687)

Amin Doostmohammadi

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

Acknowledgements

SM is grateful for the generous support of the Rosenfeld Foundation fellowship at the Niels Bohr Institute, University of Copenhagen. SM, GR, and JA acknowledge support for this research provided by US ARO funding through the Multidisciplinary University Research Initiative (MURI) grant no. W911NF-19-1-0245. AD acknowledges funding from the Novo Nordisk Foundation (grant no. NNF18SA0035142 and NERD grant No. NNF21OC0068687), Villum Fonden grant no. 29476, and the European Union via the ERC-Starting Grant PhysCoMeT. The authors would like to thank Dr. Lakshmi Balasubramaniam and Prof. Benoît Ladoux (Institut Jacques Monod, University Paris City), Guanming Zhang and Prof. Julia M Yeomans (The Rudolf Peierls Centre for Theoretical Physics, University of Oxford), Prof. Jörn Dunkel (Mathematics Department, MIT), and Prof. M Cristina Marchetti (Department of Physics, University of California Santa Barbara) for helpful discussions. The authors are also grateful for the comments and the feedback provided by the anonymous reviewers.

Version history

Preprint posted: August 18, 2021 (view preprint)
Received: August 3, 2022
Accepted: March 22, 2023
Version of Record published: April 18, 2023 (version 1)

Copyright

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.

Metrics

2,061

views
264

downloads
6

citations

Views, downloads and citations are aggregated across all versions of this paper published by eLife.

Download links

A two-part list of links to download the article, or parts of the article, in various formats.

Downloads (link to download the article as PDF)

Article PDF

Open citations (links to open the citations from this article in various online reference manager services)

Mendeley

Cite this article (links to download the citations from this article in formats compatible with various reference manager tools)

Siavash Monfared
Guruswami Ravichandran
José Andrade
Amin Doostmohammadi

(2023)

Mechanical basis and topological routes to cell elimination

eLife 12:e82435.

https://doi.org/10.7554/eLife.82435

Share this article

Cite this article

Cell extrusion in a 3D representation of a confluent cell layer.

Nematic and hexatic disclinations govern cell extrusion.

Topological, rather than geometrical, route to cell extrusion.

Temporal build-up of mechanical stress before extrusion events.

Spatial localization of mechanical stress leading to extrusion events.

Establishing the correlation between extrusion events and half-integer nematic defects via hypothesis testing.

Results of the Kolmogorov–Smirnov (KS) tests for various cell–cell to cell–substrate adhesion ratio (Ω=ωcc/ωcw) and for half-integer topological defects.

Lewis's relationship for extruding cells.

Coordination and area statistics of extruding cells.

Oscillations of out-of-plane shear stress for extruding cells relative to non-extruding cells.

Defect statistics.

Probability density of (a) +1∕2 and (b) −1∕2 defect density distributions.

Probability density of pairwise distance between (a) +1∕2 and +1∕2 defects, (b) −1∕2 and −1∕2 defects, and (c) +1∕2 and −1∕2 defects.

Extrusion phase diagram.

Energy spectra for three sample simulations with (a) Ω=0.2, (b) Ω=0.4, and (c) Ω=0.6.

Temporal evolution of nematic order parameter (a) and hexatic order parameter (b) for various relative cell–cell adhesions, Ω.

In the absence of active forces, cells tend to equilibrate into a hexagonal lattice.

Author details

Siavash Monfared

Contribution

For correspondence

Competing interests

Guruswami Ravichandran

Contribution

Competing interests

José Andrade

Contribution

Competing interests

Amin Doostmohammadi

Contribution

For correspondence

Competing interests

Downloads (link to download the article as PDF)

Open citations (links to open the citations from this article in various online reference manager services)

Cite this article (links to download the citations from this article in formats compatible with various reference manager tools)

Categories and tags

Further reading

Results of the Kolmogorov–Smirnov (KS) tests for various cell–cell to cell–substrate adhesion ratio ( $Ω = ω_{cc} / ω_{cw}$ ) and for half-integer topological defects.

Energy spectra for three sample simulations with (a) $Ω = 0.2$ , (b) $Ω = 0.4$ , and (c) $Ω = 0.6$ .

Temporal evolution of nematic order parameter (a) and hexatic order parameter (b) for various relative cell–cell adhesions, $Ω$ .