Mechanical basis and topological routes to cell elimination
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 phasefield model of a threedimensional 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 halfinteger, nematic, defects for extruding a cell. We unify our findings by accessing threedimensional 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 threedimensional phasefield model for cell monolayers and use this to investigate the relationship between singlecell 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 cellcell and cellsubstrate 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.sa0Introduction
The ability of cells to selforganize 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 EtienneManneville, 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 twodimensional systems, hindering indepth exploration of intrinsically threedimensional 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 threedimensional. Thus, studying the underlying mechanisms necessitates access to both inplane and outofplane 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; BlanchMercader 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 (ArmengolCollado 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 twodimensions (Loewe et al., 2020), 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 domelike 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 trefoilshaped halfinteger 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 twodimensional models of the cell layers, fundamental questions about the threedimensional phenomenon of cell extrusion and its connection to the interplay between cellgenerated forces at the interface between cells and the substrate, with multicellular force transmission across the cell layer, remain unanswered.
In this article, we explore threedimensional collective cell migration in cellular monolayers. Based on largescale 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 threedimensional 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 selfpropulsion 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 selfpropulsion forces associated with frontrear 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_{0} is the initial cell radius. Figure 1b displays the outofplane normalized velocity profile, ${\overrightarrow{\stackrel{~}{v}}}_{z}=\left(\overrightarrow{v}\left(\overrightarrow{x}\right)\cdot {\overrightarrow{e}}_{z}\right)/{v}_{z}^{\text{max}}$ where ${v}_{z}^{\text{max}}$ is the maximum value of the velocity component in ${\overrightarrow{e}}_{z}$ direction in the displayed crosssection 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 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 $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 Figure 2a) show the continuous emergence of halfinteger (±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 selfpropulsion, 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; BlanchMercader et al., 2018), discrete models of selfpropelled rods (Bär et al., 2020; Meacock et al., 2021), and recently proposed continuum model of polar active matter (Amiri et al., 2022).
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 cometshaped and 1/2 trefoilshaped topological defects. To quantify this, Figure 2c and d display the probability density of the normalized minimum distance ${\stackrel{~}{d}}_{\text{min}}^{\pm 1/2}={d}_{\text{min}}^{\pm 1/2}/{R}_{0}$ between an extruding cell and ±1/2 topological defects in the interval $\stackrel{~}{t}\in [{\stackrel{~}{t}}_{e}5.625,{\stackrel{~}{t}}_{e}+0.625]$, where $\stackrel{~}{t}=t/{\tau}_{0}$ is the normalized time, ${\tau}_{0}=\xi {R}_{0}/\alpha $, $\xi $ corresponds to cell–substrate friction, $\alpha $ denotes the strength of polarity force, and ${\stackrel{~}{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, $\mathrm{\Omega}={\omega}_{\text{cc}}/{\omega}_{\text{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}_{\text{min}}^{\pm 1/2}5{R}_{0}\phantom{\rule{veryverythickmathspace}{0ex}}\left(=40\right)$. 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 hypothesistesting 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 halfinteger 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}_{\text{sim}})$, where ${n}_{\text{sim}}=29,000$ is total number of time steps. As an example, Figure 3a shows probability density function for ${\stackrel{~}{d}}_{\text{min}}^{+1/2}$ for $\mathrm{\Omega}=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 simulationbased extrusion events are uncorrelated with the halfinteger 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 phasefield 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 (${\stackrel{~}{d}}_{\text{min}}=0$) averaged over the time interval, $\stackrel{~}{t}\in [{\stackrel{~}{t}}_{e}5.625,{\stackrel{~}{t}}_{e}+0.625]$, $\overline{z}$, for all the realizations is shown in Figure 2e, clearly exhibiting a sharp peak near $\overline{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ómezGá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 (${\overline{A}}_{z}/\overline{A}=\left(z2\right)/4$ where ${\overline{A}}_{z}$ is the average area of cell with $z$ neighbors and $\overline{A}$ is the average area of all cells) and quadratic relations (${\overline{A}}_{z}/\overline{A}={\left(z/6\right)}^{2}$) and find the agreement poor, as shown for the case of ${\omega}_{\text{cw}}=0.0025$ and $\mathrm{\Omega}=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 thresholdbased 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 setup 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 threedimensional spheroids are considered. Since we conjecture that both topological defect and disclinationmediated extrusion mechanisms are closely linked with stress localization, we characterize the inplane and outofplane stresses associated with the simulated monolayer. We compute a coarsegrained stress field (Christoffersen et al., 1981; Li et al., 2022), ${\sigma}_{ij}=\left(1/(2{V}_{\text{cg}})\right){\sum}_{m\in {V}_{\text{cg}}}\left({\overrightarrow{T}}_{i}\left({\overrightarrow{x}}_{m}\right)\otimes {\overrightarrow{e}}_{j}^{n}+{\overrightarrow{T}}_{j}\left({\overrightarrow{x}}_{m}\right)\otimes {\overrightarrow{e}}_{i}^{n}\right)$, where ${\overrightarrow{x}}_{0}$ represents the center of the coarsegrained volume, ${V}_{\text{cg}}={\mathrm{\ell}}_{\text{stress}}^{3}$, corresponding to coarsegrained length ${\mathrm{\ell}}_{\text{stress}}$ and unit vector ${\overrightarrow{e}}_{i}^{n}=\left({\overrightarrow{x}}_{0}{\overrightarrow{x}}_{m}\right)/{\overrightarrow{x}}_{0}{\overrightarrow{x}}_{m}$. Herein, the stress fields are computed using ${\mathrm{\ell}}_{\text{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 ${\stackrel{~}{\sigma}}^{\text{iso}}\left(\overrightarrow{x}\right)={\sigma}^{\text{iso}}\left(\overrightarrow{x}\right)/{\sigma}_{\text{max}}^{\text{iso}}$ and outofplane shear ${\stackrel{~}{\sigma}}_{\text{xz}}\left(\overrightarrow{x}\right)={\sigma}_{\text{xz}}\left(\overrightarrow{x}\right)/{\sigma}_{\text{xz}}^{\text{max}}$, where ${\sigma}_{\text{max}}^{\text{iso}}$ and ${\sigma}_{xz}^{\text{max}}$ are the maximum values in their corresponding fields (see Figure 4a and b). We observe a high, outofplane, 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 4c shows the evolution of spatially averaged normalized isotropic stress for extruding cell $i$, ${\overline{\stackrel{~}{\sigma}}}_{i}^{\text{iso}}\left(\stackrel{~}{t}\right)={\u27e8{\sigma}^{\text{iso}}(\overrightarrow{x},\stackrel{~}{t})\u27e9}_{\overrightarrow{x}\in {\mathcal{R}}_{i}}/{\u27e8{\sigma}^{\text{iso}}(\overrightarrow{x},\stackrel{~}{t})\u27e9}_{\overrightarrow{x}\in \mathcal{R}}$, demonstrating a clear stress buildup, followed by a drop near $\stackrel{~}{t}{\stackrel{~}{t}}_{e}=0$, as a cell detaches the substrate and loses contact with other cells, where ${\stackrel{~}{t}}_{e}$ is detected by our stressindependent criterion, $\mathcal{R}={\bigcup}_{i=1}^{N}{\mathcal{R}}_{i}$ and ${\mathcal{R}}_{i}$ is the domain associated with cell i, ${\mathcal{R}}_{i}:=\{\overrightarrow{x}{\varphi}_{i}\left(\overrightarrow{x}\right)\ge 0.5\}$.
Similarly, Figure 4d displays the spatially averaged normalized outofplane shear stress, ${\overline{\stackrel{~}{\sigma}}}_{xz}^{i}\left(\stackrel{~}{t}\right)={\u27e8{\sigma}_{xz}(\overrightarrow{x},\stackrel{~}{t})\u27e9}_{\overrightarrow{x}\in {\mathcal{R}}_{i}}/{\u27e8{\sigma}_{xz}(\overrightarrow{x},\stackrel{~}{t})\u27e9}_{\overrightarrow{x}\in \mathcal{R}}$, 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 nonextruding 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 outofplane shear stress has parallels with the phenomenon of plithotaxis, where it was shown that cells collectively migrate along the orientation of the minimal inplane intercellular shear stress (Tambe et al., 2011). In this context, based on the association of cell extrusion events with regions of high outofplane shear stress, we conjecture that high shear stress concentration hinders collective cell migration with cell extrusion providing a mechanism to reestablish the status quo. This is also consistent with the observation we made earlier about large oscillations in ${\overline{\stackrel{~}{\sigma}}}_{xz}^{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 inplane and outofplane mechanical stress patterns in the cell assembly. The relative strength of cell–cell to cell–substrate adhesion, $\mathrm{\Omega}$, 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 ${\stackrel{~}{d}}_{\text{min}}^{+1/2}$, $m=\u27e8{\stackrel{~}{d}}_{\text{min}}^{+1/2}\u27e9$, increases with $\mathrm{\Omega}={\omega}_{\text{cc}}/{\omega}_{\text{cw}}$ (see inset) while the peak of the probability density decreases with increasing $\mathrm{\Omega}$. At the same time, the probability of an extrusion occurring at a fivefold disclination increases with increasing $\mathrm{\Omega}$, as displayed in Figure 2e. However, nematic and hexatic order parameters do not show any clear trends with $\mathrm{\Omega}$ (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 ${\stackrel{~}{d}}_{\text{min}}^{+1/2}$ (see Figure 2d). An example for the normalized average isotropic stress field corresponding to ${\omega}_{\text{cw}}=0.0025$ and $\mathrm{\Omega}=0.4$ is shown in Figure 5c, where ${\stackrel{~}{\overline{\sigma}}}^{\text{iso}}\left(\overrightarrow{x}\right)={\overline{\sigma}}^{\text{iso}}\left(\overrightarrow{x}\right)/{\overline{\sigma}}_{\text{max}}^{\text{iso}}$, with ${\overline{\sigma}}^{\text{iso}}$ representing the average field during defect lifetime, for all nucleated defects, and ${\overline{\sigma}}_{\text{max}}^{\text{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 ${\stackrel{~}{\overline{\sigma}}}^{\text{iso}}$ and various $\mathrm{\Omega}$. Remarkably, as shown in Figure 5a, with increasing $\mathrm{\Omega}$, 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 $\mathrm{\Omega}$ is increased, as shown in Figure 5b, while no clear trend is observed for the density of halfinteger 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 $\mathrm{\Omega}$ is increased, cells collectively have a tendency towards leveraging fivefold disclinations instead of +1/2 defects for extruding an unwanted cell.
Furthermore, one may naively think that only the distance between a halfinteger nematic defect and an extrusion site is of importance. Such a view implicitly assumes the statistics we have presented (e.g., ${\stackrel{~}{d}}_{\text{min}}^{\pm 1/2}$, ${\stackrel{~}{\overline{\sigma}}}^{\text{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 threedimensional model of the collective migrationmediated 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 $\mathrm{\Omega}={\omega}_{\text{cc}}/{\omega}_{\text{cw}}$ is increased, (ii) more spread stress concentration around a +1/2 defect with increasing $\mathrm{\Omega}$, 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 outofplane 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 reestablish 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 Ecadherin 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 threedimensional 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 liquidtosolid transition, as well as probing the impact of threedimensionality and outofplane 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 ${\overrightarrow{e}}_{n}\phantom{\rule{veryverythickmathspace}{0ex}}\left(={\overrightarrow{e}}_{z}\right)={\overrightarrow{e}}_{x}\times {\overrightarrow{e}}_{y}$ and periodic boundaries in both ${\overrightarrow{e}}_{x}$ and ${\overrightarrow{e}}_{y}$, where $({\overrightarrow{e}}_{x},{\overrightarrow{e}}_{y},{\overrightarrow{e}}_{z})$ constitute the global orthonormal basis (Figure 1). Cells are initiated on a twodimensional 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 threedimensions using a phasefield, ${\varphi}_{i}={\varphi}_{i}\left(\overrightarrow{x}\right)$. The interior and exterior of cell i corresponds to ${\varphi}_{i}=1$ and ${\varphi}_{i}=0$, respectively, with a diffuse interface of length $\lambda $ connecting the two regions and the midpoint, ${\varphi}_{i}=0.5$, delineating the cell boundary. A threedimensional 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:
where $\mathcal{F}$ 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}\phantom{\rule{veryverythickmathspace}{0ex}}\left(=\left(4/3\right)\pi {R}_{0}^{3}\right)$, such that cells – each initiated with radius R_{0} – are compressible. Additionally, $\kappa $ and $\omega $ capture repulsion and adhesion between cell–cell (subscript $cc$) and cell–substrate (subscript $cw$), respectively. Moreover, $\gamma $ sets the cell stiffness and µ captures cell compressibility and ${\varphi}_{w}$ denotes a static phasefield representing the substrate. This approach resolves the cellular interfaces and provides access to intercellular forces. The dynamics for field ${\varphi}_{i}$ can be defined as:
where $\mathcal{F}$ is defined in Equation (1), and ${\overrightarrow{v}}_{i}$ is the total velocity of cell $i$. To resolve the forces generated at the cellular interfaces, we consider the following overdamped dynamics for cells:
where ${\overrightarrow{t}}_{i}$ denotes traction as defined for Bayesian Inversion Stress Microscopy in Saw et al., 2017, $\xi $ is substrate friction, and ${\overrightarrow{F}}_{i}^{\text{sp}}=\alpha {\overrightarrow{p}}_{i}$ represents selfpropulsion forces due to polarity, constantly pushing the system outofequilibrium. In this vein, $\alpha $ characterizes the strength of polarity force and $\mathbf{\Pi}}^{\text{int}}=\left(\sum _{i}\left(\delta \mathcal{F}/\delta {\varphi}_{i}\right)\right)\mathbf{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 frontrear 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
where ${\theta}_{i}\in [\pi ,\pi ]$ is the angle associated with polarity vector, ${\overrightarrow{p}}_{i}=(\mathrm{cos}{\theta}_{i},\mathrm{sin}{\theta}_{i},0)$, and $\eta $ is the Gaussian white noise with zero mean, unit variance, ${D}_{r}$ is rotational diffusivity, $\mathrm{\Delta}{\theta}_{i}$ is the angle between ${\overrightarrow{p}}_{i}$ and ${\overrightarrow{t}}_{i}$, and positive constant $J$ sets the alignment time scale. It is worth noting that the selfpropulsion forces, ${\overrightarrow{F}}_{i}^{\text{sp}}$, associated with cell polarity, ${\overrightarrow{p}}_{i}$, act inplane but can induce outofplane components in force and velocity fields as a cell described by ${\varphi}_{i}\left(\overrightarrow{x}\right)$ deforms in threedimensions (see Equation 3).
We perform largescale 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 ${\omega}_{\text{cw}}\in \{0.0015,0.002,0.0025\}$ and vary the cell–substrate to cell–cell adhesion ratio in the range $\mathrm{\Omega}={\omega}_{\text{cc}}/{\omega}_{\text{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}_{\text{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 hypothesistesting 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, ${\overline{d}}_{\text{min}}={d}_{\text{min}}/{R}_{0}$ between the extrusion location and the nearest halfinteger 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}_{\text{sim}})$. Then, we quantify the minimum distance ${\overline{d}}_{\text{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.
Simulation parameters
We perform largescale 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\mu m$, velocity $\sim 20\mu m/h$, and average pressure of ∼100 Pa, measured experimentally in MDCK monolayers (Saw et al., 2017), corresponding to $\mathrm{\Delta}x\sim 0.5\mu m$, $\mathrm{\Delta}t\sim 0.1s$, and $\mathrm{\Delta}F\sim 1.5$ nN for force. In this study, we set the cell–substrate adhesion strength ${\omega}_{\text{cw}}\in \{0.0015,0.002,0.0025\}$ and vary the cell–substrate to cell–cell adhesion ratio in the range $\mathrm{\Omega}={\omega}_{\text{cc}}/{\omega}_{\text{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 ${\kappa}_{\text{cc}}=0.5$, ${\kappa}_{\text{cw}}=0.15$, $\xi =1$, $\alpha =0.05$, $\lambda =3$, $\mu =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):
where ${\overline{A}}_{z}$ is the average area of cells with $z$ neighbors and $\overline{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.
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:
Outofplane shear, ${\sigma}_{xz}$, for nonextruding cells
The fluctuations in outofplane shear, ${\overline{\stackrel{~}{\sigma}}}_{xz}^{i}\left(\stackrel{~}{t}\right)={\u27e8{\sigma}_{xz}(\overrightarrow{x},\stackrel{~}{t})\u27e9}_{\overrightarrow{x}\in {\mathcal{R}}_{i}}/{\u27e8{\sigma}_{xz}(\overrightarrow{x},\stackrel{~}{t})\u27e9}_{\overrightarrow{x}\in \mathcal{R}}$ for an extruding cell $i$ normalized by the maximum outofplane shear for all nonextruding cells in the same temporal window, as shown in Appendix 1—figure 4.
Halfinteger defect statistics
We have computed the lifetime for a halfinteger 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 halfinteger 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.
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), ${\omega}_{cw}$, affects the extrusion rate. The results show that increasing cell–substrate adhesion leads to less extrusion events, while the ratio $\mathrm{\Omega}={\omega}_{cc}/{\omega}_{cw}$ does not seem to play a significant role on the likelihood of an extrusion event occurring.
Energy spectra
We calculated energy spectra for different cellcell adhesion strengths, which suggests different powerlaw regimes, as shown in Appendix 1—figure 9. The kinetic energy spectrum, ${\widehat{E}}_{v}=\frac{1}{2}\u27e8{\widehat{v}}_{i}\left(k\right){\widehat{v}}_{i}\left(k\right)\u27e9$, where ${\widehat{v}}_{i}$ is the Fourier transforms of the velocity field, and ${\stackrel{~}{E}}_{v}={\widehat{E}}_{v}\left(k\right)/{\widehat{E}}_{v}^{\text{max}}\left(k\right)$. Furthermore, $\stackrel{~}{k}=k/\left(2\pi /{R}_{0}\right)$.
$p$ atic order
We computed the $p$ atic order parameter, associated with a liquid crystal exhibiting $p$ fold rotational symmetry (Nelson and Halperin, 1979), ${\psi}_{p}^{i}=\frac{1}{{N}_{n}^{i}}{\sum}_{j}^{{N}_{n}^{i}}\mathrm{exp}(pi{\theta}_{ij})$, where ${N}_{n}^{i}$ is the number of neighbors for cell $i$, ${\theta}_{ij}$ is the angle between vector ${\overrightarrow{r}}_{ij}$, connecting cell $i$ and neighboring cell $j$, and ${\overrightarrow{e}}_{x}$. Lastly, $p=2$ for nematic and $p=6$ for hexatic phases. The mean of the absolute value, $\overline{{\psi}_{2}}$ and $\overline{{\psi}_{6}}$ for various $\mathrm{\Omega}={\omega}_{cc}/{\omega}_{cw}$ is shown in Appendix 1—figure 10.
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, $\overline{{\psi}_{6}}$ displayed in Appendix 1—figure 11b, plateauing at $\overline{{\psi}_{6}}=1$ indicative of perfect hexatic order.
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

Observations on the social behaviour of cells in tissue cultureExperimental Cell Research 6:293–306.https://doi.org/10.1016/00144827(54)901767

BookPhysical Models of Cell Motility. Biological and Medical Physics, Biomedical EngineeringCham: Springer International Publishing.https://doi.org/10.1007/9783319244488

Selfpropelled rods: Insights and perspectives for active matterAnnual Review of Condensed Matter Physics 11:441–466.https://doi.org/10.1146/annurevconmatphys031119050611

Turbulent dynamics of epithelial cell culturesPhysical Review Letters 120:208101.https://doi.org/10.1103/PhysRevLett.120.208101

Forces driving epithelial wound healingNature Physics 10:683–690.https://doi.org/10.1038/nphys3040

Free energy of a nonuniform systemI. Interfacial Free Energy. Journal of Chemical Physics 28:258–267.https://doi.org/10.1063/1.1744102

Physical models of collective cell motility: From cell to tissueJournal of Physics D 50:11.https://doi.org/10.1088/13616463/aa56fe

Why we need mechanics to understand animal regenerationDevelopmental Biology 433:155–165.https://doi.org/10.1016/j.ydbio.2017.09.021

A micromechanical description of granular material behaviorJournal of Applied Mechanics 48:339–344.https://doi.org/10.1115/1.3157619

Free energies of isolated five and sevenfold disclinations in hexatic membranesPhysical Review. E, Statistical Physics, Plasmas, Fluids, and Related Interdisciplinary Topics 53:2551–2559.https://doi.org/10.1103/physreve.53.2551

BookThe Physics of Liquid Crystals. Number 83 in The International Series of Monographs on PhysicsClarendon Press.

Single and collective cell migration: The mechanics of adhesionsMolecular Biology of the Cell 28:1833–1846.https://doi.org/10.1091/mbc.E17030134

DefectMediated morphologies in growing cell coloniesPhysical Review Letters 117:048102.https://doi.org/10.1103/PhysRevLett.117.048102

Spontaneous shear flow in confined cellular nematicsNature Physics 14:728–732.https://doi.org/10.1038/s4156701800997

Tethering, crumpling, and melting transitions in hexatic membranesEurophysics Letters 13:441–446.https://doi.org/10.1209/02955075/13/5/011

Theory of defectmediated morphogenesisScience Advances 8:15.https://doi.org/10.1126/sciadv.abk2712

Characterization of the interface between normal and transformed epithelial cellsNature Cell Biology 11:460–467.https://doi.org/10.1038/ncb1853

Mechanical regulation of epithelial tissue homeostasisPhysical Review X 11:031029.https://doi.org/10.1103/PhysRevX.11.031029

Lewis’ law revisited: The role of anisotropy in sizetopology correlationsNew Journal of Physics 16:015024.https://doi.org/10.1088/13672630/16/1/015024

Mechanobiology of collective cell behavioursNature Reviews. Molecular Cell Biology 18:743–757.https://doi.org/10.1038/nrm.2017.98

BookTheory of elasticityIn: Landau LD, Lifshitz EM, editors. Number Vol. 7 in Course of Theoretical Physics. Pergamon Press. pp. 7–44.

Tissue crowding induces caspasedependent competition for spaceCurrent Biology 26:670–677.https://doi.org/10.1016/j.cub.2015.12.072

Solidliquid transition of deformable and overlapping active particlesPhysical Review Letters 125:038003.https://doi.org/10.1103/PhysRevLett.125.038003

Chiral active hexatics: Giant number fluctuations, waves, and destruction of orderPhysical Review Letters 125:238005.https://doi.org/10.1103/PhysRevLett.125.238005

Bacteria solve the problem of crowding by moving slowlyNature Physics 17:205–210.https://doi.org/10.1038/s41567020010706

SoftwareSiavashmonfared.github.io, version swh:1:rev:544e446b755bb5949e3e2623fb18f38a98f22c7eSoftware Heritage.

Emergence of active nematic behavior in monolayers of isotropic cellsPhysical Review Letters 122:048004.https://doi.org/10.1103/PhysRevLett.122.048004

Dislocationmediated melting in two dimensionsPhysical Review B 19:2457–2484.https://doi.org/10.1103/PhysRevB.19.2457

A mechanical instability in planar epithelial monolayers leads to cell extrusionBiophysical Journal 118:2549–2560.https://doi.org/10.1016/j.bpj.2020.03.028

Hexatic phase in a model of active biological tissuesSoft Matter 16:3914–3920.https://doi.org/10.1039/d0sm00109k

Sustained oscillations of epithelial cell sheetsBiophysical Journal 117:464–478.https://doi.org/10.1016/j.bpj.2019.06.013

Morphologies of compressed active epithelial monolayersThe European Physical Journal. E, Soft Matter 44:99.https://doi.org/10.1140/epje/s1018902100094x

Defects in flexible membranes with crystalline orderPhysical Review. A, General Physics 38:1005–1018.https://doi.org/10.1103/physreva.38.1005

Tumour cell invasion: An emerging role for basal epithelial cell extrusionNature Reviews. Cancer 14:495–501.https://doi.org/10.1038/nrc3767

Collective cell guidance by cooperative intercellular forcesNature Materials 10:469–475.https://doi.org/10.1038/nmat3025

Topological turbulence in the membrane of a living cellNature Physics 16:657–662.https://doi.org/10.1038/s4156702008419

Active nematic defects and epithelial morphogenesisPhysical Review Letters 129:098102.https://doi.org/10.1103/PhysRevLett.129.098102

Topological and geometrical quantities in active cellular structuresThe Journal of Chemical Physics 150:164108.https://doi.org/10.1063/1.5085766

Multiphase field models for collective cell migrationPhysical Review. E 104:054410.https://doi.org/10.1103/PhysRevE.104.054410
Decision letter

Pierre SensReviewing Editor; Institut Curie, CNRS UMR168, France

Jonathan A CooperSenior Editor; Fred Hutchinson Cancer Center, United States
Our editorial process produces two outputs: (i) public reviews designed to be posted alongside the preprint for the benefit of readers; (ii) feedback on the manuscript for the authors, including requests for revisions, shown below. We also include an acceptance summary that explains what the editors found interesting or important about the work.
Decision letter after peer review:
Thank you for submitting your article "Mechanical Basis and Topological Routes to Cell Elimination" for consideration by eLife. Your article has been reviewed by 3 peer reviewers, and the evaluation has been overseen by a Reviewing Editor and Jonathan Cooper as the Senior Editor. The reviewers have opted to remain anonymous.
The reviewers found that the 3D phase field model you developed is a very valuable improvement over existing 2D models, but that your claims of how extrusion is linked to topological parameters was insufficiently justified.
The reviewers have discussed their reviews with one another, and the Reviewing Editor has drafted this to help you prepare a revised submission.
Essential revisions:
1) Substantiate the claim that the threedimensional phase field model is crucial for understanding cell extrusion as compared to 2D models. This should include in depth analysis of necessarily 3D parameters, such as basal and apical cell surfaces.
2) Provide a better statistics of the topological parameters and their fluctuations in normal (non extruding) cells to better assess their change during extrusion.
3) Improve the characterisation of the extrusion points, e.g. by improving the representation and discussion of the stress (Figure 4).
4) Provide a pointbypoint answer to the criticism of all three reviewers and appropriately amend the paper.
Reviewer #1 (Recommendations for the authors):
1. The authors say "The similarities between cellular systems and liquid crystals featuring both nematic order (Saw et al., 2017; Kawaguchi et al., 2017; Duclos et al., 2018; BlanchMercader et al., 2018; Tan et al., 2020; Zhang et al., 2021) and hexatic order (Pasupalak et al., 2020; Maitra et al., 2020; Hoffmann et al., 2022) with the two phases potentially coexisting (ArmengolCollado et al., 2022) and interacting provide a fresh perspective for understanding cellular processes." Do the authors also find that the two phases coexist? I find this rather difficult to understand.
2. The 1/2 defects are defined using a nematic order parameter (I guess the positions of the defects are identified as the positions of zeros in this order parameter and the sign is assigned by looking at the orientation field around that). However, I think it will make the work more selfcontained and easier to understand if the authors explicitly explain how they obtain a nematic order parameter field from a map of projected 2d cell shapes.
3. I am guessing that the 1/6 defects are defined purely using the coordination number of the cells (or do the authors construct a hexatic field and use the zeros of that to identify the hexatic disclinations?). But the probability density of average coordination number from Figure 2e seems to be peaked at 5 rather than 6. Therefore, I don't understand the rationale for considering hexatic defects (or hexatic order). Also, would it be possible to define a hexatic order parameter the same way the authors (presumably) defined a nematic order parameter (and director field) and obtain the defects explicitly as singularities in this field?
4. The authors say "Figure 4(c) shows the evolution of spatially averaged normalized isotropic stress for extruding cell…. demonstrating a clear stress build up, followed by a drop near $\stackrel{~}{t}{\stackrel{~}{t}}_{e}=0$". Maybe I am missing something, but I do not see this drop. Instead, the mean (red line) increases monotonically.5. The authors say that the outofplane shear stress "prior to extrusion exhibits oscillations with large magnitudes relative to the mean field." However, since they don't present any data on the usual fluctuation of the outofplane shear stress (i.e. for nonextruding cells), the reader has no way to judge whether this is atypical. Surprisingly, the standard deviation seems to be larger at earlier times (i.e., away from the extrusion event). Could the authors compare the standard deviation of the shear stress fluctuations of extruding and nonextruding cells?
6. If the earlier point is not clarified, it is not really clear to me what new information is obtained from a 3d model of the cell layer that couldn't be obtained from a 2d phase field model. (I do understand that their criteria for identifying extrusion is only available in the 3d model, but in a 2d model, one could choose a different criterion  shrinkage of the (2d projected) area of a cell below a critical value or cell overlap for instance  which would probably lead to similar results).
7. The cells in the authors' model move in the direction of the cell polarisation. It would be useful if the authors could comment how this cell polarisation is determined.
8. In the same vein, the cell polarisation doesn't directly affect the cellular structure (i.e., there is no term involving the polarisation in the free energy in Eq. 1). Is that assumption correct?
9. The authors say "Interestingly, the association of cell extrusion events with regions of high outofplane shear stress has parallels with the phenomenon of plithotaxis, where it was shown that cells collectively migrate along the orientation of the minimal inplane intercellular shear stress (Tambe et al., 2011). In this context, based on the association of cell extrusion events with regions of high outofplane shear stress, we conjecture that high shear stress concentration hinders collective cell migration with cell extrusion providing a mechanism to reestablish the statusquo." I don't see clear evidence of high outofplane shear stress in 4b. Figure 4b has two extrusion sites, one of which certainly displays high outofplane shear stress, but the other, not so much. Could the authors quantify their claim that extrusion events are statistically associated with high outofplane shear stress?
10. The authors claim that the probability of extrusion at "nematic and hexatic disclinations " changes depending on cellcell and cellsubstrate adhesion. Is this a secondary consequence of the structure of the cell layer changing depending on those parameters (i.e., going from a predominantly nematiclike organisation to a more hexatic organsiation)?
11. If I understand correctly, the cells are extruded at 1/2 disclinations as well, which I find puzzling.
Reviewer #2 (Recommendations for the authors):
I would like to point out:
– The identification of defects, as shown in Figure 2, sensitively depends on various thresholds. The cell with 4 neighbors in Figure 2b could easily be also considered as a cell with 8 neighbors. How sensitive are the results with respect to the thresholds used?
– Why is the approach to find correlations between defects and extruding cells different for positional and orientational defects? I would find it more natural to also average over various runs and time intervals to identify such correlations for the positional defects.
– In the conclusion it is argued that negative Gaussian curvature cannot form due to the rigid substrate. This somehow indicates that the basal side is considered, for which I understand this argument. For the apical side I don't! As I assume that the experimental data in Saw et al. show the apical side, I wonder what the relation is?
Reviewer #3 (Recommendations for the authors):
This work explores the linkage between extrusion and topological defects in cell monolayers. To better understand this linkage, the authors used cuttingedge numerical simulations that were developed by some of the authors in Mueller et al. 2019. I have some major concerns regarding some theoretical analyses and interpretations of results (see major point 1 and minor points). I cannot recommend the publication of the present manuscript before these concerns are addressed.
1. Linkage between topological defects and extrusion events.
1.1. In line 104 to 108, it is suggested that the state in Figure 1a is an active nematic turbulent state based on the nucleation and dynamics of nematic defects. However, an active turbulent state is also characterized by specific statistical features of flows. Are the flows in simulation compatible with an active nematic turbulent state?
1.2. In line 117124, it is suggested that extrusions in Saw et.al. 2017 correlate with the position of + and 1/2 defects. Unless I am mistaken, their observations showed spatiotemporal correlations only for +1/2 defects. Can the authors comment on this point? If this is the case, how do the authors interpret the fact that 1/2 defects also correlate with cell extrusion in their simulations? Is the mechanism for cell extrusion in 1/2 defects of mechanical origin? Adding the average isotropic stress near 1/2 defects would help to clarify this point.
1.3. In lines 124128, the authors define a time interval around an extrusion event, but I could not understand the reason for the choice of the lower and upper values (5.625 and 0.625). How did the authors choose these values? In the Figure 4, the stress buildup occurs within a time interval of less 1 unit of time, and in Figure 8, the change in the coordination number and the area seem to occur over time intervals of 10100 units of time. Can the authors discuss the separation of these time scales? Another relevant time scale is the mean lifetime of nematic topological defects in their simulations, can discuss show the distribution of lifetimes of nematic defects and discuss how this is related to the other time scales?
1.4. The authors observed that extrusion preferentially occur near + and 1/2 defects. However, in Figure 2a, there are two extrusion events and tens of defects, which seems to indicate that many defects are nonfunctional because they do not generate an extrusion event. Is it the case? Can the authors provide an estimate of the fraction of defects that are nonfunctional? Along these lines, a question that remains unaddressed is whether an extrusion event can favor the nucleation of a pair of nematic defects. Can the authors comment on this?
1.5. In Figure 2a it seems that the typical separation between nematic defects is approximately 6 cell sizes. In the case that extrusion events occur randomly, I would have expected d_min to be below 3. However, in Figure 3a, d_min can reach 25 cell sizes. Can the author comment on these differences? Can the authors include the distributions of density of nematic defects in their simulation?
Besides, can the authors add numerical details in Methods on how were the extrusion events generated in the hypothesistesting approach? For example, were cells pulled upwards with a constant force or at a constant speed?
[Editors' note: further revisions were suggested prior to acceptance, as described below.]
Thank you for resubmitting your work entitled "Mechanical Basis and Topological Routes to Cell Elimination" for further consideration by eLife. Your revised article has been evaluated by Jonathan Cooper (Senior Editor) and a Reviewing Editor. The referees did acknowledge that the manuscript has been significantly improved, but there remain a number of serious criticism that must be answered.
Model assumptions: Many possible 3D features that could be responsible for complex 3D cell shapes are not included in the model. This needs to be explicitly discussed in the manuscript. These include:
1) You present a 3D model and claim that it is essential to understand cell extrusion, but most of the analysis is done on 2D quantities. The results regarding the isotropic and outofplane shear stress are certainly interesting but do not constitute such an explanation for the extrusion process. Furthermore, a number of assumptions of the model relate to 2D features. For instance, the polarity appears to be a 2D quantity defined by an angle (probably in the xy plane, this should be specified). This seems to exclude active outofplane relative movement between cells which could participate in the extrusion process in real systems. Is there any reason for the RHS of Eq. 3. to be a twodimensional vector or is the cell velocity not strictly in the plane? This should be clarified, and the limitations of the model clearly discussed
2) Furthermore, the magnitude of the polarity vector is set to unity, which precludes fluctuations of activity among cells, a feature that is often associated with local rearrangement processes such as extrusion. This is a strong assumption that is not discussed.
2) Extrusion has often been described as the result of different mechanical properties in the basal and apical sides on the epithelia. You say in your rebuttal that this is just one element of the 3D complexity on which you choose not to focus. It is acceptable not to include this in a model in order to focus on other effects, however, this possibility should be mentioned as a relevant process in the manuscript.
3) Scutoids. It does not appear to be clearly demonstrated that scutoids are a feature of curved epithelia only, and it is quite possible that such arrangement might be relevant to the extrusion process even in flat epithelia. While a model such as the present one that does not include this possibility is certainly valuable, the text should not give the impression that they are irrelevant.
Analysis:
4) Most of the analysis is done on projected quantities, but how the projection is obtained is not clearly explained. What does the 2D representation (such as Figure 2 for instance) exactly show? Is it a cut through the epithelia at fixed z, or some kind of maximum intensity projection? This must be specified.
5) in Figure 2a, it seems that the director field (red bars) around blue +1/2 defects correspond to 1/2 defects, and viceversa the director field around green 1/2 defects correspond to +1/2 defects. This feature is clearer for defects that are far from others. If this is indeed the case, it should be corrected.
6) In the new Figure 12, the probability density of pairwise distance between defects is presented. It is unclear whether this distance corresponds to the minimal distance between pairs of +1/2 defects or rather the distance between halfinteger defects. To compare with the results from Figure 2 the former seems more appropriate. Please comment on this point and include the former distribution if necessary.
7) Some new results presented in the revised version are confusing, or insufficiently discussed. For instance the new Figure 6 in the rebuttal. If one holds cellsubstrate adhesion fixed, the mean number of extrusions has a clear maximum at a particular cellcell adhesion. However, no other quantity that the authors present seems to show a similar maximum. Maybe the origin of this is hiding in the plot of nematic defect density for different adhesions (Figure 5 of the rebuttal). If so, that plot should be improved and both that and the Figure 6 of the rebuttal should be brought to the main text, as this seems an important feature of the problem
Discussion:
8) Threshold: In your response and in the manuscript (l.192), you seem to suggest that your method does not suffer from the arbitrariness of setting a threshold to defined extrusion events, but it appears (l.99) that extrusion also involves thresholding (of the cell vertical displacement – l.99) in your case. How would this threshold influence the results?
https://doi.org/10.7554/eLife.82435.sa1Author response
Essential revisions:
1) Substantiate the claim that the threedimensional phase field model is crucial for understanding cell extrusion as compared to 2D models. This should include in depth analysis of necessarily 3D parameters, such as basal and apical cell surfaces.
2) Provide a better statistics of the topological parameters and their fluctuations in normal (non extruding) cells to better assess their change during extrusion.
3) Improve the characterisation of the extrusion points, e.g. by improving the representation and discussion of the stress (Figure 4).
4) Provide a pointbypoint answer to the criticism of all three reviewers and appropriately amend the paper.
We thank the reviewers for constructive feedback and the Reviewing Editor for summarizing the areas for improvement. What follows is our pointbypoint response to the three reviewers. In summary:
We show that the localization of isotropic stresses in an extruding cell, a major point we make in our manuscript, is a threedimensional phenomenon and cannot be captured in twodimensional approaches (see comment 5 by reviewer 1). Furthermore, we also show the significant contrast in amplitude of outofplane shear oscillations in extruding cells relative to nonextruding cells (see comment 4 by reviewer 1). Also, we address comments about basal and apical properties in our studied model and provide additional data regarding the role of substrate properties on extrusion rate (see comment 3 by reviewer 2).
We have added more statistics to quantify halfinteger defects, including probability densities for defect life time (see comment 3 by reviewer 3), defect density (see comment 9 by reviewer 1), and distance between defects (see comment 5 by reviewer 3) to further substantiate the claims made in our manuscript.
We have clarified some misunderstanding regarding characterization of extrusion events and their links with nematic and hexatic defects (see comment 5 by reviewer 3, comment 6 by reviewer 3, comment 17 by reviewer 3, comment 6 by reviewer 3, comment 8 by reviewer 1, comment 9 by reviewer 1).
Reviewer #1 (Recommendations for the authors):
1. The authors say "The similarities between cellular systems and liquid crystals featuring both nematic order (Saw et al., 2017; Kawaguchi et al., 2017; Duclos et al., 2018; BlanchMercader et al., 2018; Tan et al., 2020; Zhang et al., 2021) and hexatic order (Pasupalak et al., 2020; Maitra et al., 2020; Hoffmann et al., 2022) with the two phases potentially coexisting (ArmengolCollado et al., 2022) and interacting provide a fresh perspective for understanding cellular processes." Do the authors also find that the two phases coexist? I find this rather difficult to understand.
The existence of bothtypes of defects in our simulations, namely: (i) halfinteger defects in a director field based on the elongation axes of cells and (ii) the associated 5fold and 7fold disclinations in hexagonal packing of the cells suggests the coexistence of nematic and hexatic phases. To show this more concretely, we have now computed the p−atic order parameter, associated with a liquid crystal exhibiting p−fold rotational symmetry [1], ${\psi}_{p}^{i}=\frac{1}{{N}_{n}^{i}}\sum _{j}^{{N}_{n}^{i}}exp({\text{pi\theta}}_{\text{ij}})$ where $N}_{n}^{i$ is the number of neighbors for cell i, θ_{ij} is the angle between vector $\overrightarrow{r}}_{\text{ij}$, connecting cell i and neighboring cell j, and $\overrightarrow{e}}_{x$. Lastly, p = 2 for nematic and p = 6 for hexatic phases. The mean of the absolute value, $\overline{\psi}}_{2$ and $\overline{\psi}}_{6$ for various Ω = ω_{cc}/ω_{cw} is now added to the Appendix and the coexistence of nematic and hexatic order are evident from these measures.
Actions taken:
Figure 1 is added to Appendix. It details the temporal evolution of the mean nematic and hexatic order parameters.
Line 244245: However, nematic and hexatic order parameters do not show any clear trends with 78 varying cellcell adhesion strength Ω (see Figure 15 in Appendix).
2. I am guessing that the 1/6 defects are defined purely using the coordination number of the cells (or do the authors construct a hexatic field and use the zeros of that to identify the hexatic disclinations?). But the probability density of average coordination number from Figure 2e seems to be peaked at 5 rather than 6. Therefore, I don't understand the rationale for considering hexatic defects (or hexatic order). Also, would it be possible to define a hexatic order parameter the same way the authors (presumably) defined a nematic order parameter (and director field) and obtain the defects explicitly as singularities in this field?
The five and sevenfold disclinations in the hexatic arrangement are identified based on their coordi88 nation numbers. We believe there is a misunderstanding here. Figure 2(e) is the average coordination number over a temporal window relative to extrusion, $\overline{z}$, for extruded cells, not all cells. The peak in Figure 2e is therefore indication of the higher probability of extrusion near fivefold disclinations. Figure 2(e) should be contrasted to Figure 5(b), which shows the coordination number for all cells and simulation time steps, z, and various relative cellcell adhesion strengths, Ω = ω_{cc}/ω_{cw}. As shown in Figure 5(b), when all cells are considered the probability of a cell having six neighbors, i.e. z = 6, is highest in all cases. We have added the following to the main text for further clarification:
Actions taken:
Figure 2. caption: The probability density of average coordination number ¯z for an extruding cell during $\stackrel{~}{t}=(t/{t}_{0})\in [{\stackrel{~}{t}}_{e}2.5,{\stackrel{~}{t}}_{e}+0.3125]$, where $\stackrel{~}{t}}_{e$ denotes extrusion time, τ_{0} = ξR_{0}/α and for varying cellsubstrate to cellcell adhesion ratios Ω.
Figure 5 caption: Probability density function for the number of neighbors z and various Ω, for all cells and simulation time steps.
Line 167176: To quantify the relation between extrusion events and the disclinations, the probability density of the coordination number of an extruding cell $({\stackrel{~}{d}}_{min}=0)$ averaged over the time interval, $\stackrel{~}{t}\in \left[{\stackrel{~}{t}}_{e}5.625,{\stackrel{~}{t}}_{e}+0.625\right],\stackrel{~}{z},$ for all the realizations is shown in Figure 2(e), clearly exhibiting a sharp peak near $\stackrel{~}{z}=5$. The coordination number is determined based on the interactions of cells (see Appendix) and this property is independent of apical or basal considerations (Kaliman et al., 2021), unlike bent epithelia. In the case of curved epithelia, a geometrical structure called scutoids can arise as a general feature with different properties associated with apical and basal sides (Galvez et al., 2018). In our setup, the asymmetric interactions of cells with apical and basal sides are captured by varying the strength of cellsubstrate adhesion. In our simulations, increasing cellsubstrate adhesion leads to lower extrusion events (see Figure 13 in Appendix).
3. The authors say "Figure 4(c) shows the evolution of spatially averaged normalized isotropic stress for extruding cell…. demonstrating a clear stress build up, followed by a drop near $\stackrel{~}{t}{\stackrel{~}{t}}_{e}=0$. Maybe I am missing something, but I do not see this drop. Instead, the mean (red line) increases monotonically.
As each extruding cell detaches from the substrate and loses contact with the surrounding cells, the forces due to its interactions with other cells vanish and a drop in isotropic stress is observed. To illustrate this, Author response image 1 displays one clear example with extended temporal window relative to Figure 4(c) in the main text. For further clarification we have now added explicitly the following to the main text:
Actions taken:
Line 218218: as a cell detaches the substrate and loses contact with other cells
4. The authors say that the outofplane shear stress "prior to extrusion exhibits oscillations with large magnitudes relative to the mean field." However, since they don't present any data on the usual fluctuation of the outofplane shear stress (i.e. for nonextruding cells), the reader has no way to judge whether this is atypical. Surprisingly, the standard deviation seems to be larger at earlier times (i.e., away from the extrusion event). Could the authors compare the standard deviation of the shear stress fluctuations of extruding and nonextruding cells?
To show this more clearly we have plotted the evolution of outofplane shear stress, $\overline{\stackrel{~}{\sigma}}}_{\text{xz}}^{i$ for extruding cell i normalized by the maximum value of the shear stress, $\stackrel{~}{\overline{\stackrel{~}{\sigma}}}}_{\text{xz}$, for all nonextruding cells during the same temporal window, as shown in Appendix 1—figure 9. The yaxis displays the contrast, more than an order of magnitude enhancement of the fluctuations for extruding cells, even when normalized by the maximum value within the same temporal window and for all nonextruding cells. We note that this further emphasizes the importance of the access to 3D stress fields, which is not possible with 2D models and current experimental characterization of the traction forces on 2D substrates.
Actions taken:
We have added Appendix 1—figure 9 to the Appendix. It shows the temporal evolution of the outofplane shear stress for an extruding cell i, normalized by the maximum outofplane shear of nonextruding cells, within the same temporal window.
Line 224225: ...a stark departure from their nonextruding counterparts (see Figure 9 in Appendix).
5. If the earlier point is not clarified, it is not really clear to me what new information is obtained from a 3d model of the cell layer that couldn't be obtained from a 2d phase field model. (I do understand that their criteria for identifying extrusion is only available in the 3d model, but in a 2d model, one could choose a different criterion  shrinkage of the (2d projected) area of a cell below a critical value or cell overlap for instance  which would probably lead to similar results).
We are confident that the previous point is now clarified.
Nevertheless, we believe there is a misunderstanding here as we do not use any criterion – critical values or thresholds – to induce an extrusion event. Rather, an extrusion event generically emerges from our simulations and depends on the cellcell and cellsubstrate interactions.
We respectfully do not agree that similar results can be observed in a 2D model. As clarified in the response to the previous comment we clearly have shown the importance of having access to three dimensional fields and we have discussed thoroughly mechanisms such as buckling of hexatic defects which we would not be able to identify if we did not have access to threedimensional fields. Furthermore, the isotropic stress buildup shown in Figure 4(c) in the main text is a threedimensional quantity, i.e. σ^{iso} = (1/3)(σ_{xx} + σ_{yy} + σ_{zz}). A similar buildup is not observed in twodimensions, i.e. σ^{iso} = (1/2)(σ_{xx} + σ_{yy}), as shown in Author response image 2. Moreover, a decrease in cell area does not necessarily mean an extrusion event as it is wellestablished that projected cell area can indeed significantly fluctuate in 158 cell monolayers [2, 3]. It is then immediately clear that setting an arbitrary threshold based on the cell area reduction can lead to false extrusion events in 2D models.
6. The cells in the authors' model move in the direction of the cell polarisation. It would be useful if the authors could comment how this cell polarisation is determined.
This is explained in Materials and methods. In particular Eq. (4) outlines the dynamics of cell polarisation and immediately after this equation the polarity vector for cell i, $\overrightarrow{p}}_{i$ is defined. We have added a sentence in the main text to clarify this.
Actions taken:
Line 9696: see Methods for polarisation dynamics.
7. In the same vein, the cell polarisation doesn't directly affect the cellular structure (i.e., there is no term involving the polarisation in the free energy in Eq. 1). Is that assumption correct?
Yes, the cell polarisation is not part of the free energy of a cell but it drives the system out of equilibrium and can influence the modes of collective selforganization of the cells and their interactions. Please see Materials and methods section for justification of the choice of polarisation mechanism based on previous (qualitative and quantitative) comparison with existing experimental results.
8. The authors say "Interestingly, the association of cell extrusion events with regions of high outofplane shear stress has parallels with the phenomenon of plithotaxis, where it was shown that cells collectively migrate along the orientation of the minimal inplane intercellular shear stress (Tambe et al., 2011). In this context, based on the association of cell extrusion events with regions of high outofplane shear stress, we conjecture that high shear stress concentration hinders collective cell migration with cell extrusion providing a mechanism to reestablish the statusquo." I don't see clear evidence of high outofplane shear stress in 4b. Figure 4b has two extrusion sites, one of which certainly displays high outofplane shear stress, but the other, not so much. Could the authors quantify their claim that extrusion events are statistically associated with high outofplane shear stress?
We believe this is now fully clarified based on our response to comment 4. Figure 4(b) is just a snapshot in time. The fields (and the colormap) are normalized and the ”hotspots” are only at the two marked extrusion sites. Since there is a time gap between these two extrusion events, one extrusion site may appear to have a higher outofplane shear relative to the other site. However, compared to the whole domain, both sites have much larger values (please note the colorbar).
The statistics backing our claim is clearly shown in Figure 4(d). Please note the scaling of the yaxis (×10^{4}) for the normalized stress. Furthermore, please see our response to comment 4 above.
9. The authors claim that the probability of extrusion at "nematic and hexatic disclinations " changes depending on cellcell and cellsubstrate adhesion. Is this a secondary consequence of the structure of the cell layer changing depending on those parameters (i.e., going from a predominantly nematiclike organisation to a more hexatic organisation)?
In Appendix 1—figure 11(b) we show that increasing relative cellcell to cellsubstrate adhesion strength Ω results in a higher probability for a cell to be at a fivefold disinclination. For completion, we have also included the halfinteger defect densities for various Ω in the appendices. As shown, there are no clear trends with Ω. The density is defined as the number of ±1/2 defects normalized by the domain area, i.e. a = L_{x} × L_{y}. In addition, see our response to comment 1, above, where we compute mean nematic $\left({\overline{\psi}}_{2}\right)$ and hexatic $\left({\overline{\psi}}_{6}\right)$ order parameters, where we also do not observe significant variation of ${\overline{\psi}}_{2}$ and ${\overline{\psi}}_{6}$ with varying the cellcell adhesion strength. Instead, we show in Appendix 1—figure 11(a) that changing cellcell adhesion strength results in significant changes in the distribution of isotropic stress. This is the 208 potential explanation which we discuss thoroughly in the main text.
Actions taken:
We have added Appendix 1—figure 11 to the Appendix. It displays the probability density of halfinteger nematic defect density for different Ω = ω_{cc}/ω_{cw}.
Line 260261: while no clear trend is observed for the density of halfinteger defects (see Figure 11 in Appendix).
10. If I understand correctly, the cells are extruded at 1/2 disclinations as well, which I find puzzling.
Yes, that is correct and consistent with Saw et al., Nature, 2017 where they observed extrusion events associated with both +1/2 and −1/2 defects, with a relatively stronger association for +1/2 defects. This text is from Saw et al., Nature, 2017: “We found that extrusion events were strongly correlated with the positions of a subset of +1/2 defects (and less so with 1/2 defects) (Figure 1f, Extended Data Figure 1eh, Methods)”.
Actions taken:
Line 7374: …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).
Reviewer #2 (Recommendations for the authors):
I would like to point out:
– The identification of defects, as shown in Figure 2, sensitively depends on various thresholds. The cell with 4 neighbors in Figure 2b could easily be also considered as a cell with 8 neighbors. How sensitive are the results with respect to the thresholds used?
We determine the number of neighbors based on how many cells a cell is interacting with. We specifically chose this criterion as opposed to typical voronoi tessellation since during our simulations and in regions with high number of extrusion events in temporal proximity the confluency can be lost. In that case, methods based in voronoi tessellation would overcount the neighbors. To determine interactions among cells, we use the overlap between the phasefields with a cutoff value as detailed now in the appendices.
Please note that there is only one threshold which determines the extent of the phasefield domain.
Actions taken:
We have added the following to Appendix:
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} ∩ {ϕ_{j}ϕ_{j} > 0.25} ≠ ∅ (1)
– Why is the approach to find correlations between defects and extruding cells different for positional and orientational defects? I would find it more natural to also average over various runs and time intervals to identify such correlations for the positional defects.
We believe this is a misunderstanding. In both cases, we look at the probability densities associated with relevant parameters. Since the halfinteger defect cores associated with an extrusion event is not perfectly centered on an extruding cell, we use $\stackrel{~}{d}}_{min$ to characterize this link. In the case of hexatic disclinations, it depends on the number of neighbors a cell is interacting with and thus an extruding cell could itself be a fivefold disclination for which ${\stackrel{~}{d}}_{min}=0$. Furthermore, similar to halfinteger nematic defects (Figure 2(c)(d)), the hextic defects in Figure 2(e) also correspond to four realization of the parametric space we probed, as also highlighted in the caption for Figure 2.
Actions taken:
Line 167176: To quantify the relation between extrusion events and the disclinations, the probability density of the coordination number of an extruding cell $({\stackrel{~}{d}}_{min}=0)$ averaged over the time interval, $\stackrel{~}{t}\in \left[{\stackrel{~}{t}}_{e}5.625,{\stackrel{~}{t}}_{e}+0.625\right],\stackrel{~}{z},$ for all the realizations is shown in Figure 2(e), clearly exhibiting a sharp peak near $\stackrel{~}{z}=5$. The coordination number is determined based on the interactions of cells (see Appendix) and this property is independent of apical or basal considerations (Kaliman et al., 2021), unlike bent epithelia. In the case of curved epithelia, a geometrical structure called scutoids can arise as a general feature with different properties associated with apical and basal sides (Galvez et al., 2018). In our setup, the asymmetric interactions of cells with apical and basal sides are captured by varying the strength of cellsubstrate adhesion. In our simulations, increasing cellsubstrate adhesion leads to lower extrusion events (see Figure 13 in Appendix).
– In the conclusion it is argued that negative Gaussian curvature cannot form due to the rigid substrate. This somehow indicates that the basal side is considered, for which I understand this argument. For the apical side I don't! As I assume that the experimental data in Saw et al. show the apical side, I wonder what the relation is?
In our work, the cells are on a rigid substrate – that would be the basal side. The apical side is free and cells are allowed to detach from the substrate and move outofplane. Furthermore, we only simulate a monolayer. Thus, our arguments about hexatic defects buckling and creating positive and negative Gaussian curvatures are only relevant to the basal side. The experimental data analyses in Saw et al. 2017 is based on top view of the monolayer, while all the stress measurements are based on the traction forces on the basal side. There is no reason to believe in such a setting, a flat cell layer on a relatively rigid monolayer, there will be a significant difference in projected cell shape, whether viewed apically or basally. This is further established to be the case experimentally for a flat monolayer of epithelial cells [4]. However, this may not be the case in curved epithelia, as shown by GomezGalvez et al., 2018, where Scutoids could become relevant.
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 cellsubstrate interactions. We have added a new phase diagram showing how changing cellsubstrate adhesion (basal), ω_{cw}, affects the extrusion rate.
Actions taken:
We have added Figure 6 to Appendix. It displays the mean and standard deviation of the number of extruded cells for varying cellsubstrate adhesion strength, ω_{cw}, and relative cellcell adhesion strength, Ω = ω_{cc}/ω_{cw}. The results correspond to four realizations of this parametric space.
Line 170172: The coordination number is determined based on the interactions of cells (see Appendix) and this property is independent of apical or basal considerations (Kaliman et al., 2021), unlike bent epithelia.
Line 174176: In our setup, the asymmetric interactions of cells with apical and basal sides are captured by varying the strength of cellsubstrate adhesion. In our simulations, increasing cell substrate adhesion leads to lower extrusion events (see Figure 13 in Appendix).
Reviewer #3 (Recommendations for the authors):
This work explores the linkage between extrusion and topological defects in cell monolayers. To better understand this linkage, the authors used cuttingedge numerical simulations that were developed by some of the authors in Mueller et al. 2019. I have some major concerns regarding some theoretical analyses and interpretations of results (see major point 1 and minor points). I cannot recommend the publication of the present manuscript before these concerns are addressed.
1. Linkage between topological defects and extrusion events.
1.1. In line 104 to 108, it is suggested that the state in Figure 1a is an active nematic turbulent state based on the nucleation and dynamics of nematic defects. However, an active turbulent state is also characterized by specific statistical features of flows. Are the flows in simulation compatible with an active nematic turbulent state?
Since characterizing the flow features is not the focus of the current study and could lead to deviation from the main message, we now change the terminology to “defect chaos”. Nevertheless, we acknowledge the referee’s suggestion and we have now performed calculation of energy spectra for different cellcell adhesion strengths, which suggests different powerlaw regimes (see Appendix 1figure 14). However, since there is not large enough span of wave numbers to concretely conclude a powerlaw scaling, we prefer to defer the analyses of the flow to future studies.
Actions taken:
Line 110110: (see Figure 14 in Appendix for energy spectra characterization).
We have added Appendix 1figure 14 to the Appendix. It characterizes energy spectra for various relative cellcell adhesion strengths, Ω = ω_{cc}/ω_{cw}.
1.2. In line 117124, it is suggested that extrusions in Saw et.al. 2017 correlate with the position of + and 1/2 defects. Unless I am mistaken, their observations showed spatiotemporal correlations only for +1/2 defects. Can the authors comment on this point? If this is the case, how do the authors interpret the fact that 1/2 defects also correlate with cell extrusion in their simulations? Is the mechanism for cell extrusion in 1/2 defects of mechanical origin? Adding the average isotropic stress near 1/2 defects would help to clarify this point.
Yes, that is correct and consistent with Saw et al., Nature, 2017 where they observed extrusion events associated with both +1/2 and −1/2 defects, with a relatively stronger association for +1/2 defects. This text is from Saw et al., Nature, 2017: “We found that extrusion events were strongly correlated with the positions of a subset of +1/2 defects (and less so with 1/2 defects) (Figure 1f, Extended Data Figure 1eh, Methods)”. The mechanism for extrusion close to −1/2 defects in experiments is suggested to be related to the compression patterns that form not exactly at the defect core, but at a distance from it, as is evident from the average isotropic stress measured in experiments and in the phase field model.
Actions taken:
Line 7374: …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).
1.3. In lines 124128, the authors define a time interval around an extrusion event, but I could not understand the reason for the choice of the lower and upper values (5.625 and 0.625). How did the authors choose these values? In the Figure 4, the stress buildup occurs within a time interval of less 1 unit of time, and in Figure 8, the change in the coordination number and the area seem to occur over time intervals of 10100 units of time. Can the authors discuss the separation of these time scales? Another relevant time scale is the mean lifetime of nematic topological defects in their simulations, can discuss show the distribution of lifetimes of nematic defects and discuss how this is related to the other time scales?
The time axis for Figure 4(c)4(d) are normalized, i.e. t˜= t/τ_{0} where τ_{0} = ξR_{0}/α, ξ is friction, R_{0} is the initial cell radius and α is the strength of selfpropulsion forces. However, the time axis for Figure 8 is not normalized but it does correspond to exactly the same temporal window. We have now made sure that this is represented consistently to show that there are no separation of time scales. This temporal window is based on the first moment of the defect lifetime distributions, which we have now added for completion (see Appendix 1figure 10). As shown, the first moment of defect lifetimes ⟨t˜^{±1/2}⟩ ≈ 5.69.
Actions taken:
We have updated Figure 8 with normalized time axes, t˜.
Line 130132: This temporal window is chosen based on the first moment of a defect’s lifetime distribution (see Figure 10 in the Appendix).
We have added Appendix 1figure 10 to the Appendix. It displays the probability density of halfinteger defect lifetimes.
1.4. The authors observed that extrusion preferentially occur near + and 1/2 defects. However, in Figure 2a, there are two extrusion events and tens of defects, which seems to indicate that many defects are nonfunctional because they do not generate an extrusion event. Is it the case? Can the authors provide an estimate of the fraction of defects that are nonfunctional? Along these lines, a question that remains unaddressed is whether an extrusion event can favor the nucleation of a pair of nematic defects. Can the authors comment on this?
We disagree with the view that a defect not in the vicinity of an extrusion event is “nonfunctional”, keeping in mind that Figure 2(a) is just a snapshot in time for one of the cases we considered. We know well that these defects are interacting, attracted and repelled by each other. As a consequence, those interactions affect the modes of collective selforganization of the cells and consequently their mechanical interactions.
By performing a laser ablation induced cell removal, it was shown experimentally in Saw et al. 2017 that extrusion does not favor the nucleation of a pair of nematic defects. We have clarified this point in the revised manuscript.
Actions taken:
Line 139142: 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).
1.5. In Figure 2a it seems that the typical separation between nematic defects is approximately 6 cell sizes. In the case that extrusion events occur randomly, I would have expected d_min to be below 3. However, in Figure 3a, d_min can reach 25 cell sizes. Can the author comment on these differences? Can the authors include the distributions of density of nematic defects in their simulation?
Respectfully, there are a number of misunderstandings here, which we address one by one: (i) the typical distance between two nematic defects is much larger than a typical distance between an extrusion event and a halfinteger nematic defect, (ii) determining what is a typical distance requires statistics and cannot be determined from one simulation snapshot as suggested here re. Figure 2(a). We now provide this data for completeness. (iii) A typical $\overline{d}}_{min$ based on the probability density functions (PDF) shown in Figure 2(c) and Figure 2(d) is ≈ 1.5 cell size. This is based on the peak of these PDFs. (iv) Figure 3(a) corresponds to a specific hypothesis that we outlined and falsified in the text regarding the randomness of the proximity of extrusion sites and halfinteger nematic defects. The case where $\overline{d}}_{min$ reaches 25 cell sizes corresponds to the case where extrusions were generated using a poisson point process. In the same Figure 3(a), the results corresponding to our simulations (solid black line) peaks near 1.5 cell size. Hope these points clarify any misunderstanding of the data we present in our manuscript.
We have also added the density of defects in the appendices as requested.
Actions taken:
See our response to reviewer 1, comment 9 for the density of defects (Figure 5).
Line 136138: …at a much smaller distance relative to a typical distance between two defects (see Figure 12 in Appendix).
We have added Figure 10 to Appendix. It characterizes the distance between halfinteger nematic defects during the simulations.
Besides, can the authors add numerical details in Methods on how were the extrusion events generated in the hypothesistesting approach? For example, were cells pulled upwards with a constant force or at a constant speed?
The goal of this hypothesis testing approach is to show that the correlation between +1/2 defects and the extrusion events is not random. To provide a null hypothesis we generated extrusion event locations. Specifically, we use a poisson point process to randomly generate (x,y) coordinate of an extrusion location. Furthermore, we use a uniform distribution to place this random extrusion event at a random time step. Then, we carry out our analysis based on the nematic defects in that simulation time frame. Thus, there are no physical criterion for inducing an extrusion, e.g. constant force or speed. We note that the same approach is used in the experiments for hypothesis testing to show that the correlations observed between defects and extrusions are significant compared to defects correlations with random point within the tissue.
Actions taken:
Line 149149: …randomly generate positions for extrusion events a…
[Editors' note: further revisions were suggested prior to acceptance, as described below.]
The referees did acknowledge that the manuscript has been significantly improved, but there remain a number of serious criticism that must be answered first.
Model assumptions: Many possible 3D features that could be responsible for complex 3D cell shapes are not included in the model. This needs to be explicitly discussed in the manuscript. These include:
1) You present a 3D model and claim that it is essential to understand cell extrusion, but most of the analysis is done on 2D quantities. The results regarding the isotropic and outofplane shear stress are certainly interesting but do not constitute such an explanation for the extrusion process. Furthermore, a number of assumptions of the model relate to 2D features. For instance, the polarity appears to be a 2D quantity defined by an angle (probably in the xy plane, this should be specified). This seems to exclude active outofplane relative movement between cells which could participate in the extrusion process in real systems. Is there any reason for the RHS of Eq. 3. to be a twodimensional vector or is the cell velocity not strictly in the plane? This should be clarified, and the limitations of the model clearly discussed.
Polarity in our current model is frontrear polarity and it only acts in plane. However, the inplane active selfpropulsion force associated with cell polarity can induce forces and velocities with nonzero outofplane components, complementing forces due to other interactions e.g. to cellcell adhesion, as a cell described by ${\varnothing}_{i}\left(\overrightarrow{\chi}\right)$ deforms in threedimensions. The velocity field is a threedimensional vector field, as shown in Figure 1(b) and thus includes outofplane components. Furthermore, we have established before (see our response to Reviewer 1, comment 5 of our previously submitted rebuttal), we define pressure as a threedimensional quantity $p=(\frac{1}{3})({\sigma}_{\text{xx}}+{\sigma}_{\text{yy}}+{\sigma}_{\text{zz}})$ rather than a twodimensional one $p=(\frac{1}{2})({\sigma}_{\text{xx}}+{\sigma}_{\text{yy}})$. These points have been clarified in the text.
Actions taken:
Line 353356: It is worth noting that the selfpropulsion forces, $\overrightarrow{F}}_{i}^{\text{sp}$, associated with cell polarity, $\overrightarrow{p}}_{i$, 50 acts inplane but can induce outofplane components in force and velocity fields as a cell described 51 by ${\varnothing}_{i}\left(\overrightarrow{\chi}\right)$ deforms in threedimensions (see Equation (3)).
Line 351: ${\overrightarrow{p}}_{i}=(\mathrm{cos}{\theta}_{i},\mathrm{sin}{\theta}_{i}.0)$
Line 96: frontrear
Line 347: frontrear
2) Furthermore, the magnitude of the polarity vector is set to unity, which precludes fluctuations of activity among cells, a feature that is often associated with local rearrangement processes such as extrusion. This is a strong assumption that is not discussed.
As discussed in the text, the activity in our model is due to selfpropulsion forces $\left({\overrightarrow{F}}_{i}^{\text{sp}}\right)$ associated with cell polarity $\overrightarrow{p}}_{i$. The dynamics of polarisation is given in Equation (4) which includes a Gaussian white noise (η) and a rotational diffusivity constant (D_{r}) leading to fluctuations in the direction of polarity and thus selfpropulsion forces via $\overrightarrow{F}}_{i}^{\text{sp}}=\alpha {\overrightarrow{p}}_{i$, where α represents the strength of cell motility.
Actions taken:
Line 353356: It is worth noting that the selfpropulsion forces, $\overrightarrow{F}}_{i}^{\text{sp}$, associated with cell polarity, $\overrightarrow{p}}_{i$, acts inplane but can induce outofplane components in force and velocity fields as a cell described by ${\varnothing}_{i}\left(\overrightarrow{\chi}\right)$ deforms in threedimensions (see Equation (3)).
3) Extrusion has often been described as the result of different mechanical properties in the basal and apical sides on the epithelia. You say in your rebuttal that this is just one element of the 3D complexity on which you choose not to focus. It is acceptable not to include this in a model in order to focus on other effects, however, this possibility should be mentioned as a relevant process in the manuscript.
We now have explicitly mentioned this possibility as a topic of interest for future studies.
Actions taken:
Line 277278: 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.
4) Scutoids. It does not appear to be clearly demonstrated that scutoids are a feature of curved epithelia only, and it is quite possible that such arrangement might be relevant to the extrusion process even in flat epithelia. While a model such as the present one that does not include this possibility is certainly valuable, the text should not give the impression that they are irrelevant.
We have updated our description of scutoids in the text, first suggested by one of the reviewers to be included in our manuscript, and have simplified it as to avoid any impression about its relevance or irrelevance in flat geometry.
Actions taken:
Line 173: geometrical structures called scutoids that have been identified in curved epithelial tubes (GómezGálvez et al., 2018).
Analysis:
5) Most of the analysis is done on projected quantities, but how the projection is obtained is not clearly explained. What does the 2D representation (such as Figure 2 for instance) exactly show? Is it a cut through the epithelia at fixed z, or some kind of maximum intensity projection? This must be specified.
The projected fields (Figures 2(a), 4(a)(b), 5(c)) correspond to the basal side, i.e. the layer immediately on top of the substrate. This choice is based on (a) most of the experimental analyses are typically done on the basal side (see e.g. Balasubramaniam et al., 2021 [1]) and (b) we are interested in the interplay of cellsubstrate and cellcell adhesions.
Actions taken:
Figure 4 caption: projected into xy− plane (z = 0, i.e. the basal side)
Figure 5 caption: projected into xy− plane with z = 0 i.e. the basal side
Line 107: (z = 0, i.e. the basal side)
6) in Figure 2a, it seems that the director field (red bars) around blue +1/2 defects correspond to 1/2 defects, and viceversa the director field around green 1/2 defects correspond to +1/2 defects. This feature is clearer for defects that are far from others. If this is indeed the case, it should be corrected.
We thank the reviewer for noticing this. We have updated the figure with correct symbols/legend.
Actions taken:
See updated Figure 2(a).
7) In the new Figure 12, the probability density of pairwise distance between defects is presented. It is unclear whether this distance corresponds to the minimal distance between pairs of +1/2 defects or rather the distance between halfinteger defects. To compare with the results from Figure 2 the former seems more appropriate. Please comment on this point and include the former distribution if necessary.
This is explained in the Appendix (lines 619621). Figure 12 shows the probability density of distance between halfinteger defects, specifically the pairwise distance between +1/2 and +1/2 defects (Figure 12(a)) and −1/2 and −1/2 defects (Figure 12(b)). We have also added an additional plot for the pairwise distance between +1/2 and −1/2 defects displaying a similar behavior (Figure 12(c)).
Actions taken:
Figure 12(c) has been added.
Figure 12 in the Appendix caption: 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.
8) Some new results presented in the revised version are confusing, or insufficiently discussed. For instance the new Figure 6 in the rebuttal. If one holds cellsubstrate adhesion fixed, the mean number of extrusions has a clear maximum at a particular cellcell adhesion. However, no other quantity that the authors present seems to show a similar maximum. Maybe the origin of this is hiding in the plot of nematic defect density for different adhesions (Figure 5 of the rebuttal). If so, that plot should be improved and both that and the Figure 6 of the rebuttal should be brought to the main text, as this seems an important feature of the problem
We thank the reviewer for their thoughts on the newly presented figures. The mentioned maximum is only present for one value of cellsubstrate adhesion (ω_{cw} = 0.002). Furthermore, this maximum becomes less pronounced considering the standard deviation in extrusions as shown in Figure 6(b). Lastly, the probability densities shown in Figure 5 of the rebuttal correspond to data sorted by the relative cellcell adhesion Ω and not cellsubstrate adhesion ω_{cw}.
Discussion:
9) Threshold: In your response and in the manuscript (l.192), you seem to suggest that your method does not suffer from the arbitrariness of setting a threshold to defined extrusion events, but it appears (l.99) that extrusion also involves thresholding (of the cell vertical displacement – l.99) in your case. How would this threshold influence the results?
There is a misunderstanding here as we have been very careful with this in the text. To reiterate, our method does not suffer from the arbitrariness of setting a threshold to defined extrusion events. The extrusion events emerge from the collective interactions of the cells and involve a cell detaching from the substrate. l.192 refers to what we use to detect (≠ define) an extrusion event as we postprocess the simulation results. This is a sufficient criteria as we have verified it through visualization of the simulations and a more computationally expensive method of the intersection of cell domains with the substrate.
References
[1] Lakshmi Balasubramaniam, Amin Doostmohammadi, Thuan Beng Saw, Gautham Hari Narayana Sankara Narayana, Romain Mueller, Tien Dang, Minnah Thomas, Shafali Gupta, Surabhi Sonam, Α S. Yap, Yusuke Toyama, Ren´eMarc M`ege, Julia M. Yeomans, and Benoˆıt Ladoux. Investigating the nature of active forces in tissues reveals how contractile cells can form extensile monolayers. Nature Materials, February 2021.
https://doi.org/10.7554/eLife.82435.sa2Article and author information
Author details
Funding
Multidisciplinary University Research Initiative (W911NF1910245)
 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. W911NF1910245. 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 ERCStarting 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.
Senior Editor
 Jonathan A Cooper, Fred Hutchinson Cancer Center, United States
Reviewing Editor
 Pierre Sens, Institut Curie, CNRS UMR168, France
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
© 2023, Monfared et al.
This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.
Metrics

 1,921
 Page views

 236
 Downloads

 2
 Citations
Article citation count generated by polling the highest count across the following sources: Crossref, PubMed Central, Scopus.
Download links
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)
Further reading

 Microbiology and Infectious Disease
 Physics of Living Systems
Microsporidia are eukaryotic, obligate intracellular parasites that infect a wide range of hosts, leading to health and economic burdens worldwide. Microsporidia use an unusual invasion organelle called the polar tube (PT), which is ejected from a dormant spore at ultrafast speeds, to infect host cells. The mechanics of PT ejection are impressive. Anncaliia algerae microsporidia spores (3–4 μm in size) shoot out a 100nmwide PT at a speed of 300 μm/s, creating a shear rate of 3000 s^{1}. The infectious cargo, which contains two nuclei, is shot through this narrow tube for a distance of ∼60–140 μm (Jaroenlak et al, 2020) and into the host cell. Considering the large hydraulic resistance in an extremely thin tube and the lowReynoldsnumber nature of the process, it is not known how microsporidia can achieve this ultrafast event. In this study, we use Serial BlockFace Scanning Electron Microscopy to capture 3dimensional snapshots of A. algerae spores in different states of the PT ejection process. Grounded in these data, we propose a theoretical framework starting with a systematic exploration of possible topological connectivity amongst organelles, and assess the energy requirements of the resulting models. We perform PT firing experiments in media of varying viscosity, and use the results to rank our proposed hypotheses based on their predicted energy requirement. We also present a possible mechanism for cargo translocation, and quantitatively compare our predictions to experimental observations. Our study provides a comprehensive biophysical analysis of the energy dissipation of microsporidian infection process and demonstrates the extreme limits of cellular hydraulics.

 Computational and Systems Biology
 Physics of Living Systems
The adaptive dynamics of evolving microbial populations takes place on a complex fitness landscape generated by epistatic interactions. The population generically consists of multiple competing strains, a phenomenon known as clonal interference. Microscopic epistasis and clonal interference are central aspects of evolution in microbes, but their combined effects on the functional form of the population’s mean fitness are poorly understood. Here, we develop a computational method that resolves the full microscopic complexity of a simulated evolving population subject to a standard serial dilution protocol. Through extensive numerical experimentation, we find that stronger microscopic epistasis gives rise to fitness trajectories with slower growth independent of the number of competing strains, which we quantify with powerlaw fits and understand mechanistically via a random walk model that neglects dynamical correlations between genes. We show that increasing the level of clonal interference leads to fitness trajectories with faster growth (in functional form) without microscopic epistasis, but leaves the rate of growth invariant when epistasis is sufficiently strong, indicating that the role of clonal interference depends intimately on the underlying fitness landscape. The simulation package for this work may be found at https://github.com/nmboffi/spin_glass_evodyn.