Mechanical basis and topological routes to cell elimination

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 remain 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 five-fold, 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.


INTRODUCTION
The ability of cells to self-organize and to collectively migrate drives numerous physiological processes including morphogenesis [1,2], epithelial-mesenchymal transition [3], wound healing [4], and tumor progression [5].Advanced experimental techniques have linked this ability to mechanical interactions between cells [6][7][8].Specifically, cells actively coordinate their movements through mechanosensitive adhesion complexes at the cellsubstrate interface and cell-cell junctions.Moreover, cellcell and cell-substrate adhesions seem to be coupled [9], 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 integrityare 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 [10].This process can get initiated through apoptotic signaling [10], oncogenic transformation [11], and overcrowding of cells [12][13][14] or induced by replication stress [15].Most importantly, cell extrusion plays an important role in developmental [16], homeostatic [13,17] and pathological processes [18], 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 [19][20][21][22][23][24] and hexatic order [ [25][26][27][28][29] with the two phases potentially coexisting [30] and interacting provide a fresh perspective for understanding cellular processes.The five-fold disclinations in hexatic arrangement of cells are numerically shown to favor overlaps between the cells in two-dimensions [31], potentially contributing to the cell extrusion in threedimensions.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 [29,32,33].Moreover, in cellular monolayers, cometand trefoil-shaped half-integer topological defects, cor-responding to +1/2 and −1/2 charges, respectively, are prevalent [34,35].These are singular points in cellular alignment that mark the breakdown of orientational order [36].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 [19].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 [10].Nevertheless, it is not clear if 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 five-fold disclinations, and (ii) the interplay of cell-cell and cellsubstrate 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.

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 Fig. 16 in SI for an example).As we introduce self-propulsion forces associated with front-rear cell polarity (see Methods for polarisation dynamics), the system is pushed away from its equilibrium hexagonal configuration, resulting in defects manifested as five-fold and seven-fold disclinations, as shown in Fig. 2(b).Fig. 1(a) 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 0 is the initial cell radius.Fig. 1(b) displays the out-of-plane normalized velocity profile, ṽz = ( v ( x) • e z ) /v max z where v max z is the maximum value of the velocity component in 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.In order to probe the possible mechanical routes to cell extrusion, we begin with 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 xy−plane (z = 0, i.e. the basal side) and identify topological defects as the singularities in the orientation field.The results (example snapshot in Fig. 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 Fig. 14 in SI for energy spectra characterization).It is noteworthy that unlike previous studies of active nematic behavior in 2D cell layers [37,38], 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 inline with experimental observations in cell monolayers [19,22], discrete models of self-propelled rods [39,40], and recently proposed continuum model of polar active matter [41].
Remarkably, in accordance with experimental observations [19], 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, Figs.2(c)-(d min /R 0 between an extruding cell and ±1/2 topological defects in the interval t ∈ [ te − 5.625, te + 0.625], where t = t/τ 0 is the normalized time, τ 0 = ξR 0 /α, ξ corresponds to cell-substrate friction, α denotes the strength of polarity force and te is the (normalized) extrusion time.This temporal window is chosen based on the first moment of a defect's lifetime distribution (see Fig. 10 in SI).The data in Figs.2(c)-(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 (≈ 1.5R 0 ), at a much smaller distance relative to a typical distance between two defects (see Fig. 12 in SI), and falls off to nearly zero for d ±1/2 min 5R 0 (= 40).Furthermore, laser ablation experiments have established that an induced extrusion event does not favor the nucleation of a pair of nematic defects [19].
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 ∼ U (1, n sim ), where n sim = 29, 000 is total number of time steps.As an example, Fig. 3(a) shows probability density function for d+1/2 min 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 rejects this (see Tab.I and Fig. 6 in SI) and thus falsify the hypothesis that simulation based extrusion events are uncorrelated with the half-integer nematic defects.
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 five-fold and seven-fold disclinations (see Fig. 2(b)).To quantify the relation between extrusion events and the disclinations, the probability density of the coordination number of an extruding cell ( dmin = 0) averaged over the time interval, t ∈ [ te −5.625, te +0.625], z, for all the realizations is shown in Fig. 2(e), clearly exhibiting a sharp peak near z = 5.The coordination number is determined based on the interactions of cells (see SI) and this property is independent of apical or basal considerations [42], unlike geometrical structures called scutoids that have been identified in curved epithelial tubes [43].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 Fig. 13 in SI).
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, i.e. con-fluent state -is a linear one and it is due to F. Lewis [44] with other types of relations, e.g.quadratic, proposed since his work [45].We compare our simulation results against both the linear ( Āz / Ā = (z − 2) /4 where Āz is the average area of cell with z neighbors and Ā is the average area of all cells) and quadratic relations ( Āz / Ā = (z/6) 2 ) and find the agreement poor, as shown for the case of ω cw = 0.0025 and Ω = 0.4 in Fig. 3(b) (see also Figs. 7 and 8 in SI).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 the such relations to not be valid [46,47].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 Fig. 8 in SI).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.[48]), 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 five-fold disclinations: The occurrence of disclinations in a flat surface produces local stress concentration [49].Generally, it is energetically favorable to bend a flat surface, rather than to compress or to stretch it [50].Thus, the local stress concentration can lead to a five-fold (positive Gaussian curvature) or a seven-fold (negative Gaussian curvature) disclination [51,52].In our set-up and given that we consider a rigid substrate, five-fold 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 disclinationmediated extrusion mechanisms are closely linked with stress localization, we characterize the in-plane and outof-plane stresses associated with the simulated monolayer.We compute a coarse-grained stress field [53,54], where x 0 represents the center of the coarse-grained volume, V cg = 3 stress , corresponding to coarse-grained length stress and unit vector e n i = ( Herein, the stress fields are computed using stress = R 0 /4. For the example simulation snapshot displayed in Fig. 1(a), at the onset of two extrusion events, we visualize  of-plane, shear stress concentration (Fig. 4(b)) as well as tensile and compressive stress pathways (Fig. 4(a)) reminiscent of force chains in granular systems [55].
Figure 4(c) shows the evolution of spatially averaged normalized isotropic stress for extruding cell i, σiso i t = σ iso x, t x∈Ri / σ iso x, t x∈R , demonstrating a clear stress build up, followed by a drop near t− te = 0, as a cell detaches the substrate and loses contact with other cells, where te is detected by our stress-independent criterion, R = N i=1 R i and R i is the domain associated with cell i, R i := { x|φ i ( x) ≥ 0.5}.Similarly, Fig. 4(d) displays the spatially averaged normalized out-of-plane shear stress, σi xz t = σ xz x, t x∈Ri / σ xz x, t x∈R prior to a cell extrusion and for all extrusion events in our simulations, i.e. 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 Fig. 9 in SI).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 [56].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 σi xz prior to an extrusion event, for all extruding cells (Fig. 4(d)).

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 cellcell to cell-substrate adhesion, Ω, further alter the likelihood of an extrusion event being associated with a +1/2 defect or a five-fold disclination.This is clearly observed in Fig. 2(d), that show the first moment of the distribution for d+1/2 min , m = d+1/2 min , 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 five-fold disclination increases with increasing Ω, as displayed in Fig. 2(e).However, nematic and hexatic order parameters do not show any clear trends with Ω (see Fig. 15 in SI).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 life time, in a square domain of size L ≈ 1.5R 0 centered on the defect location and accounting for its orientation, where L is chosen based on the peak in d+1/2 min (see Fig. 2(d)).An example for the normalized average isotropic stress field corresponding to ω cw = 0.0025 and Ω = 0.4 is shown in Fig. 5(c), where σiso ( x) = σiso ( x) /σ iso max , with σiso representing the average field during defect life time, for all nucleated defects, and σiso max is the maximum of the average field.This is in agreement with experimental measurements on epithelial monolayers [9,19], 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 σiso and various Ω.Remarkably, as shown in Fig. 5(a), with increasing Ω, the peak of the probability density function decreases and the shoulders become wider, i.e. the stress localization becomes more spread.At the same time, it is worth noting that the probability for the occurrence of a five-fold disclination increases as Ω is increased, as shown in Fig. 5(b), while no clear trend is observed for the density of halfinteger defects (see Fig. 11 in SI).Therefore, the more spread localized stress is more likely to only clear the lower energetic barrier associated with buckling of a fivefold disclination [51,52] -forming a positive Gaussian curvature -as opposed to cells with six neighbors.Furthermore for a single disclination, this energy is higher for a seven-fold disclination [57] and in our case the rigid substrate defies any attempts by a seven-fold 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 five-fold disclinations instead of +1/2 defects for extruding an unwanted cell.
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.d±1/2 min , σ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 Figs.5(a) and 5(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 cellcell 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 -leverag-ing 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 [19]; and (ii) cell extrusion at five-fold 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 five-fold 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 five-fold 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 five-fold disclination (forming a positive Gaussian curvature) as opposed to cells with six neighbors as well as seven-fold disclinations.In the latter case, in addition to higher energy barrier, the rigid substrate denies a seven-fold 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 [9].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 threedimensional 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 co-existence 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 e n (= e z ) = e x × e y and periodic boundaries in both e x and e y , where ( e x , e y , e z ) constitute the global orthonormal basis (Fig. 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 selfpropulsion 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 ( 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 [37,[58][59][60] is considered with additional contributions to account for cell-cell and cell-substrate adhesions: where F contains a contribution due to the Cahn-Hilliard free energy [61] which stabilizes the cell interface, followed by a soft constraint for cell volume around V 0 = (4/3) πR 3 0 , such that cells -each initiated with radius R 0 -are compressible.Additionally, κ and ω capture repulsion and adhesion between cell-cell (subscript cc) and cell-substrate (subscript cw), 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: where F is defined in Eq. ( 1) and 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: where t i denotes traction as defined for Bayesian Inversion Stress Microscopy in [19], ξ is substrate friction and F sp i = α 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 = ( i − (δF/δφ i )) 1.While only passive interactions are considered here, active nematic interactions can be readily incorporated in this framework [9,37].To complete the model, the dynamics of front-rear cell polarity is introduced based on contact inhibition of locomotion (CIL) [62,63] by aligning the polarity of the cell to the direction of the total interaction force acting on the cell [64,65].As such, the polarization dynamics is given by: where θ i ∈ [−π, π] is the angle associated with polarity vector, 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 p i and t i , and positive constant J sets the alignment time scale.It is worth noting that the self-propulsion forces, F sp i , associated with cell polarity, p i , acts in-plane but can induce out-of-plane components in force and velocity fields as a cell described by φ i ( x) deforms in three-dimensions (see Eq.( 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 ∈ {0.0015, 0.002, 0.0025} and vary the cell-substrate to cell-cell adhesion ratio in the range Ω = ω cc /ω cw ∈ {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 [9] (see SI). .

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, dmin = 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 ∼ U (1, n sim ).Then, we quantify the minimum distance dmin as we have done for our simulations.The results are shown in Fig. 6.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 Table I, the results of the KS test falsifies this hypothesis, i.e. simulation based extrusion events are uncorrelated with the topological defects.
Table I. 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.

Figure 1 .
Figure 1.Cell extrusion in a 3D representation of a confluent cell layer.(a): A representative simulation snapshot (cell-substrate adhesion ωcw = 0.0025 and relative cellcell adhesion Ω = ωcc/ωcw = 0.4) of a three-dimensional cell monolayer.Two cells are visibly extruding.(b): A crosssection (dotted yellow line A − A ) of the cell monolayer highlighting the two extruding cells via the normalized outof-plane velocity (ṽz = ( v • ez) /v max z ), where v max z is the maximum value of the vz component of the velocity field v in the shown cross-section.
) display the probability density of the normalized minimum distance d±1/2 min = d ±1/2

Figure 2 .
Figure 2. Nematic and hexatic disclinations govern cell extrusion.A representative analysis corresponding to the configuration shown in Fig. 1(a) and projected into xy− plane (z = 0, i.e. the basal side).(a): a coarse-grained director field with coarse-graining length of one cell size dir.= R0 and +1/2 (filled circles with the line indicating orientation) and −1/2 (three connected lines with 3-fold symmetry) nematic defects.(b): Number of neighbors z for each cell, including fivefold and seven-fold 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, dmin = dmin/R0, for varying cell-cell to cell-substrate adhesion ratios Ω for (c) −1/2 and (d) +1/2 topological defects (inset: distribution mean m = d+1/2 min vs.Ω).(e): The probability density of average coordination number z for an extruding cell during t = (t/τ0) ∈ [ te − 2.5, te + 0.3125], where te denotes extrusion time, τ0 = ξR0/α and for varying cell-substrate to cell-cell adhesion ratios Ω.The data in (c)-(e) corresponds to four different realizations.

Figure 3 .
Figure 3. 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, d+1/2 min , 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.Āz is the average area for cells with z neighbors and Ā is the average area of all cells.The colorbar indicates simulation time step and the data correspond to the case ωcw = 0.0025 and Ω = 0.4.

Figure 4 .
Figure 4. Temporal build-up of mechanical stress before extrusion events.A representative analysis corresponding to the configuration shown in Fig. 1(a) and projected into xy− plane (z = 0, i.e. the basal side).(a): Normalized isotropic stress field σiso ( x) = σ iso ( x) /σ iso max where σ iso max is the maximum value of the isotropic stress field and (b) normalized shear stress field, σxz ( x) = σxz ( x) /σ max xz , where σ max xz is the maximum value of σxz ( x) field.The symbol + denotes the center of mass for two extruding cells.(c): cell (spatially) averaged normalized isotropic stress σiso i t = σ iso x, t x∈R i / σ iso x, t x∈R and (d) shear stress σi xz t = σxz x, t x∈R i / σxz x, t x∈R for an extruding cell i during t = (t/τ0) ∈ [ te − 2.5, te + 0.3125], where te denotes extrusion time and τ0 = ξR0/α.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, the red line shows the mean and the standard deviation of the normalized stresses.

Figure 5 .
Figure 5. Spatial localization of mechanical stress leading to extrusion events.(a): Probability density function for the normalized, ensemble average of isotropic stress field, σiso , projected into xy-plane with z = 0 i.e. 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, σiso ( x) = σiso ( x) /σ iso max , with σiso representing the average field during defect life time, for nucleated defects, and σiso max is the maximum of the average field, for the case ωcw = 0.0025 and Ω = 0.4.

Figure 8 .
Figure 8.(a) The average 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 manuscript.(b) Percent of extrusion events with a given z -characterizing the change in number of interacting neighbors at extrusion time te and te − 2.5.The data corresponds to all simulations and the four distinct realizations considered in manuscript.(c) The temporal evolution of area for extruding cells normalized with the area at te − 2.5, for one of the realizations.

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

Figure 13 .
Figure 13.(a) mean and (b) standard deviation of the number of extrusions for a range of values for Ω = ωcc/ωcw and cell-substrate adhesion, ωcw over the range of 29,000 times steps and for four realizations.

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

Figure 16 .
Figure 16.In 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).