Two types of critical cell density for mechanical elimination of abnormal cell clusters from epithelial tissue

Recent technological advances in high-resolution imaging and artificial modulation of genetic functions at different times and regions have enabled direct observations of the formation and elimination of abnormal cell populations. A recent trend in cell competition research is the incorporation of cell mechanics. In different tissues and species, abnormal cells developing in epithelial tissues are mechanically eliminated by cell contraction via actomyosin accumulation at the interface between normal and abnormal cells. This mechanical cell elimination process has attracted attention as a potential universal defense mechanism. Here, we theoretically examined the conditions for mechanical elimination of growing abnormal cell populations. Simulations and mathematical analyses using a vertex dynamics model revealed two types of critical cell density associated with mechanical elimination of abnormal cell clusters. One is a subtype of homeostatic density, in which the frequencies of spontaneous mechanical cell elimination and proliferation are balanced, even if no explicit dependence of proliferation or apoptosis on the cell density is assumed. This density is related to the mechanical stability of a single cell. The other is density related to mechanical stability as a cell population under external pressure. Both density types are determined by tissue mechanical properties. In solid tissues, the former type is reached first as the intensity of interfacial contraction increases, and it functions as a critical density. On the other hand, the latter type becomes critical when tissues are highly fluid. The derived analytical solution explicitly reveals the dependence of critical contractile force and density on different parameters. We also found a negative correlation between the proliferation rate of abnormal cells and the likelihood of the abnormal cell population expanding by escaping elimination. This is counterintuitive because in the context of cell competition, fast-growing cell populations generally win. These findings provide new insight into, and interpretation of, the results from experimental studies.


Introduction
Epithelial tumors begin with the emergence of a small population of abnormal cells (e.g., via mutations in tumor suppressor genes), and eventually these abnormal cells infiltrate large areas of the body via growth and increase the risk of malignancy due to further mutations. The formation of a population from a single abnormal cell is a microscopic and random event, which has rendered research of this process difficult. Advances in technologies that enable artificial induction of gene abnormalities, as well as enhanced imaging techniques, have finally made it possible to observe these processes, which have become well studied in recent years, especially from the viewpoint of cell competition [1,2].
In the field of cell competition, the study of related chemical signal transduction systems was originally the primary focus [3][4][5][6][7][8][9][10][11]; however, cell mechanics have recently become relevant as well [12][13][14][15][16][17][18][19][20][21][22][23][24][25]. For instance, cultured Madin-Darby canine kidney II (MDCK) cells are hypersensitive to compaction after knockdown of the polarity gene scribble, and their interaction with wild-type cells causes compaction of knockdown cells and elevation of p53 levels in them, resulting in cell death [16,26]. When apoptosis is density dependent, the density threshold determines the tolerant density specific to each cell population, and mixing proliferating cell populations possessing different thresholds will result in survival of the population with the highest threshold. This phenomenon can be explained by homeostatic cell density/pressure, a physics concept proposed by Basan et al. [27]. In addition, misspecification of a subset of cells in the wing imaginal disc of Drosophila can trigger actomyosin accumulation at the borders between normal and abnormal cells, leading to contraction at the cell interface [14,28]. When the population size of abnormal cells is small, the cells are eliminated by cell death via this contraction, and when the population size is intermediate, a cyst is created [14]. The contraction of apoptotic cells with an internal ring-like structure composed of F-actin has also been observed [29,30]. Furthermore, in the context of organ development, increased density causes live-cell extrusion, without apoptosis, from the apical or basal side of an epithelial sheet [19,20]. In this manner, mechanical elimination processes such as mechanically driven apoptosis and live-cell extrusion have been observed in many different species and tissues and have attracted attention as a potential mechanism of homeostasis or protection against cancer.
Mechanical elimination of abnormal cells has been investigated both experimentally and theoretically. While experimental research focuses on the examination of actual specific cells or tissues [19], one of the advantages of a theoretical approach is that it allows comprehensive examination of elimination conditions under various situations and understanding of the overall perspective of the elimination process. In our previous study, motivated by the above findings, using numerical simulations of a vertex dynamics model and mathematical analyses, we derived the general conditions for mechanically eliminating a single abnormal cell via contraction of the cell edge. The minimum force required for elimination (termed critical contractility) was given as a function of different quantities such as the mechanical and geometrical cell parameters on a targeted abnormal cell [31].
In contrast, in more realistic situations, abnormal cells continue to proliferate to form a population. Thus, in the present study, we numerically and analytically examined the conditions necessary for elimination of growing abnormal cell populations via contraction at the interface between normal and abnormal cells. Our analysis revealed two types of critical density in abnormal cell populations for achieving mechanical elimination of the populations. The first is a subtype of the homeostatic density previously introduced by Basan et al. [27] and is related to the mechanical stability of a single cell. The other is related to the mechanical stability as a cell population. The values of both densities are determined by tissue physical properties (e.g., bulk and shear moduli), and the smaller of the two is the critical density above which abnormal cell clusters are completely eliminated. In addition, the analytical solution we derived also provides a good estimate of the critical interfacial contractility required for eliminating abnormal cell clusters with different physical properties.

A vertex dynamics model for epithelial tissue growth simulation
As stated in the Introduction, the mechanical elimination of abnormal cells has been reported under different situations. Here, we focus on the contraction force at the boundary between normal and abnormal cells as a mechanism for eliminating abnormal cell populations ( Fig 1A) and examined the dependence of the elimination success rate on different factors using a 2D vertex dynamics model (see Methods), which describes the apical dynamics of epithelial cells [32][33][34][35][36]. Although the elimination of abnormal cells from the epithelium is a 3D process, here we focused on the conditions for elimination at the apical surface because shrinkage of the apical surface is initially observed in many cases of mechanical cell elimination [29].
In a vertex dynamics model, the physical properties of a cell population or tissue are determined by the coefficients of line tension (Λ) and perimeter elasticity (Γ); more quantitatively, the bulk (K) and shear (G) moduli are given as functions of Λ and Γ (see Methods and also Fig  1B (left) for contours of K and G in (Λ-Γ) space). The values of K and G are negatively correlated in a wide range of regions in (Λ-Γ) space (Fig 1B, right). For smaller values of Λ and Γ, the tissue fluidity is higher [34,37,38]. As examples, according to the values of Λ and Γ estimated in previous works, MDCK cell monolayers are more fluid (K = 3.2 and G = 0.71) [29], whereas epithelial cells in the wing disc of Drosophila are less fluid, i.e., more solid (K = 2.01 and G = 1.27) [34]. In the following simulations and mathematical analysis, we examine the elimination conditions for tissues with different physical properties.
We considered two types of cells, normal and abnormal. Unless noted otherwise, both cell types were assumed to have the same physical parameters Λ and Γ, except at the interface between normal and abnormal cells. Following our previous study on the mechanisms of mechanical elimination of a single abnormal cell [31], edge contraction was reflected by Λ; the coefficient between normal and abnormal cells was set as μ-fold that between normal cells or between abnormal cells, where μ � 1 (Fig 1C).
It was assumed that the two cell types exhibit different cell cycle durations (τ N and τ A , respectively), with the cell cycle being shorter in an abnormal cell (τ A ) than a normal cell (τ N ) (see Methods regarding cell cycle modeling). Unless otherwise stated, the cell cycle of normal cells was set at infinity (i.e., the situation in which only abnormal cells proliferate was considered). Since abnormal cells generally continue to grow due to dysregulated proliferation, we also assumed that τ A is independent of abnormal cell density [39][40][41]. Note that we included some randomness in the cell cycle time to avoid synchronous cell division [37].
We considered two scenarios for the simulations. Scenario 1 is a somewhat artificial and simplified situation in which an abnormal cell cluster with N θ cells arises at a certain moment within a normal cell population, and immediately after its appearance, the interfacial contractile force comes into play; i.e., the coefficients of the interfacial edges were set to μΛ (Fig 1C,  top). The abnormal cell clusters were positioned to form a circle around the center of the normal tissue. After onset of the contractile force, the abnormal cells begin to divide. The reason  (right) The relationship between K and G, where the values calculated at 100 points in the gray region of (Λ-Γ) space are plotted. Red closed circles in the left and right panels (labeled as (i)-(ix)) show the parameter sets used for the simulations. (C) Two types of scenarios for the simulations. Scenario 1 is a artificial/simplified situation, in which an abnormal cell cluster appears at a certain moment within a normal cell population, and immediately after its appearance, the interfacial contraction μΛ comes into play (top). Scenario 2 is a more realistic situation, in which a single abnormal cell appears in a large population of normal cells and grows via cell division. Once the number of abnormal cells reach a certain size (N θ ), interfacial contraction μΛ occurs (bottom). See the text for details. (D) The relationship between the average stress state of the normal cells at the interface and the number of abnormal cells N (black lines). The gray regions show the standard deviation. σ 1 and σ 2 are the maximum (tensile, tangential to the interface) and minimum (compressive, perpendicular to the interface) principal stresses, respectively. In particular, σ 2 and σ 1 − σ 2 show a clear dependence on N. See [37] for the calculation of cell stress tensor.
https://doi.org/10.1371/journal.pcbi.1010178.g001 for evaluating the dynamics under this scenario is that, as we will see later, the densities (or average sizes) of normal and abnormal cells at the onset of the contractile force are the same, which simplifies the assumptions made when deriving the analytical solution for the elimination condition. Importantly, the basic logic behind the elimination dynamics examined in this simpler scenario, also holds in the more realistic Scenario 2 described next. Scenario 2 is a more realistic situation in which a single abnormal cell randomly arises within a large population of normal cells (e.g.,~1,000) and multiplies via cell division. Once the number of abnormal cells reaches N θ , contractile forces occur at the interface between the abnormal cells and surrounding normal cells (Fig 1C, bottom). This is under the assumption that surrounding normal cells can detect the appearance or existence of abnormal cells with some delay in the detection process, which should be correlated with the size of the abnormal cell population N θ . Thus, we dealt with N θ as a parameter and examined the system dynamics for different values of N θ . Furthermore, as shown in Fig 1D, the stress state of normal cells at the interface shows a clear correlation with the number of abnormal cells, specifically between the minimum principal stress (compressive stress) σ 2 and the number of abnormal cells N, and between the difference in principal stress σ 1 − σ 2 and N; note that σ 1 − σ 2 is strongly correlated with cell shape [37] (see also [37,42] for the calculation of cell stress tensor). Thus, Scenario 2 can also be interpreted as a scenario in which the surrounding normal cells recognize the existence of abnormal cell populations by the change in tissue stress, the magnitude of which is a function of N; in fact, myosin activity is higher at the cell edge under tissue stress during tissue development (e.g., [43][44][45]). Under Scenario 2, the abnormal cell density at the onset of contractility is higher than that in the surrounding normal tissue, and as will be seen, the mathematical analysis of the elimination condition is a bit more complicated. In the following sections, the simulation results under Scenario 1 are shown unless otherwise noted.

Two types of phase diagrams for elimination success/failure
We performed numerical simulations for cells with different values of K(Λ,Γ) and found two qualitatively different elimination mechanisms, depending on the tissue physical property K. In the case of tissues with lower fluidity (e.g., K = 2), three typical time courses for the number of abnormal cells were observed with the change in relative contractility μ (Fig 2A and 2B): (a) for smaller values of μ, elimination failed, and the abnormal cell population continued growing (Fig 2A, top), (b) for sufficiently large values of μ, complete elimination was successful (Fig 2A,  bottom), and (c) for intermediate values of μ, dynamic equilibrium or growth suspension was observed, in which the abnormal cell cluster maintained a certain size due to a balance between increasing cell number via proliferation and decreasing cell number via spontaneously occurring cell elimination due to mechanical instability (Fig 2A, middle). Regarding the growth suspension phase, finite time simulations cannot exactly determine if this equilibrium really exists in the limit of t ! 1. In this study, when the total number of abnormal cells was smaller than a certain threshold over 100 cell cycles of an abnormal cell (i.e., 100τ A ), the resultant time course was classified as this phase. As shown in Fig 2B, the frequency of this phase slowly decreases with the simulation time, but in some simulations, this dynamic equilibrium state was observed to last a very long time (e.g., more than 10,000τ A ). Biologically, suppressing the expansion of abnormal cell populations over a long period is critical for individual survival, so here we adopted the above classification with three phases, not two, i.e., elimination success/ failure. Furthermore, in the growth suspension phase, the size of the abnormal cell cluster remains almost constant for a long time and its size is dependent on the initial cluster size N θ at the onset of interfacial contraction (Fig 2C, top), but interestingly, the cell density of the cluster is nearly constant independent of N θ (denoted by ρ � ) (Fig 2C, bottom). Importantly,    suggesting that ρ � is the critical density associated with mechanical elimination of abnormal cell clusters. As shown later, ρ � is almost equal to the mechanical homeostatic cell density ρ 1 (see also Figs 3C and 4B). On the other hand, for the tissue with higher fluidity (F), 1=� r is an upwardly convex function at smaller μ values, without a plateau density, before jumping to zero at a higher value of μ. The cell density just before the jump, denoted as ρ �� , provides another critical density for mechanical elimination different from ρ 1 (see Fig 4C and the text for details). (G) The dependence of elimination success on interfacial contractility μΛ for different sets of mechanical parameters; the frequencies of elimination failure (red), growth suspension (black), and elimination success (green). Each curve was obtained using the same Hill function fitting explained above. It should be noted that the phase diagrams for elimination success/failure were similar for the sets of mechanical parameters yielding the same bulk modulus value K(Λ,Γ). Each vertical broken line in the phase diagram shows the critical contractility μ 2 Λ (black) or μ 1 Λ (red) obtained from analytical solutions (see the text and Fig 4D). the value of ρ � was independent of the simulation time (S1 Fig). Plotting the average area of the abnormal cell (1=� r) against μ for a given N θ (see Methods for the definitions of ρ � and � r) showed that 1=� r is a downwardly convex function for smaller values of μ. A range of μ values in which � r is constant (i.e., � r = ρ � ) is observed, and 1=� r jumps to zero at a larger μ (Fig 2D), suggesting that the density in the growth suspension phase ρ � is a critical density associated with mechanical elimination of abnormal cell clusters.
On the other hand, when the tissue has higher fluidity (e.g., K = 3.2), the profile of � r is qualitatively different from that when the tissue has lower fluidity; 1=� r is an upwardly convex function for smaller μ values and decreases continuously with increasing μ values, finally jumps to zero at a certain value of μ without entering the growth suspension phase (Fig 2E and 2F). We denote the cell density just before the jump as ρ �� , which appears to be another critical density associated with the mechanical elimination of abnormal cell clusters.
We examined the dynamics for different sets of mechanical parameters (Λ,Γ) and found that the phase diagrams for elimination success/failure were similar for the sets of mechanical parameters giving the same value of bulk modulus K(Λ,Γ) ( Fig 2G). This indicates that tissue physical properties determine which density, i.e., ρ � or ρ �� , is critical for the mechanical elimination of abnormal cell clusters. It should be noted that none of the simulation results changed qualitatively when a small amount of noise was added to the system (e.g., the one obeying the Ornstein-Uhlenbeck process [46] to the coefficient of line tension Λ).

Mechanical homeostatic cell density
In the growth suspension phase for tissues with lower fluidity, a balance between the increase in the number of cells due to cell division and the decrease due to mechanical cell elimination was observed, while keeping a unique density ρ � for a given set of (Λ,Γ). This balance is considered to be equivalent to the concept of homeostatic density, defined as the cell density with a steady state between cell division and apoptosis, introduced earlier by Basan et al. [27]. To date, signaling pathways for density-dependent regulation of cell division or apoptosis have been reported in some experimental studies [15][16][17]19,20,47], and the threshold of this regulation is considered to determine the homeostatic density. Importantly, in our model, we did not explicitly introduce such a density threshold (note that the T2 threshold is much smaller than 1/ρ � or 1/ρ �� and is not related to the focal density). In the growth suspension phase, as our previous study suggested, cell elimination is thought to occur spontaneously due to mechanical instability via saddle-node bifurcation when mechanical conditions allowing each cell to exist stably with a finite size are violated due to a density increase caused by cell division [31]. In fact, when cells continue to proliferate within a fixed wall (Fig 3A, left), they reach a plateau density specific to their mechanical parameters regardless of the size of the restricted space (Fig 3A, right). We call this mechanical homeostatic cell density ρ 1 (Λ,Γ), as a subtype of homeostatic density defined previously. We confirmed that the value of ρ 1 was in close agreement with the density in the growth suspension phase ρ � , i.e., ρ � � ρ 1 (Fig 2D).
At present, it is difficult to determine the exact value of ρ 1 from the physical properties, but it is possible to estimate it by statistical regression as follows. According to Lewis's law [48], the average size of a polygonal cell is highly correlated with its number of sides. In a proliferating cell population in the growth suspension phase or within a fixed wall, an increase in cell number due to proliferation is balanced by mechanical cell elimination, and triangular and/or quadrilateral cells, which generally have smaller areas, are more likely to be eliminated compared with polygonal cells with a greater number of sides. For example, in the case of (Λ,Γ) = (0.12, 0.04), which corresponds to the physical properties of epithelial cells in the Drosophila wing disc [34], the distribution of the number of sides in simulations indicates that the frequency of quadrilaterals in the growth suspension phase (or growth within a fixed wall) is significantly smaller than that for a growing tissue composed of a single cell type under the free boundary condition (Fig 3B). Previously, we showed analytically that the conditions under which each cell can exist in a mechanically stable manner are determined by the number of sides of the focal cell and the density of surrounding cells [31]; the analytical solution of the  upper limit of surrounding cell density below which a regular polygon with i edges can exist (denoted by ρ MCE(i) , where MCE stands for mechanical cell elimination) was obtained. The relationship ρ MCE(i) < ρ MCE(j) (i < j) holds, and the value of ρ MCE(i) depends on the physical parameters (Λ,Γ) (Fig 3C). In the case of (Λ,Γ) = (0.12, 0.04), since the rate of elimination of quadrilaterals and the rate of proliferation of all abnormal cells are balanced under the growth suspension phase, ρ 1 is expected to be bounded by the density at which pentagons, but not quadrilaterals, can exist stably, and the inequality ρ MCE(4) < ρ 1 < ρ MCE(5) actually holds for ρ 1 value obtained by the simulations (Fig 3C). For another parameter set, the elimination rate of triangles and the proliferation rate of all abnormal cells are balanced, and the inequality ρ MCE (3) < ρ 1 < ρ MCE(4) holds true (e.g., (Λ,Γ) = (0.06, 0.055)) ( Fig 3C). Interestingly, across different physical parameter sets, the values of ρ 1 and ρ MCE(i) (i = 3, 4, 5) are strongly correlated with each other (Fig 3D), and the following multiple linear regression gives a very good fit (Fig 3E): where c i are constants. Since ρ MCE(i) can be derived analytically for a given set of (Λ,Γ), ρ 1 can be predicted with high accuracy.
Regarding the critical density ρ �� in tissues with higher fluidity, ρ �� was unexpectedly different from ρ 1 (i.e., ρ �� 6 ¼ ρ 1 ) (Fig 2F), suggesting that there is another critical cell density other than mechanical homeostatic density associated with the mechanical elimination of abnormal cell clusters. In a later section, we will derive approximate analytical solutions for elimination conditions by mathematical analysis and will clarify the differences in the elimination mechanism depending on the physical properties of tissues and how the values of ρ � , ρ �� , and corresponding critical contractility are estimated.

Scaling of phase diagrams
Before addressing why the phase diagram differs depending on tissue physical properties, and how to determine the values of critical cell density and critical contractility necessary to achieve abnormal cell cluster elimination, we examined the dependence of phase diagram on the initial cluster size N θ . To achieve elimination success, the pressure difference between abnormal and normal tissues generated by interfacial tension derived from the relative contractility Δμ (= μ−1) must overcome the expansion of abnormal cells via proliferation. According to Laplace's law, in the case of the vertex dynamics model, this pressure difference between the two tissues (denoted by ΔP) is given by ΛΔμ/R, where R is the radius of the abnormal cell population (Fig 3F). Assuming that the square root of the number of abnormal cells at the onset of the contraction force, N θ , corresponds to R, the phase diagram for elimination success/failure is expected to scale with LDm= ffi ffi ffi ffi ffiffi N y p . Fig 3G and 3H show the superimposed diagrams obtained by simulations for different N θ values after rescaling the horizontal axis, demonstrating that the scaling relationship holds well. This scaling relationship was confirmed for different physical parameter sets (Fig 3G for a solid tissue and Fig 3H for a fluid tissue).

Mathematical analysis for the derivation of analytical solutions
To better understand the observations from numerical simulations, we attempted to derive the approximate analytical solutions for elimination conditions. Let us consider the following simplified situations (Fig 4A): (i) The abnormal cell cluster is composed of N cells and has a rotationally symmetric circular shape with density ρ. All abnormal cells have the same area a = 1/ρ, do not proliferate, and are in mechanical equilibrium under a given contractility μΛ.   There are two types of characteristic cell density that can be critical for the mechanical elimination of abnormal cell clusters; one is mechanical homeostatic density (ρ 1 , red horizontal broken line), and the other is related to mechanical stability as a population (ρ 2 , gray). For tissues with lower fluidity (B), the inequality ρ 2 > ρ 1 holds. Thus ρ 1 is reached first when μ increases and functions as a critical density for mechanical elimination, where the corresponding contractility is denoted by μ 1 Λ. On the other hand, for tissues with higher fluidity (C), ρ 2 < ρ 1 holds, and ρ 2 functions as the critical density, for which the corresponding contractility is μ 2 Λ. (D) Dependency of the cell area at equilibrium and its local stability on interfacial contractility for different sets of mechanical parameters under Scenario 1 (thick solid/dotted curves). The red and black vertical dotted lines represent μ 1 Λ and μ 2 Λ, respectively, and are the same as those in Fig  2G, showing that the derived analytical solution explains the simulation results well. In theory, when ρ 2 > ρ 1 , a growth (ii) At the interface, n abnormal cells are in contact with the surrounding normal tissue, and they are congruent hexagons with a perimeter of 6k ffi ffi ffi a p , where k is a constant. Each edge length between abnormal and normal cells (r) is k ffi ffi ffi a p . All normal cells are axial symmetric with respect to the radial axis through the center of the abnormal cell cluster, and they have the same area A N and perimeter L N . ffi ffi ffi a p 3 À 3q 1 ffi ffi ffi a p þ 2q 2 ¼ 0; ð2AÞ Since Eq (2) is the cubic equation of ffi ffi ffi a p ð¼ 1= ffi ffi ffi r p Þ, the area or density of abnormal cells can be solved analytically for a given contractility μΛ. When Eq (2) has finite positive real solutions, the abnormal cell cluster can exist; otherwise, the cluster disappears.
To calculate the roots of this equation, let us make the following additional assumptions about the quantities k, n, A N , and L N : (iv) The shape of each interfacial abnormal cell is close to a regular hexagon, k � ffi ffi ffi ffi ffi ffi ffi ffi ffi 2 ffi ffi ffi 3 p p =3.
(v) Regarding n, the relationship n � ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffiffi 2 ffi ffi ffi 3 p pN p holds under assumption (iv) and for large N (see Methods). Here, we adopted N θ as the value of N.
(vi) The area and perimeter in the ground state in Eq (6) (see Methods) are good candidates for those of surrounding normal tissues, A N and L N . Note that we focused only on the parameter region in which the ground state is given as the regular hexagonal packing [34,49]. Then, the relationship 49]. Fig 4B and 4C show the dependency of the cell area a (= 1/ρ) at the equilibrium of Eq (2) and its local stability on μΛ for a solid tissue with (Λ,Γ) = (0.12, 0.04) and a fluid tissue with (Λ,Γ) = (0.01, 0.025), respectively (see also Methods). With the change in μ, the number of equilibria varies. For a wide range of mechanical parameters, there is a bifurcation point at which positive real roots disappear (the black closed circle in Fig 4B); more rigorously, the bifurcation point appears when the tissue has a physical property for which q 1 > 0 is satisfied in the cubic equation, Eq (2). We denote the area/density and relative contractility at this suspension phase is expected. (E) The results of a similar analysis under Scenario 2; (bottom) phase diagram from the simulations and (top) analytical solutions. When the tissue fluidity is lower, the differences between ρ 1 and ρ 2 and between μ 1 Λ and μ 2 Λ are more marked. In (B) and (C), the parameter sets (Λ,Γ) = (0.12, 0.04) and (Λ,Γ) = (0.01, 0.025) were used, respectively. In (D), the nine parameter sets shown in Fig 1B were used. In (E), the three parameter sets (Λ,Γ) = (0.01, 0.025), (0.06, 0.035), and (0.12, 0.04) were used. All results were calculated for N θ = 100.
https://doi.org/10.1371/journal.pcbi.1010178.g004 bifurcation point as a 2 (or ρ 2 ) and μ 2 Λ, respectively, and those values can be given analytically as follows: Unlike mechanical homeostatic density (ρ 1 ), which is determined only by the cell mechanical parameters (Λ,Γ), ρ 2 depends on n (/ ffi ffi ffi ffi N p ); note that A N can be approximated as a function of (Λ,Γ) according to the above condition (vi). However, when the cluster size of abnormal cells is large enough, the effect of the term that includes n can be neglected, and ρ 2 is also determined solely by the cell mechanical parameters. Importantly, for tissues with higher fluidity (or larger bulk modulus K), this density with respect to mechanical stability as the cell population ρ 2 is in close agreement with ρ �� observed in the simulations under Scenario 1 and is lower than the mechanical homeostatic cell density ρ 1 ; i.e., ρ 2 � ρ �� < ρ 1 (Figs 2F and 4C). In contrast, for the tissues with lower fluidity, ρ 2 was larger than, and not consistent with, ρ � (� ρ 1 ) (Figs 2D and 4B).
These results can explain the qualitative difference in the phase diagram depending on the tissue physical properties observed in the simulations (Figs 2G and 4D). For a given physical property, ρ 1 (Λ,Γ) and ρ 2 (Λ,Γ) are determined. Since the physical property determines which of these is smaller, the one that is reached first differs when μ increases; when ρ 2 < ρ 1 , ρ 2 is the critical density for the elimination of abnormal cell clusters, and when ρ 2 > ρ 1 , mechanical homeostatic cell density ρ 1 becomes the critical density. Intuitively, when tissue fluidity is higher, the bulk modulus is larger and the shear modulus smaller, rendering it easier for individual cells to deform but more difficult for them to alter their volume. Therefore, when the abnormal tissue is compressed, the mechanical stability as a cell population becomes the criterion for elimination. For this reason, ρ 2 can be interpreted as the critical density. On the other hand, when tissue fluidity is lower, it is easier for individual cells to change their volume. Thus, when the interfacial contractility increases, the mechanical stability of individual cells becomes critical for elimination. For this reason, ρ 1 , which is strongly correlated with ρ MCE(i) (i = 3, 4, 5), as shown in Fig 3D, is the critical density.
The values of critical density and contractility necessary to achieve mechanical elimination of abnormal cell clusters differ depending on the scenario. Fig 4E shows the phase diagram obtained from simulations under Scenario 2 and the analytical solutions; it should be noted that under Scenario 2, condition (vi) was somewhat of an overestimation for the values of A N and L N probably because the abnormal cell density at the onset of the interfacial contraction was higher than that under Scenario 1. Thus we adopted the following modified condition (vi') for calculation of ρ 2 and μ 2 Λ: A N þ Gðn À pÞ nk ffi ffi ffi ffiffi a 0 p L N À a 0 À 6Gðn þ pÞ n À pL nk ffi ffi ffi ffiffi where a 0 is the mean area of abnormal cells just before the onset of interfacial contraction, and its value from the simulations was used. In addition, A N ¼ ð ffi ffi ffi 3 p =24ÞL 2 N was also assumed as in condition (vi).
As shown in Fig 4E, the critical density/contractility predicted by the analytical solution was in good agreement with the simulation. Expectedly, the contraction force required for elimination was smaller under Scenario 2 than under Scenario 1 for a given set of (Λ,Γ). One of the major differences in the results under Scenario 1 versus Scenario 2 is the gap between ρ 1 and ρ 2 and the disappearance of the bifurcation point in Scenario 2 for solid tissues (Fig 4E,  right). This can be understood by Eqs (2B) and (3A); under Scenario 2, A N is smaller and ρ 2 larger; especially, when A N is less than 6Γ(1+π/n), q 1 becomes negative, and the bifurcation point disappears.
On the other hand, the force balance equation can also explain the scaling of the phase diagram observed in the simulations. Assuming that n (/ ffi ffi ffi ffi N p ) is sufficiently large (specifically, n�π), for simplicity, Eq (2) can be rewritten as follows: Thus, the value of LDm= ffi ffi ffi ffi ffiffi N y p determines the success or failure of mechanical elimination of abnormal cells.

Effects of the relative proliferation rate and the difference in physical properties between normal and abnormal cells
Here, we show some other findings obtained from the simulations, in which Scenario 2 was used because it reflects more realistic situations in a biological sense.
(i) Abnormal cells are eliminated more easily when their relative proliferation rate is higher. In the above simulations, it was assumed that the normal cell cycle time was sufficiently long and negligible compared with that of the abnormal cell cycle time. Here, we explicitly considered the proliferation of normal cells and examined how the elimination success rate depends on the relative growth rate of abnormal cells to normal cells, r = (1/τ A )/(1/τ N ) = τ N /τ A . Note that we modulated the relative proliferation rate by changing the cell cycle duration of normal cells while leaving that of abnormal cells unchanged. Fig 5A shows the phase diagram for r = 1000, 100, and 10 for (Λ,Γ) = (0.12, 0.04). Interestingly, when surrounding normal cells proliferate, and the relative proliferation rate of abnormal cells decreases, the abnormal cells become difficult to eliminate. That is, slower proliferation of abnormal cells promotes their growth as a population. This might be counterintuitive but can be explained as follows. When a cell divides, its daughter cells have on average fewer sides compared with the mother cell [50,51]. For example, a hexagonal cell divides into two pentagonal daughter cells in general (Fig 5B). Simultaneously, the two cells adjacent to the divided cell exhibit an increase in their number of edges ( Fig 5B). As stated, the number of cell edges and the average cell size are positively correlated, referred to as Lewis's law [48]. Therefore, the size of the adjacent cells that obtained new edges tends to increase. When normal cells proliferate at a high frequency (i.e., when the relative growth rate of abnormal cells is low), division of normal cells at the interface increases the number of sides and the size of the adjacent abnormal cells; this leads to an increased likelihood of outward expansion of the abnormal cell population and also acts as a perturbation at the interface, rendering it difficult to maintain the growth suspension phase (Fig 5B). Conversely, when normal cells do not proliferate, or their proliferation is sufficiently slow, abnormal cells lose their edges and areas at interfaces via cell division and are pushed inward by normal cells, resulting in higher density. As mentioned earlier, the condition for stable existence of each polygonal cell depends on the density or the average size of its surrounding cells. In high-density cell populations, individual cells are sensitive to fluctuations in size due to the T1 process and/or division of surrounding cells, and the probability of them losing their stable equilibrium and disappearing via saddle-node bifurcation increases [31]. In a related study of the number of cell sides, Tsuboi et al. addressed the issue of which cell population fills the gap when cell death occurs at the interface between two proliferating cell populations [52]. Thus, for abnormal cells, rapid growth is not always an advantage.
(ii) Effects of differences in physical properties between normal and abnormal cells. Above, we analyzed cases in which abnormal and normal cells have the same physical properties. However, there are known cases in which the physical properties differ between mutant and normal cells [26,53,54]. Thus, we investigated the effects of different physical properties between normal and abnormal cells on the elimination of abnormal cell clusters. Fig 5C shows a phase diagram for the case when abnormal cells are more fluid (Λ A ,Γ A = 0.12, 0.01) and the surrounding normal cells more solid (Λ N ,Γ N = 0.12, 0.04). The critical contractility μΛ required for the elimination of abnormal cell clusters was somewhat smaller than (but not much different from) that in the case in which normal cells have the same physical properties as abnormal cells. This difference can be explained as follows: when normal cells are solid, the large shear modulus renders the T1 process less likely to occur at the interface, and the abnormal cells are trapped within the normal cells, increasing their density and facilitating their elimination. On the other hand, Fig 5D shows the phase diagram when abnormal cells are solid (Λ A ,Γ A = 0.12, 0.04) and normal cells are more fluid (Λ N ,Γ N = 0.12, 0.01). The major difference from the case in which both cell types are solid is that the growth suspension phase was not observed. This is because the T1 process occurs more readily when the surrounding normal cells are fluid due to a smaller shear modulus, thus acting as a perturbation to the topology PLOS COMPUTATIONAL BIOLOGY of interfacial abnormal cells and rendering it difficult to maintain the size of an abnormal cell population constant for longer. This also makes it more difficult to confine abnormal cells, leading to lower cell densities and an increase in the critical contractility μΛ required for elimination success, compared with the case in which normal cells have the same physical properties as normal cells (i.e., solid).

Discussion
In this study, motivated by experimental observations, we theoretically evaluated the conditions necessary for the elimination of a growing abnormal cell population by mechanical contraction at the interface between normal and abnormal cell populations. The main result of this study is the numerical and analytical demonstration that there are two types of critical densities for mechanical elimination of abnormal cell populations, and both depend on tissue physical properties. One of these critical cell densities, ρ 1 , is related to the homeostatic density previously introduced by Basan et al. [27]. However, we did not explicitly assume dependence of the abnormal cell proliferation rate on density, and the ability of the abnormal cell population to maintain a constant density ρ 1 even as it continues to proliferate is attributed to the spontaneous mechanical elimination of each cell. In this sense, ρ 1 is different from the original definition of homeostatic density, which explicitly assumes the dependence (i.e., threshold) of cell proliferation or death on density, and here we refer to ρ 1 as mechanical homeostatic cell density, the value of which is determined by cell mechanical parameters only. The other critical cell density, ρ 2 , is a quantity related to the mechanical stability as a population, and it satisfies the cubic equation as shown in the force balance equation, Eq (2). Note that Eq (2) has the same form as that of the force balance equation for a single cell with interfacial contractility between adjacent cells, which was analyzed in our previous study [31]. For a given physical property of cells, when ρ 2 > ρ 1 , ρ 1 is first reached as μ increases, providing the critical cell density required for mechanical elimination of the abnormal cell population, and vice versa. From the theoretical side, how the concepts of critical density developed here can be meaningfully generalized to a tensorial quantity (such as homeostatic stress tensor [55]) remains an important issue. We also found counterintuitive properties, such as a negative correlation between the proliferation rate of abnormal cells and the likelihood of their expansion by escaping elimination. As might be obvious, in some cases of cell competition, the fast-growing cell populations win [21]. However, our results indicate that when contractile forces act on the interface between two cell populations (i.e., μ > 1), rapid growth is not always an advantage. This would be an interesting point that can be verified experimentally. For example, in some experimental systems used to evaluate cell competition, it would be feasible to examine the effect on the elimination efficiency of abnormal cells when their relative growth rate is decreased by increasing the growth rate of normal cells, as we did in our simulation. In addition, it may be possible to demonstrate experimentally that the phase diagram for elimination success/failure can vary depending on the tissue physical properties, as shown by our simulations (Fig 2B and 2E). For example, the stiffness of human mesenchymal stem cells (hMSCs) depends on the substrate rigidity [56], and the stiffness of MDCK cells in the confluent state depends on the number of cells seeded [57]. Furthermore, it is also possible to induce local actomyosin activity and to control its intensity, in principle, using recent optogenetic techniques, which enable drawing a phase diagram by changing the contractility at the interface between normal and abnormal tissues.
The inability to describe the 3D effects of the cell elimination process is a limitation of our model. In our 2D model, the elimination of abnormal cells is indicated by their disappearance from the apical plane. However, as Bielmeier et al. reported experimentally, when the size of an abnormal cell population increases, it does not lose its apical surface but instead invaginates along the apico-basal plane and is eliminated after forming a cyst [14]. Extending our model to a 3D version might reveal conditions for mechanical cell elimination via apical exclusion and invagination followed by cyst formation, respectively.
In the current study, we assumed that the surrounding normal cells recognize the emergence of abnormal cells and that the contractile force at the interface between normal and abnormal cells increases once the abnormal cells reach a certain population size (in the case of Scenario 2). In real biological situations, some aspects of the signal transduction system involved in cell competition are known, although the recognition mechanisms are not wholly understood. Changes in the mechanical as well as chemical environment caused by the emergence of abnormal cells may provide a clue to their recognition. The possibility of mechanical recognition of abnormal cells would be an interesting issue both theoretically and experimentally to explore in the future.

A vertex dynamics model
In brief, in vertex dynamics models, each cell is represented as a polygon formed by linking several vertices, and the motion of each vertex obeys the following equations (non-dimensionalized version) [37]: where r i is the positional vector of the i-th vertex. The first term of the energy function U indicates an area constraint, where A α is the area of cell α. The second term represents the line tension and cell-cell adhesion at the edges between cell α and cell β, where l αβ is the edge length, and Λ αβ is the coefficient of line tension. The third term represents the perimeter elasticity, where L α is the perimeter of cell α, and Γ α is the coefficient of perimeter elasticity. Λ and Γ determine the physical properties of a cell population or tissue. Quantitatively, the bulk or shear modulus (denoted by K or G, respectively) is an index characterizing the physical properties of the tissue and, in the case of the above vertex dynamics model, is given as follows [34,49]: where l g (Λ,Γ) and A g (Λ,Γ) are the cell edge length and cell area at the ground state, respectively, with ground-state networks defined as the polygon configurations minimizing the energy function U for given physical parameters (Λ,Γ) within an infinite plane [34,49] (see Table 1 for the parameter values used in this study).
As a consequence of push-pull dynamics among cells in growing tissues, cellular rearrangement occurs. This rearrangement is driven by edge reconnections (T1 process) and occurs when the edge length is less than the T1 threshold θ T1 [34,37]. Cell elimination occurs simply by removing a cell with an area less than the T2 threshold θ T2 [31]. We mainly used y T1 = � l = y T2 = � A = 0.01, where � l is the average length of normal cells and � A is the average cell area of normal cells. We confirmed that the changes in the values of θ T1 and θ T2 had minor effects on our results; specifically, the relationships ρ � � ρ 1 and ρ �� � ρ 2 still hold, although the mechanical homeostatic cell density (ρ 1 ) changed slightly depending on θ T1 and θ T2 (please see below or the text for the definitions of ρ 1 , ρ 2 , ρ � , and ρ �� ) (S2 and S3 Figs).
In the simulations, the Euler method was used with the time step Δt = 0.0001, corresponding to 1/50000 of the cell cycle time, which is sufficiently small based on a previous study by Kursawe et al. [58]. As for boundary conditions, we tested both free and fixed boundary conditions and confirmed that the simulation results are qualitatively the same under both conditions; in this study, we show the simulation results under the free boundary condition.

Model for the cell division
We introduced a clock representing the cycle for each cell. The clock progresses linearly with time. When the timer reaches a threshold (τ A for abnormal cells and τ N for normal cells), the cell enters the mitotic phase of division. To prevent the synchronization of cell division, we included ± 20% randomness (uniformly distributed) into the cell cycle time. In most simulations, we used τ A = 5 and τ N = infinity. As shown in S4 Fig, the intensity of the interfacial contractility necessary for eliminating abnormal cell clusters depends somewhat on the cell cycle time, but the phase diagram for elimination success/failure converges to a certain profile with increasing cell cycle time. For tissues with lower fluidity, the relationship ρ � � ρ 1 holds well. For tissues with higher fluidity, the difference between ρ �� and ρ 2 becomes larger for larger τ A , but both ρ �� and ρ 2 still show close values. During the mitotic phase, the target cell area increases linearly with time and thus the first term of Eq (6B) is replaced by (A α − (1 + bt)) 2 / 2, where b is the growth rate of the target area, and we used b = 10. When the actual area of the focal cell (A α ) doubles, the cell divides with an axis through its center, and the cell cycle clock is reset to zero. A recent study of the cell cycle suggested that some epithelial cells divide by adding a constant volume between consecutive division events [59]. Thus, we also tested this model, called the adder model; specifically, a cell divides when the actual area of the cell (A α ) becomes A α + δA, where δA is a fixed increment, rather than doubling (i.e., 2 A α ). We confirmed that this change in the rule of cell division has a very minor impact quantitatively on our results (S5 Fig). Regarding cell division orientation, we performed simulations under the following two cases: (i) the orientation is determined randomly from the uniform distribution, and (ii) the division plane is perpendicular to the cell longitudinal axis (i.e., Hertwig's law [60]). Since there was no qualitative difference in the results between these cases, we present simulation results for case (i).

Definitions of characteristic densities
Since multiple characteristic densities appear in this paper, we summarize them here. For each simulation run, the density of an abnormal cell population at time t, ρ(t), was defined as rðtÞ � NðtÞ= P i A i ðtÞ, where A i (t) and N(t) are the area of the i-th abnormal cell and the total number of abnormal cells, respectively. For each of these time profiles of density (S1B Fig), � r, defined by � r � min t E t A ½rðtÞ�, was used as a quantity representing the time average of the density after the initial response to the interfacial contractile force, where E τA [] represents a moving average with a τ A time window. In the case of elimination success, 1=� r was defined as zero.
In the simulations for tissues with lower fluidity, a growth suspension phase was observed when the contractile force at the interface between abnormal and normal cells was moderate. In the growth suspension phase, the density was constant and depends on the cell mechanical parameters. This density observed in the simulation is called ρ � ; mathematically, 1/ρ � is defined as the population average of 1=� r for each sample path over all simulation runs showing the growth suspension phase, i.e., 1/ρ � � <1=� r> G.S. We showed that ρ � was in close agreement with mechanical homeostatic density (denoted as ρ 1 ), which is defined as the plateau density when cells continue to proliferate within a fixed wall. The value of ρ 1 depends on cell mechanical parameters but is independent of the size of the restricted space.
In the simulations for tissues with higher fluidity, another characteristic density related to mechanical elimination of abnormal cell clusters was observed, which we called ρ �� . ρ �� was defined as the density under the critical contractile force (i.e., at the jump point). We found this density to be in good agreement with the density ρ 2 , which is responsible for the mechanical stability of the abnormal cells as a population, where the value of ρ 2 was derived analytically.

Derivation of force balance equation Eq (2) and stability of equilibria
We considered a simplified model for deriving the approximate analytical solutions for elimination conditions (see the main text for the detailed settings and Fig 4A). The radial forces applied to the two interfacial vertices V 1 and V 2 (denoted by F 1 and F 2 , respectively) of a boundary abnormal cell are given as follows: F 1 ¼ F 1 e 1 ¼ ½ð1 À aÞr cos y 1 À 2mL sin y 1 À 12Gr sin y 1 À ð1 À A N Þr cosy 1 þ L þ 2GL N À 2GL N siny 1 �e 1 ; ð8AÞ F 2 ¼ F 2 e 2 ¼ ½ð1 À aÞr cos y 2 À L þ 12Gr sin y 2 À 12Gr À ð1 À A N Þr cos y 2 þ 2mL sin y 2 þ 2GL N sin y 2 �e 2 ; ð8BÞ where e 1 and e 2 are the unit vectors in the radial directions at V 1 and V 2 , respectively (see Fig  4A and the main text for the definitions of a, r, R, Λ, Γ, μ, A N , L N , θ 1 , and θ 2 ). When the force is balanced in each vertex within the abnormal cell population, the radial component (along the direction of e 1 ) of the net force (denoted by F) applied to the focal abnormal cell on the interfacial side with normal cells is calculated as Using the sum-to-product formulas, the following relationships hold: cos y 1 þ cos y 2 ¼ 2B cos sin y 1 À sin y 2 ¼ 2B sin where B = cos((θ 1 +θ 2 )/2). In addition, for sufficiently large n and small θ 1 +θ 2 , the approximations cos(π/(2n)) � 1, sin(π/(2n)) � π/(2n), sin 2 (π/(2n)) � 0, and B � 1 hold. Then, the righthand side of Eq (9) becomes difference between 1/ρ � and 1/ρ 1 for different cell cycle durations for tissues with lower fluidity. ρ 1 was very close to ρ � independently of cell cycle durations. (F) The difference between 1/ ρ �� and 1/ρ 2 for different cell cycle durations for tissues with higher fluidity. For larger τ A , the difference was larger (but not significant). This is thought to be because the density of abnormal cells becomes slightly smaller when there is sufficient time for area relaxation after division. In (