Active Deformable Cells Undergo Cell Shape Transition
Associated with Percolation of Topological Defects

A population of interacting self-propelled agents, often referred to as “active matter”, exhibits collective movement without external cues, as exemplified by migrating cell populations Marchetti et al. (2013) or flocks of birds Cavagna and Giardina (2014); Marchetti et al. (2013). Even with only excluded volume interactions, various phenomena characteristic to non-equilibrium active systems have been discovered, including giant density fluctuations Fily and Marchetti (2012), motility-induced phase separation Cates and Tailleur (2015), and active jamming Henkes et al. (2011). The shapes of the particles can alter the collective order; the motility-induced phase separation observed in spherical particles Cates and Tailleur (2015) vanishes in slightly non-spherical particles Großmann et al. (2020). A self-propelled rod exhibits the local orientational order Bär et al. (2020), and more complex collective patterns appear in more complex shapes Denk et al. (2016); Spellings et al. (2015); Wensink et al. (2014); Bär et al. (2020).

For the collective motions of cells, individual cells exhibit significant deformability, which was not considered in conventional models of active matter Fily and Marchetti (2012); Cates and Tailleur (2015); Henkes et al. (2011); Großmann et al. (2020); Bär et al. (2020); Denk et al. (2016); Spellings et al. (2015); Wensink et al. (2014); Bär et al. (2020). This deformability of the cells manifests in a densely packed situation and can drastically alter the tissue scale rheological properties, such as elasticity or viscosity, by deforming the cells themselves or facilitating cell rearrangements, and is closely linked to biological functions in embryogenesis, wound healing, and cancer invasion. The population of deformable active particles needs to be addressed to understand the underlying principle behind tissue organization and homeostasis.

Epithelial cells forming a confluent monolayer are well described by the cell-vertex model, in which the cells are densely packed with no gaps, and their shapes are approximated by polygons. The vertex model was recently used to explain the experimentally observed rigidity transition Angelini et al. (2011) from the a glassy phase, where the movements of cells are nearly frozen , to a fluidic phase, where frequent cell rearrangements result in tissue-scale collective motion Bi et al. (2015). Despite the constant density with a volume fraction of one, the vertex model undergoes a rigidity transition, where the deformability of cells triggers a cascade of cell rearrangements. The transition occurs when the target cell shape index $q$ defined by $\mbox{(cell perimeter)}/\sqrt{\mbox{(cell area)}}$ exceeds $q=3.81$. The self-propelled Voronoi model (SPV) that is almost equivalent to the vertex model, including the self-propelled force Bi et al. (2016) also showed the same rigidity transition at $\langle q\rangle=3.81$, where $\langle q\rangle$ is the average value of the shape index for each cell. These studies suggest that the rigidity transition is conditioned by the geometric parameter $q$: which is also supported experimentally Park et al. (2015).

On the other hand, it is yet to be well understood how cells other than epithelial cells, such as mesenchymal cells, for whom polygonal approximation is inappropriate, behave at high density. Although the rigidity transition in epithelial cells was initially studied in the context of epithelial-mesenchymal transition (EMT), experimental evidence suggests that EMT does not correspond to rigidity transition Mitchel et al. (2020). If so, it is important to clarify how the mesenchymal cell population acquires fluidity and mobility that often appear in biological processes and how they differ from the transition in epithelial cells Bi et al. (2015, 2016). To simulate cells with non-polygonal shapes, the phase-field model with self-propulsion was used to explore cells at high density Loewe et al. (2020), and it was shown that a rigidity transition similar to the vertex model occurs merely by considering the excluded volume interaction. A model of foam-like deformable cells with cell adhesion was also investigated Kim et al. (2021), where strong cell adhesion mediates the rigidity transition. For both studies, the rigidity transition occurred at a similar $\langle q\rangle$ value in the vertex model. However, these models can only handle a smaller number of cells than the vertex model, and it remains to be determined whether a mesenchyme-specific phase exists.

Here, we propose the Fourier contour cell model that allows us to handle up to $10^{4}$ deformable cells on a single CPU. The cells in the model have a round shape when isolated, whereas, in a high-density situation, the model can describe drastic cell deformation caused by the balance between the cell surface tension and excluded volume interactions. We demonstrate that a novel fluid phase specific to mesenchyme-like cells appears through excluded volume interactions and self-propulsion of each cell. The observed phase transitions were accompanied by the percolation of topological defects, providing a new perspective on mesenchymal cell dynamics that is experimentally verifiable.

Deformable cells interacting in two-dimensional space are considered here. We propose the Fourier contour cell model in which the cell contour is expressed by polar coordinates $R(\theta)$ (Fig. 1a). The cell contour $R(\theta)$ for $i$-th cell centered at the origin is expressed by a Fourier expansion up to $M$-th order, as follows:

$\displaystyle R^{i}(\theta)=R_{0}\sum_{n=0}^{M}\left[a_{n}^{i}\cos n(\theta-\theta^{i})+b_{n}^{i}\sin n(\theta-\theta^{i})\right]~{},$

(1)

where $\{a_{n}^{i}\}$ and $\{b_{n}^{i}\}$ are Fourier coefficients, $R_{0}$ is a constant parameter with radius dimensions, and $\theta^{i}$ denotes the orientation of the $i$-th cell, which corresponds to the self-propulsion direction of the cell. In this study, we adopt $M=6$. Considering the constraints that the cell area $\pi R_{0}^{2}$ is constant and that the center of the polar coordinates coincides with the cell centroid, $a_{0}^{i}$, $a_{1}^{i}$ and $b_{1}^{i}$ are determined as $a_{0}^{i}=\sqrt{1-\sum_{n=1}^{M}[(a_{n}^{i})^{2}+(b_{n}^{i})^{2}]/2}$ and $a_{1}^{i}=b_{1}^{i}=0$, respectively. Furthermore, we impose a constraint $\sum_{n=1}^{M}\sqrt{(a_{n}^{i})^{2}+(b_{n}^{i})^{2}}<\sqrt{2/3}$ to avoid self-crossing of the contour (see Supplemental Materials). The variables for describing $i$-th cell are composed of the centroid position $\mbox{\boldmath$r$}_{c}^{i}$, orientation of the self-propulsion $\theta^{i}$, and Fourier coefficients to describe shapes $\{a_{n}^{i}\}$ and $\{b_{n}^{i}\}$ ($n=2,3,\ldots,M$).

As the simplest interaction, only the excluded volume effect is considered. The interaction Hamiltonian $\mathcal{H}_{int}$ between $i$-th and $j$-th cells is given as an energetic penalty against the overlapped area. To compute this overlap, the field representation of the cell, originally used to represent the elliptical shape of the cell Großmann et al. (2020), was applied to describe deformable cells. The interior region of $i$th cell is described by $\phi^{i}(\mbox{\boldmath$r$})\sim 1$ (Fig. 1b) using $\phi^{i}(\mbox{\boldmath$r$})$ calculated as follows:

$\displaystyle\phi^{i}(\mbox{\boldmath$r$})=\frac{1}{2}+\frac{1}{2}\tanh\left(\frac{R^{i}\left(\theta(\mbox{\boldmath$r$},\mbox{\boldmath$r$}_{c}^{i})\right)-\Delta^{i}(\mbox{\boldmath$r$},\mbox{\boldmath$r$}_{c}^{i})}{\epsilon/2}\right)~{}~{}~{}$

(2)

where $\theta(\mbox{\boldmath$r$},\mbox{\boldmath$r$}_{c}^{i})$ is the angle between $\mbox{\boldmath$r$}-\mbox{\boldmath$r$}_{c}^{i}$ and $x$-axis, and $\Delta^{i}(\mbox{\boldmath$r$},\mbox{\boldmath$r$}_{c}^{i})=|\mbox{\boldmath$r$}-\mbox{\boldmath$r$}_{c}^{i}|$. With small $\epsilon$, function $\phi^{i}(\mbox{\boldmath$r$})$ sharply rises to $1$ inside the cell (i.e., $\Delta^{i}<R^{i}$) and drops to $0$ outside the cell (Fig. 1b). Interface width between $\phi^{i}=0$ and $\phi^{i}=1$ domains is controlled by $\epsilon$, which is considered infinitesimal below. The interaction Hamiltonian is then given by $\mathcal{H}_{int}=\sum_{i<j}\int d\mbox{\boldmath$r$}\phi^{i}\phi^{j}$ where we set the coefficient to unity. We also incorporate the energy for penalizing interface length associated with the membrane tension and obtain the total Hamiltonian $\mathcal{H}=\mathcal{H}_{int}+\eta\int_{-\pi}^{\pi}\sqrt{(R^{i})^{2}+(R^{i}{}^{\prime})^{2}}d\theta$, where $\eta$ is the tension parameter, and $R^{i}{}^{\prime}$ denotes the derivative of $R^{i}$ with respect to $\theta$. The second term does not depend on $\mbox{\boldmath$r$}_{c}$ and $\theta^{i}$.

The time-evolution equation for $\mbox{\boldmath$r$}^{i}_{c}$ is derived from the functional derivative of $\mathcal{H}$ with the self-propulsion term $\mbox{\boldmath$v$}^{i}=(v_{0}\cos\theta^{i},v_{0}\sin\theta^{i})$, which is often introduced for representing active motion of cells Fily and Marchetti (2012); Cates and Tailleur (2015); Henkes et al. (2011); Großmann et al. (2020); Bär et al. (2020); Denk et al. (2016); Spellings et al. (2015); Wensink et al. (2014); Bär et al. (2020):

$\displaystyle\dot{\mbox{\boldmath$r$}}_{c}^{(i)}$

$\displaystyle=$

$\displaystyle\mbox{\boldmath$v$}^{i}+\frac{\delta\mathcal{H}}{\delta\phi^{i}}\nabla\phi^{i}=\mbox{\boldmath$v$}^{i}+\sum_{j\neq i}\int d\mbox{\boldmath$r$}\phi^{j}|\nabla\phi^{i}|\frac{\nabla\phi^{i}}{|\nabla\phi^{i}|}$

(3)

By taking the sharp interface limit $\epsilon\to 0$, $|\nabla\phi^{i}|$ becomes the surface delta function $\delta_{s}(\mbox{\boldmath$r$}-R^{i})$ (see supplemental text), and the second term is replaced by an integral along the cell contour $\sum_{j\neq i}\oint ds\phi^{j}\mbox{\boldmath$n$}^{i}$ where $\mbox{\boldmath$n$}^{i}$ denotes the normal vector of the contour. By using $\mbox{\boldmath$e$}_{\theta}=(\cos\theta,\sin\theta)$ and $\mbox{\boldmath$e$}_{\perp}=(-\sin\theta,\cos\theta)$, $\dot{\mbox{\boldmath$r$}}_{c}^{(i)}$ is given by

$\displaystyle\dot{\mbox{\boldmath$r$}}_{c}^{i}$

$\displaystyle=$

$\displaystyle\mbox{\boldmath$v$}^{i}+\sum_{j\neq i}\!\int_{0}^{2\pi}\!\!\!\!d\theta\phi^{j}\!\left[R^{i}{}^{\prime}(\theta)\mbox{\boldmath$e$}_{\perp}\!\!-\!R^{i}(\theta)\mbox{\boldmath$e$}_{\theta}\right]~{}~{}~{}$

(4)

At the last integral, $\phi^{j}=1$ when the position $\mbox{\boldmath$r$}_{c}^{i}+R^{i}\mbox{\boldmath$e$}_{\theta}$ on the contour in the $\theta$-direction of the $i$-th cell is occupied by the $j$th cell, and $\phi^{j}=0$ otherwise. The time evolutions of $\theta^{i}$, $\{a_{n}\}$ and $\{b_{n}\}$ ($n=2,3,\ldots,M$) are obtained in the same manner as $\dot{\theta^{i}}\propto-\delta\mathcal{H}/\delta\theta+\sqrt{2D_{r}}\xi^{i}$, $\dot{a}_{n}\propto-\delta\mathcal{H}/\delta a_{n}$ and $\dot{b}_{n}\propto-\delta\mathcal{H}/\delta b_{n}$. Note that the tension term in the Hamiltonian affects only the equations for $\{a_{n}\}$ and $\{b_{n}\}$. For simplicity, we incorporate a noise term $\sqrt{2D_{r}}\xi^{i}$, where $\xi^{i}$ represents a normalized white Gaussian noise, only into the equation of $\dot{\theta}^{i}$. The numerical integration along the cell contour $\theta$ was performed by discretizing $0\leq\theta\leq 2\pi$ by 40 points, whereas the integration with respect to time $t$ was performed using the Euler-Maruyama method with $\Delta t=5\times 10^{-3}$. The dynamics in the centroid $\mbox{\boldmath$r$}^{i}_{c}$, orientation $\theta^{i}$ and cell shape $\{a_{n}^{i},b_{n}^{i}\}$ are described by a set of ordinary differential equations (ODEs): where the two-dimensional numerical integration in Eq. (3) is reduced to a one-dimensional integration, which significantly reduces the computation time. Figure 1(c) shows a representative simulation snapshot of $N=10^{4}$ active deformable particles (see also VIDEO1 in Supplemental Material), which may not be feasible using the phase-field method Loewe et al. (2020); Löber et al. (2015); Nonomura (2012).

Using the proposed model, we study active deformable particles densely packed with volume fraction $\rho=0.95$ and interact with each other solely through excluded volume interactions. The simulation box size was set as $L_{x}:L_{y}=2:\sqrt{3}$ and the particles were initially positioned on a perfect hexagonal lattice with random $\theta^{i}$. Although up to $N=10^{4}$ particles were handled, calculations were performed mainly with $N=1024$ when searching the parameter space.

Representative snapshots are shown in Figs. 2(a)-(c). At a small self-propulsion $v_{0}$ and a large tension $\eta$, the system exhibits a solid phase (Fig. 2a) where the particles form a hexagonal lattice and hardly show rearrangements (particle trajectories are shown in the inset of Fig. 2a; see also VIDEO2 in Supplemental Material). Accordingly, the mean square displacement (MSD) measured from the trajectories of the cell centroids exhibits a saturation curve (Fig. S1; $\eta$ = 0.08, 0.1). The particles in this phase were nearly circular without a large deformation (Fig. 2a). At large $v_{0}$ and large $\eta$, the system is in the fluid phase (Fig. 2b), where the particle positions are frequently rearranged (Fig. 2b inset; see also VIDEO3 in Supplemental Material) and the MSD curve exhibited linear growth, indicating diffusive dynamics of cells (Fig. S1; $\eta$ = 0.02 - 0.06). The particle shapes remained almost circular (Fig. 2b; the degree of deformation is denoted by red color). At a small $\eta$ and large $v_{0}$, where the particles are easily deformed, the third phase, which we call the soft-fluid phase, emerges (Fig. 2c; VIDEO4 in Supplemental Material). In this phase, largely deformed and relatively circular particles coexist, and they exhibit fluidic characteristics, such as frequent rearrangements (Fig. 2c, inset) and a linearly increasing MSD curve (Fig. S1; $\eta$ = 0.01).

Figure 2(d) illustrates the phase diagrams for $\eta$ and $v_{0}$. Color represents the effective diffusion coefficient $D_{\rm eff}=D_{\rm s}/D_{0}$, where $D_{0}=v_{0}^{2}/2D_{\rm r}$ is the diffusion coefficient for an isolated cell. and $\displaystyle{D_{\rm s}=\lim_{t\to\infty}\left(\langle(\Delta r)^{2}\rangle-\langle\Delta r\rangle^{2}\right)/4t}$ were estimated from the trajectories of the cell centroids. The solid phase is determined by the criterion $D_{\rm eff}<10^{-3}$, which is consistent with the MSD curves showing saturation (Fig. S1 and S2). The solid-fluid phase boundary in Fig. 2(d) indicates that the rigidity-fluidity transition can occur not only by changing the cell propulsion $v_{0}$ but also by changing the cell deformability $\eta$. For $v_{0}=0.03$, $D_{\rm eff}$ varies abruptly at the boundary, and the MSD curves (Fig. S1) change from a diffusion line ($\eta\leq 0.06$) to a saturating curve ($\eta\geq 0.08$). The phase boundary is also associated with a decrease in the hexatic order parameter $|\Psi_{6}|$ (Fig. S3) defined by $\Psi_{6}(\mbox{\boldmath$r$})=\langle\sum^{n_{i}}_{j=1}e^{i6\theta_{ij}}/n_{i}\rangle$, where $\theta_{ij}$ is the angle of the link connecting $i$-th cell’s centroid to the adjacent $j$-th cell’s centroid, $n_{i}$ is the number of adjacent pairs, and $\langle\rangle$ denotes the average over all cells.

The phase boundary between the fluid and soft-fluid is not straightforwardly determined. As these two fluid phases have distinct cell shape characteristics, the distribution of the shape index $q$ is measured (Fig. 2e-g) to characterize the phases. For a large $\eta$ (high tension), the shape distribution displays a sharp, unique peak at $q<3.81$ regardless of whether it was in the solid phase (small $v_{0}$; Fig. 2e, dashed lines) or fluid phase (large $v_{0}$; Fig. 2e, thick lines). On the other hand, for a small $\eta$ (low tension), a sharply peaked distribution at $q<3.81$ in the solid phase (small $v_{0}$; Fig. 2e, dashed lines) turns into a long-tailed distribution as $v_{0}$ increases and then into a bimodal distribution (Fig. 2e, thick lines). Changes in the distribution from unimodal to bimodal also occur without crossing the solid phase. For a fixed value of $v_{0}=0.06$, we observe the appearance of a bimodal distribution with a second peak at $q>3.81$ from the unimodal distribution with a peak at $q<3.81$ when $\eta$ decreased (Fig. 2g). This bimodal distribution clearly illustrated the coexistence of largely deformed and relatively circular particles. Accordingly, the mean $\langle q\rangle$ increases around the region with a small $\eta$ and high $v$ (Fig. S3). Based on this change in the distribution and the previously reported transition point in the vertex model Bi et al. (2015, 2016), we tentatively defined the phase boundary as the mean shape index $\langle q\rangle=3.81$ (Fig. 2d, dashed red line). This heuristic definition will be justified below.

To further examine the boundary between the fluid and soft-fluid phases, we examined behaviors of topological defects characterized by miscoordinated particles with respect to a perfect hexagonal lattice. After performing Voronoi tessellation using particle centroids, $\{\mbox{\boldmath$r$}_{c}^{i}\}$, defect particles are defined as Voronoi cells with other than six neighboring cells. Top panels in Fig. 3(a) shows the defect particles indicated by the colors for $v_{0}=0.06$ at $\eta=0.01,0.03$, and $0.1$. As $\eta$ decreases from $0.1$ to $0.01$, the number of defect particles increases, and they are connected to form a large cluster (Fig. 3(a)). Figure 3(b) demonstrates a clear hallmark of the percolation transition the largest cluster size of the connected defect particles exhibited a sharp transition between $\eta=0.02$ and $0.06$, and the transition becomes sharper with increasing system size $N$ with a fixed density. Consistently, the size distribution of the connected clusters of the defect particles exhibited a power law near the critical point $\eta=0.03$ for $v_{0}=0.06$ (Fig. 3c). The power exponent is close to the critical exponent for 2D percolation, $187/91$ Stauffer and Aharony (2018). The cluster size distributions for the other values of $\eta$ are shown in Fig. S5. In contrast to the largest cluster size, the defect ratio $\rho_{d}$, defined as the average number of defect particles/$N^{2}$, does not exhibit a significant change, such as a transition. During a decrease in $\eta$, $\rho_{d}$ increases smoothly and its curve does not change with increasing $N$ (Fig. S6). The scatterplot of the defect ratio vs. the size of the largest cluster for all examined parameters collapses to a master curve (Fig. 3d; see also Fig. S7 for the unscaled version) by applying the percolation threshold $\rho_{d\rm c}=0.5$ for two-dimensional site percolation in the triangular lattice and its scaling exponents $v=4/3$ and $\beta=5/36$. This result clearly demonstrates the existence of percolation transition at $\rho_{d\rm c}$ = 0.5. Interestingly, the parameter region with $\rho_{d}>0.5$ in the phase diagram coincided with the region with $\langle q\rangle>3.81$ (Fig. 3e; see also Fig. S4), indicating that the percolation transition separates the soft-fluid phase from the fluid phase.

One possible interpretation is that the fluid/soft-fluid transition found in the proposed model corresponds to the hexatic/fluid phase transition in Kosterlitz-Thouless- Halperin-Nelson-Young (KTHNY) theory. According to the KTHNY theory Strandburg (1988); Bernard and Krauth (2011), two-dimensional crystals display two-step melting. At the first transition, the quasi-long-range order of a crystal is destroyed by the dissociation of the thermally excited pairs of dislocations, which is a linear crystallographic defect, while the six-fold orientational order of the lattice is maintained. This phase is referred to as the hexatic phase. Subsequently, the hexatic phase undergoes the second phase transition to the isotropic fluid phase by losing the orientational order mediated by the dissociation of disclinations Strandburg (1988); Bernard and Krauth (2011), a defect in which the rotational symmetry is broken. This two-step melting was confirmed in passive Bernard and Krauth (2011) and active discs Klamser et al. (2018), as well as in the Voronoi model Li and Ciamarra (2018); Pasupalak et al. (2020). The question that arises here is whether the fluid and soft-fluid phases in the proposed model correspond to either hexatic or isotropic fluid phases in the KTHNY theory. We addressed this issue by closely examining the transition from the solid to the fluid phase. We measured the pair translational correlation function $G_{r}(\Delta x,0)$ Bernard and Krauth (2011), which is calculated using the one-dimensional cut of the two-dimensional histogram $G_{r}(\Delta\mbox{\boldmath$r$})$ of the pair distance $\Delta\mbox{\boldmath$r$}=\mbox{\boldmath$r$}_{j}-\mbox{\boldmath$r$}_{i}$ of an $i$-$j$ pair of particle centroids, and orientational correlation function $G_{6}(|\Delta\mbox{\boldmath$r$}|)$ defined by $G_{6}(|\Delta\mbox{\boldmath$r$}|)=\langle\Psi_{6}(\mbox{\boldmath$r$}_{i})\cdot\Psi_{6}^{*}(\mbox{\boldmath$r$}_{j})\rangle$.

We calculate $G_{r}(\Delta x,0)$ and $G_{6}(|\Delta\mbox{\boldmath$r$}|)$ by fixing $\eta=0.1$ and increasing $v_{0}$ from $0.01$ to $0.05$. As shown in Fig. 4(a), from $v_{0}=0.01$ to $0.03$, the translational order $G_{r}(\Delta x,0)$ decays algebraically with exponent $-1/3$, corresponding to the solid phase in KTHNY theory. $G_{r}(\Delta x,0)$ decreases exponentially for $v_{0}$ from $0.04$ to $0.05$. The orientational order $G_{6}(|\Delta\mbox{\boldmath$r$}|)$ shown in Fig. 4(b) decays more slowly than $G_{6}(|\Delta\mbox{\boldmath$r$}|)\propto|\Delta\mbox{\boldmath$r$}|^{-1/4}$ for $v_{0}$ from $0.01$ to $0.04$, corresponding to solid or hexatic phases, and decreases faster than $|\Delta\mbox{\boldmath$r$}|^{-1/4}$ for $v=0.04$ and $0.05$. These results indicate that the solid and fluid phases in the proposed model correspond to the solid and isotropic fluid phases in KTHNY theory, and the hexatic phase exists between them for $v_{0}=0.04$ and $\eta=0.1$ (Figs. 4a and b). During the transition from the fluid phase to the soft-fluid phase, $G_{r}(\Delta x,0)$ and $G_{6}(|\Delta\mbox{\boldmath$r$}|)$ decayed exponentially (Figs. 4 (c) and (d); $v_{0}=0.06$ and $\eta$ is changed), indicating that both translational and orientational orders are already broken. From these results, we conclude that both the fluid and soft-fluid phases in the proposed model are classified into the isotropic fluid phase in terms of KTHNY theory, and the soft-fluid phase is a novel phase that can be described by the percolation of topological defects but not by breaking orientational or translational orders. This scenario is summarized in the phase diagram shown in Fig. 4(e).

The present study proposes a numerical model of deformable cells based on the Fourier expansion of the cell contour, which provides a computationally efficient method for describing irregular cell shape deformations similar to mesenchymal and ameboid cells. Up to $10^{4}$ cells were handled using a single CPU, which may not be feasible in the phase-field method Loewe et al. (2020); Löber et al. (2015); Nonomura (2012) or spring-beads models Kim et al. (2021). Similar models based on the expansion of cell contours have also been proposed Ohta (2017); Yamanaka and Ohta (2014); Menzel and Ohta (2012); Itino et al. (2011); however, these models had difficulty simulating a tightly packed situation and torque effects, and thus phenomena such as glassy dynamics and rigidity transition do not occur Itino et al. (2011). Our approach is based on the field representation of cell shape $\phi^{i}$, which enables simple and efficient computation of the excluded volume effect at high density and allows us to incorporate the torque effect appropriately, which can alter collective behaviors Hiraiwa et al. (2022); Großmann et al. (2020).

The proposed model demonstrated that density-independent rigidity transition occurs by changing deformability $\eta$ solely from the excluded volume interaction and self-propulsion. Although such a transition was reported in a previous study Loewe et al. (2020), we further elucidated the emergence of two types of fluid phases: the soft-fluid phase, where the cell shapes drastically deform, and the fluid phase, where individual cells remain almost circular. The soft-fluid phase is characterized by $\langle q\rangle>3.81$ and is thus similar to the phase that appears in the cell-vertex model Bi et al. (2015, 2016), whereas the fluid phase seems to correspond to the fluidic phase in the active Brownian particle model (ABP) since the cell shapes remain circular. These two phases stem from the nature of the proposed model, in which the model displays ABP behavior in the high tension limit, while deformability plays an essential role in low tension. Interestingly, the solid-fluid transition point $\langle q\rangle=3.81$ in the cell-vertex model corresponds to the liquid-liquid transition point in the proposed model. This difference would be explained by the fact that a fluid phase with $\langle q\rangle<3.81$ cannot appear in the cell-vertex model.

The fluid/soft-fluid transition is well-characterized by the percolation of topological defects. Recently, such defects percolation has been studied in equilibrium disks and ABP systems by changing the packing fraction, and it was found that percolation occurs around the hexatic-fluid phase-transition boundary  Digregorio et al. (2022). This difference between the previous study and our study may be due to the density-dependent/independent nature of the transition or the absence/presence of deformability of cells. Our finding of percolation in the liquid phase also differs from previously reported percolations related to the rigidity transition, such as density-dependent percolation in embryo genesis Petridou et al. (2021) or stress field percolation in the epithelial monolayer Monfared et al. (2022). In our case, the analysis of the translational and orientational order parameters revealed that the fluid/soft-fluid transition does not correspond to the hexatic/isotropic fluid transition in KTHNY theory, and both phases have the characteristics of an isotropic fluid. Moreover, the hexatic phase exists around the solid/fluid transition point (see Figs. 4a and b; $\eta=0.04$ and Fig. 4e), indicating that three fluid phases can potentially appear in the active deformable cell population.

The fluidic collective motion of actual cells can be characterized using the present theoretical framework. The packing of cell populations with non-polygonal shapes appears in embryo genesis Kim et al. (2021) and in vitro experiments on cultured cells with suppressed cell-cell adhesion Balasubramaniam et al. (2021), and they potentially display fluidic collective motion by cell-softening. The percolation of topological defects can be detected under these experimental conditions and provides insights into the novel fluid/fluid transition in the cell population.

The proposed model also has the flexibility to incorporate various ingredients, such as cell size differences, cell adhesion, and cell division, which need to be addressed elsewhere. The interaction Hamiltonian is based on the overlapped area between two cells, but its extension to a harder core potential is also possible. Extension to a three-dimensional model using spherical harmonic expansions instead of Fourier series expansions can be considered. Thus, the proposed model, which enables a simple and efficient computation of deformable cells, can be a pivotal tool for studying cell populations.

We would like to acknowledge the helpful discussions with Michio Tateno, Atsushi Ikeda, Hideyuki Mizuno, Norihiro Oyama, Kyohei Takae, and Tetsuya Hiraiwa. This study was supported in part by the Japan Society for the Promotion of Science (JSPS) KAKENHI(21K03496 and 21H05793 to NS). Japan Science and Technology Agency (JST) CREST JPMJCR1923 ( NS and SI). and the Cooperative Study Program of Exploratory Research Center on Life and Living Systems (ExCELLS; program No. 20-102, 21-102 to NS)