I. Introduction

Tissue fluidity is modulated in diverse biological contexts such as embryonic development [1, 2], wound healing [3, 4], and cancer progression [5, 6]. A solid-to-fluid transition can sculpt tissues during tail elongation in developing zebrafish embryos [1, 2, 7]. This phase transition also coincides with the cell-state switching and changes in gene expressions, called Epithelial-to-Mesenchymal transition (EMT) [1, 8]. EMT has been proposed to mimic an unjamming transition in tissues [9–11]. Fluid-solid transitions were also observed in cell monolayers [12, 13], where transitions occur either with increasing cell density (due to proliferation) [13] or with increasing cellular adhesion [12]. Tissue monolayers further undergo solid-to-fluid transitions by applying external compressive stress [10, 11]. How tissue-level transitions emerge from governing cellular parameters is an important question yet to be fully understood [9, 14, 15].

Two classes of computational frameworks mainly exist to explain tissue-level transitions: Self-propelled particle (SPP) models and models that treat cells as extended objects. SPP models describe cells as soft discs [16–18] and have been used to study the density (or packing fraction) dependent [19–23] and motility-driven [12, 19, 20, 23] jamming transitions. However, SPP models can not capture the cell shape changes in a confluent tissue where there are almost no gaps between cells (i.e., packing fraction close to unity). Meanwhile, the Vertex model and its variants [24–31], Self-propelled Voronoi model [32–34], Deformable Particle model [35], Cellular Potts model [36–38], and Phase Field model [39] describe cells as extended objects accounting for the cell geometry. Cell dynamics in all these models depend on energy functions, unlike SPP models that describe the force-based motion of cell centers. For instance, 2D Vertex and Self-propelled Voronoi (SPV) models define a shape-energy function that depends on each cell’s target area (a₀) and perimeter (p₀). These models predict a density-independent solid-to-fluid transition with an increasing target shape index, so, defined as a dimensionless parameter (s₀ = p₀/). However, the increase in shape index was interpreted as an increase in cell-cell adhesion [26, 32]. This counterintuitively implies that increasing adhesion drives tissue fluidization.

In contrast, increasing cell-cell adhesion induces solidification in MDCK cell monolayers [12, 40] and organotypic cultures [41]. Also, the downregulation of adhesion molecules in the elongating zebrafish tail coincides with EMT and tissue fluidization [1, 8]. Thus, the role of intercellular adhesion on the tissue scale transitions is still debatable. We address this by introducing a 2D active force-based model of tissue monolayers, where each cell made with beads and springs represents a soft deformable object.

We combine the particle-like nature and force-based dynamics of SPP models with the spatial extension of cells as in Vertex and SPV models, eliminating the deficiencies of both types of models. For instance, our model displays spontaneous neighbor exchanges, while Vertex and SPV models implement these events in an ad hoc manner, using a threshold value on cellular edge lengths. Additionally, the interface of adjoining cells in Vertex and SPV models is represented by a shared (common) edge, making it impossible for adjoining cells to move relative to one another, thereby dictating a geometric constraint [42]. Our model removes this unphysical constraint, allowing for spontaneous fluctuations in cell shapes. Notably, our model significantly differs from topological energy-based models, like the Vertex and Voronoi models, where a single cell loses its physical meaning in isolation (see Section III in SI).

In our model, increasing cell-cell adhesion leads to liquid-to-solid transition, contrary to the predictions of Vertex and Voronoi models. Our observation is consistent with in vitro experiments on tissues [12, 41] and in vivo experiments on morphogenesis during embryonic development [1, 8]. Near the liquid-solid phase boundary, we also observe glassy dynamics with many signatures of dynamic heterogeneity. In the glassy regime, we collapse distinct exponential tails of displacement distributions onto a master curve, revealing a universal ‘diffusion length’ that explains the origin of dynamic heterogeneity. Our framework provides precise control over individual cell properties and cell-cell interaction to gain mechanistic insights into emerging tissue behavior, with potential applications in various developmental and physiological contexts.

II. Model

We adopted a force-based passive mechanical model of 2D soft grains (or cells) [43–45] instead of commonly used frameworks of energy-based Vertex or Voronoi models. A single cell is assumed to be a closed loop of beads connected by springs with stiffness K_s and natural length l₀ (Fig. 1A(top), Fig. S1). Each bead encounters tangential tension from adjacent springs and internal pressure, P. These intracellular forces capture the roles of actomyosin cortex and cytoplasmic pressure, respectively. Thus, the total intracellular force experienced by the i-th bead in a cell is:

Here l_i is the bond length between i-th and (i − 1)- th beads; while and denote the tangential and the outward normal unit vectors corresponding to the i-th bead, respectively (Fig. S1). Moreover, beads of distinct cells experience short-ranged intercellular forces consisting of two parts: a spring-like attraction, describing the cell-cell adhesion, and a spring-like repulsion, preventing cell interpenetration. Thus, intercellular forces on the i-th bead of one cell due to the interaction with the j-th bead of another cell are

The distance between two interacting beads is r_ij = |r_i,α − r_j,β|, and the corresponding unit vector is = |r_{i, α} − r_{j, β}|, where i-th and j-th beads belong to α-th and, β-th cell, respectively. Adhesive and repulsive forces are characterized by corresponding strengths (K_adh and K_rep respectively) and cut-off ranges (r_c^adh and r_c^rep).

Additionally, we incorporate a polarity vector = (cos θ_α, sin θ_α), in which direction each bead of α- th cell exerts a self-generated motility force of magnitude cv₀ (where c is the coefficient of viscous drag and v₀ is the self-propulsion speed). All beads in each cell possess the same motility direction (Fig. 1A(bottom)), but it differs from cell to cell. Together, the equation of motion governing the overdamped dynamics of the i-th bead of α-th cell becomes

Here, j-th bead belongs to another cell, interacting with the α-th cell. The unit polarity vector for the α-th cell further undergoes a rotational diffusion given by

Here, the angle θ_α defines , and ζ_α(t) is a Gaussian white noise with zero mean and variance 2D_r. We nondimensionalized the model by choosing the units of length and time as l₀ and c/K_s, respectively. Thus, four dimensionless parameters mainly determine the collective dynamics: non-dimensional intracellular pressure , adhesion strength ( = K_adh/K_s), selfpropulsion speed , and rotational noise strength . We used the model to simulate a nearly confluent tissue with 256 cells, each consisting of 50 beads. Starting with random cell centers, the tissue monolayer was simulated using Euler’s method up to 10⁷ iterations with a time step 10⁻³, under periodic boundary conditions (see SI for simulation details and Table S1 for parameter values).

FIG. 1.

III. Results

We first explored how the cell shapes and tissue state qualitatively change with the variations of cell-cell adhesion () and rotational noise strength () (see Fig. 1B). When and are low, cells are mostly elliptical and rounded shaped with noticeable intercellular gaps (Fig. 1B, bottom-left). With an increase in or , regular polygonal cell shapes (like hexagons and pentagons) coexist with rounded shapes, along with some intercellular gaps in the tissue (Fig. 1B, top-left and bottom-right). In contrast, when both and are high, almost all cells become regular hexagonal-shaped with virtually no cell-cell gaps (Fig. 1B, top-right), suggesting cell jamming.

To characterize the above phenomenon, we first evaluated the mean-squared displacement (MSD) of cell center trajectories, defined as 〈Δr²(Δt)〉 = 〈[(t₀ + Δt) − (t₀)]²〉. Here, (t) is the position vector of the a-th cell center at time t, t₀ is a reference time, Δt is the lag time, and the angular brackets denote the average over all cells and over different ensembles. We found that the cell motion is diffusive for low (see Fig. 2A). However, with the increasing , the motion becomes subdiffusive, and eventually, the MSD plateaus at large times, indicating the arrested motion of cells caged by their neighbors. Thus, the monolayer transforms from a fluid to a solid regime with increasing adhesion (when other parameters are fixed). The monolayer also gets solidified with increasing or , and decreasing (Fig. S3).

We further observed that cells spontaneously exchange their neighbors (known as T1 transition [46]) in the fluid phase (see Movie S1 and Fig. 2B). Correspondingly, the cell center trajectories appear diffusive in nature (Fig. 2B). In contrast, the trajectories look caged in the solid phase, as expected from the MSD curves (Fig. 2A). We then computed an effective diffusivity from the long-time MSD behavior [32], defined by D_eff = lim_t→∞,〈Δr²(Δt)〉/4D₀Δt, where D₀ = υ₀²/2D_r can be regarded as the self-diffusivity for an isolated single cell [32]. The order parameter, D_eff, showed a transition with increasing intercellular adhesion (Fig. 2D). We also measured the MSD exponents (β) by fitting the MSD curves at large times with the function: 〈Δr²(Δt)〉 ~ Δt^β. The MSD exponent was 1 for the fluid regime and showed a sharp drop as the tissue solidified (Fig. 2E).

Next, we explored the phase space spanned by the three dimensionless parameters, the intercellular adhesion (), the intracellular pressure (), and the rotational noise strength (). We determined the phase diagrams in the 2D parameter space by varying two of the above parameters (Fig. 3A-C). We identified the fluid phase by the cut-off D_eff > 0.001 and solid by D_eff ≤ 0.001 (cut-off values were chosen following Fig. 2D). Here, we determined the phase diagrams by keeping the cell motility () fixed, though the phase planes can also be constructed by varying (Fig. S4). The 2D phase planes can be combined to visualize a schematic 3D phase diagram (Fig. 3D). Similar to the SPV model prediction [32], Fig. 3D shows that the monolayer gets fludized with increasing persistence time scale for the rotational noise (given by 1/). But, concerning the effect of cell-cell adhesion (), our phase diagram fundamentally differs from Vertex and SPV models [32, 47], where adhesion helps fluidization in a confluent tissue. However, our prediction agrees with the speculated jamming phase diagrams derived from previous experiments [9, 41], which hypothesized that intercellular adhesion solidifies the tissue. In addition, our phase diagram crucially points out that solidification can occur with increasing intracellular pressure, consistent with the observation that high cytoplasmic pressure promotes epithelial integrity [48].

FIG. 2.

Since cells can modulate cortical tension during morphogenesis [49], we calculated the tension in cell-edges, defined as T_i = |l_i − l₀| (ignoring the multiplicative constant), where l_i is the distance between i-th and (i − 1)-th beads in a cell and l₀ is the natural spring length. Deep in the fluid phase, the tension is highly heterogeneous over the tissue (Fig. 4A), while it is homogeneous in the solid phase (Fig. 4B). We also calculated the average tension in the steady state (given by T = (∑_{M, N} T_i)/MN, where M and N are the bead and cell numbers, respectively), and found that it is lower and display higher fluctuations in the fluid phase than the solid phase (Fig. 4C, Fig. S5). Thus, tension fluctuations decrease with the solidification, as recently shown in the ‘Dynamic Vertex Model’ (DVM) with active tension dynamics [7].

FIG. 3.

However, tension fluctuation is an emerging property of our model, while in DMV, tension fluctuation was introduced via a dynamical equation for the edge tension.

We also measured the area fraction (ϕ) over time, defined as ϕ = , where A_α is the area of the α-th cell and L is the simulation box length. The area fraction substantially fluctuates in the fluid phase compared to the solid phase (Fig. 4D), although it increases marginally (from 0.83 to 0.87) from the fluid to solid phase (Fig. 4D). Thus, the tissue remains near confluent throughout the transition, similar to the measured in vivo volume fraction (around 80% − 90%) in dense zebrafish tissues [1]. We next computed the cell shape index, S = = , P_α and A_α being the perimeter and area of the a-th cell, respectively [26, 32]. Note that this shape index can be measured as an emerging property [7, 10, 26, 32, 37], while the ‘target shape index’ is a tunable parameter defined in the energy function of Vertex and Voronoi (SPV) models [26, 32]. We observed that the shape index fluctuates wildly over time in the fluid phase, but its fluctuations diminish in the solid phase (Fig. 4E).

The Vertex and SPV models [26, 32] showed that the average shape index, hSi, serves as a structural order parameter. In DVM and SPV models [7, 32], hSi increases monotonically as the tissue fluidizes. In contrast, our model reveals that hSi behaves non-monotonically with adhesion near the fluid-solid transition, especially for lower (Fig. 4F). Moreover, the dynamical order parameter, D_eff, is also non-monotonic with 〈S〉 (Fig. 4F inset), indicating that both fluid and solid phases can possess the same 〈S〉 value close to the transition. Thus, the shape index alone cannot distinguish between the tissue phases. We then examined how the shape index fluctuates across the transition, as marked by (i), (ii), and (iii) in Fig. 4F (denoting three distinct regions: deep in the fluid phase, near the phase boundary, and the solid phase, respectively). The standard deviation of S () displayed a peak near the transition boundary in region-ii (Fig. 4G). Corresponding tissue configurations are shown in Fig. 4H, and the distributions of the shape index (P(S)) are plotted in Fig. 4I. At region-ii, where σ_S has a peak, the distribution P(S) becomes more skewed and broader than other regimes (see Fig. 4I). Interestingly, elongated polygonal and rounded cell shapes coexist in region-ii (Fig. 4H(ii)), whereas cell shapes are mostly rounded in region-i and polygonal (i.e., pentagonal or hexagonal) in region-iii (Fig. 4H). Moreover, in the solid regime (region-iii), the monolayer resembles a tightly packed foam, where the normalized cell area follows a generalized k-gamma distribution as found earlier in granular material [50], foam [51], and the Vertex model [52] (see Fig. S6).

The cell shape heterogeneity that arises near the transition (Fig. 4G, H) may be linked with underlying dynamic heterogeneity in cell motion. To examine this, within a time window (t₁ to t₂), we measured the total normalized distance covered by individual cell centers, , defined for the α-th cell as , where − . As suggested by the heatmaps of Fig. 4J, different mobilities exist in region-i but in a spatially random fashion (Fig. 4J(i)). However, cells in region-iii mostly have lower mobilities (Fig. 4J(iii)). Significantly, in region-ii, fast-moving and slow-moving cells form clusters (Fig. 4J(ii)), suggesting that cells with distinct mobilities coexist with spatial correlation. This indicates that cells are dynamically heterogeneous, suggesting a glass transition, as found in colloidal glasses [53, 54] and in dense bacterial biofilm [55].

To probe the glassy versus fluidic behavior, we plotted the vector field of cell center displacements for region-i and region-ii, corresponding to lower and higher adhesion strengths, respectively (see Fig. 5A, B). The displacement vectors are random for region-i but exhibit swirling patterns indicating spatial correlations in region-ii. Also, large and small displacements coexist in the vector field of region-ii, suggesting dynamic heterogeneity similar to polycrystalline materials and glasses [53, 56–58]. Notably, similar swirling displacement fields were observed in previous models [32, 59, 60] and epithelial monolayers [3, 12, 13, 61].

FIG. 4.

Another quantity that can capture the glassy behavior is the non-Gaussian parameter defined by α₂ (Δt) = 〈Δr⁴(Δt)〉/(2〈Δr²(Δt)〉²) − 1. It is non-zero if displacement distributions deviate from Gaussian [37, 57]. In Fig. 5C, α₂ is almost zero in region-i (lower adhesion), for which diffusive cell dynamics is expected. In contrast, α₂ has a sharp non-zero peak at the lag time Δt^* in region-ii (higher adhesion), for which cell dynamics is subdiffusive. Accordingly, the displacement distribution P(Δ), measured at Δt^*, is Gaussian in region-i, but it develops broad exponential tails in region-ii (Fig. 5D, 5E). Such non-zero peaks in α₂ and corresponding nonGaussian tails in displacement distributions have been observed in inanimate glassy materials [53, 57, 58], and also in living systems like embryonic tissues [62] and cytoskeletal dynamics in human muscle cells [63].

The broad exponential tail in the displacement distribution (Fig. 5E) represents the motile cells performing random jumps between cages [53, 57, 63]. The lag time Δt^*, at which a₂ shows a peak (Fig. 5C), gives the cage breaking time scale [53, 57, 58]. Indeed, within a time window spanning At^*, a sample trajectory in region-ii displays rare and quick hopping of a cell between cages formed by its neighbors (Fig. 5F(ii)), but the trajectory looks simply diffusive in region-i (Fig. 5F(i)).

FIG. 5.

However, what could be the physical origin of the exponential tails in displacement distributions observed in the glassy regime (region-ii)? We found that the individual diffusion coefficients (D), measured from the time-averaged MSD curves, have a broad skewed distribution that has an exponential-like form (Fig. 5G). Thus, the exponential tails in displacement distributions (Fig. 5H) can stem from the cell-to-cell heterogeneity in individual diffusivities [64, 65]. Suppose individual cells in the ensemble follow independent diffusive dynamics characterized by Gaussian displacement distributions but with an exponentially distributed diffusion coefficient, D. Then, non-gaussian displacements can be achieved as the convolution of Gaussian — We can write the displacement distribution (at a lag time Δt^*) as: , where, . If we assume the distribution of D as p(D) = , where is the mean value of D averaged over all cells, we obtain the form: . Further, note that we numerically obtained individual diffusion coefficients by the definition for the 2-dimensional dynamics. Thus, for the 1D component of displacements (Δx or Δy), we take the mean diffusivity to be to consider the dimensionality. Therefore, if we consider the scaled 1D displacement , the distribution then follows an exponential distribution with exponent 1. Indeed, for various values of and , different non-Gaussian exponential tails with distinct exponents collapse onto a single master curve ( when the scaled displacement is used (see Fig. 5H, I). Thus, represents a mean diffusion length by which displacements can be normalized to obtain a universal form, elucidating the link between exponential tails and the heterogeneous diffusivities of individual cells.

IV. Discussion

Summary and relevance

Our result of tissue solidification with adhesion, though contrasts the prediction of Vertex and Voronoi models [26, 32], are in line with in vitro experiments [12, 40, 41]. This also explains how the downregulation of adhesion molecules leads to in vivo tissue fluidization during zebrafish tail development [1, 8]. Similar to the experiment [1], adhesion-induced solidification accompanies a marginal change in packing fraction in our model. We also found tissue fluidization with increasing self-propulsion speed that captures a recent experiment on motility-driven unjamming in gastrulation [2]. We further show that the average shape index alone cannot be used as a structural order parameter to distinguish the phases due to large shape fluctuations near the transition point. Also, increased adhesion results in dynamic heterogeneity and glassy behavior near the phase boundary.

Our force-based framework, which is conceptually straightforward to implement, relies on a free energy function that obeys the Young-Laplace equation at the level of a single cell in equilibrium (see Section III, SI) — this is distinct from the topological energy functions used in Vertex, Voronoi, and Cellular Potts models [26, 32, 37]. Contrasting these models, we do not assume any target shape parameters since there is no experimental support and such parameters are difficult to physically interpret. Our model exhibits spontaneous neighbor exchange (T1 transitions) unlike the Vertex/Voronoi frameworks, where T1 transitions are performed using an ad hoc criterion of a threshold value on cellular edge lengths [7, 26]. Further, intercellular gaps naturally arise in our model, unlike the Voronoi tesselation framework, where cell-cell gaps are either completely absent [26, 32] or implemented by assuming extra ‘non-physical’ vertices [7]. We also observe the spontaneous formation of multicellular rosettes (where five or more cells share a vertex) [66] (Fig. S7). In contrast, rosettes were formed in Vertex models by implementing an ad hoc ‘edge collapse’ procedure that randomly reduces edge lengths to zero [30].

Connection with other models

To compare with the jamming in passive foams, we studied our model without the self-propulsion and adhesion (υ₀ = 0, K_adh = 0). In a rigid box, we produced jammed states by inflating the cells with increasing pressure (i.e., with increasing packing fraction). We show that the tissue becomes confluent when the normalized area fraction reaches unity (ϕ/ϕ_max = 1) at an average shape index 〈S〉 ≈ 3.81 (see Fig. S8 and Section V, SI). This is similar to the crowding-dependent jamming in the Deformable Particle (DP) model [35]. Interestingly, the Vertex and SPV models [26, 32] also show a liquid-to-solid transition at 〈S〉 « 3.81.

Note that the level of confluency can be described by ϕ/ϕ_max = 1 implies full confluency) [35]. Though the absolute packing fraction cannot reach unity in our particle-based model (unlike Vertex and SPV models), we maintained a near-confluent regime, where (ϕ/ϕ_max) ~ 0.94 − 1 across the transition (corresponding to Fig. 2–4). This is similar to the in vivo experiment on zebrafish tissues that show 80%–90% confluency [1].

However, to compare with Vertex and SPV models, we performed the Voronoi tessellation of space (using the cell centers) to obtain the Voronoi polygons corresponding to the cells (Fig. S9A) and measured the shape indices of the Voronoi cells, 〈S〉_{V or} (see Section VI, SI). We found that the liquid-solid transition takes place at for low values of D_r, similar to the Vertex and SPV models [26, 32], though the transition point deviates from 3.81 for higher Dr (Fig. S9B–C). Moreover, we achieved very high confluency by increasing the rest length of cell edges (represented by Hookean springs) and found that the packing fraction hardly altered during the liquid-solid transition (ϕ ~ 0.891–0.893, i.e., (ϕ/ϕ_max) ~ 0.998 − 1), when adhesion was increased from zero (Fig. S10C-E). The average coordination number also stayed very high and almost constant (Z ~ 5.90 — 5.98) during the transition (Fig. S10F). This indicates an adhesion-dependent active glass transition in the confluent limit, not a crowding-dependent jamming. Thus, tissue solidification arises from dampening the active motility as higher adhesion leads to increased resistance to shear forces between cells. Further, measuring from Voronoi cells again showed that the transition happens at (Fig. S10G), exhibiting the behavior of Vertex and SPV models [26, 32]. Nevertheless, the actual shape indices measured from cell contours are non-monotonic across the transition (Fig. 4F, S10H).

As shown before, in the passive limit, our active model is similar to passive deformable particle (DP) models [35, 44]. However, including ‘active and noisy’ motility is crucial to observe the glass transition. Though a few studies included tunable confluency in the Voronoi framework [7, 67], our implementation is distinct and straightforward, inspired by DP models. There are also a few active models of deformable cells like ours [68, 69]. However, the intracellular pressure in these models depends on the instantaneous cell area, which is similar to vertex-based models but distinct from our constant pressure framework that aligns with in vivo measurement of osmotic pressure in tissues [70]. Recently, a computational package [71] has been developed to simulate various soft materials, including biological tissues, foams, bubbles, etc. However, this framework is overtly complex (with 23 parameters), while ours is a simple coarsegrained framework with a few parameters corresponding to experiments (such as adhesion, intracellular pressure, edge tension, etc.).

Application and future directions

Intriguingly, despite having a distinct single-cell energy function compared to the DP, Vertex, and SPV models, our model still shows the liquid-solid transitions at in certain limits. This may be tied to the argument of isostaticity based on the polygonal tiling of the 2D space, leading to an isostatic coordination number of Z_iso = 5 (corresponding to the for a pentagon) as explained in [26]. However, a detailed correspondence between the vertex-based framework and ours would be an interesting future problem. Since real biological tissues are near-confluent, our model is suitable for studying the effect of confluency on active glassy behaviors. In nearconfluent systems, it is worthwhile to study the interplay of two mechanisms of solidification: crowding-dependent jamming (in sub-confluent regimes) and adhesion-driven or tension-dependent rigidity transition (in full confluent regimes) [72]. For instance, a recent study on Xenopus development [73] shows the signatures of both mechanisms. Another work has discovered that periodic assemblies of myosin motors interconnect with the actin cortex and act like a chain of sarcomeric units in each epithelial cell [74]. Such an arrangement resembles our bead-spring loops representing the cell cortex. The authors found that the distance between sarcomeric units changes upon perturbation of the actin cortex, and it relaxes to the original state when the perturbation is removed. Our model could be suitable to understand the effect of such mechanical perturbation on cell shape and tissue rigidity. Moreover, a detailed continuum model has recently explored how T1 transitions emerge from stochastic formation and dissociation of cell-cell adhesive bonds at the level of a few cells [75]. Since T1 transitions spontaneously emerge at the tissue level in our coarsegrained model, an interesting question would be how the T1 transitions depend on confluency and adhesion.

Supplementary Information: Role of intercellular adhesion in modulating tissue fluidity

I. Model Details

A. Force inherent to an isolated cell

Motivated from previous studies [43–45], a single cell is modeled as a closed loop of beads connected through elastic springs with stiffness constant K_s and natural spring length l₀ (Fig. S1). In our simulation, each cell is made up of 50 beads. Each bead experiences tangential tension forces from the two adjacent neighbor springs, and an internal pressure force, Pl₀, directed (outward) normal to the line tension of the given spring. (Fig. S1).

FIG. S1.

In a coarse-grained view, these intracellular forces describe the roles of the actomyosin cortex and cytoplasmic pressure that give the cell its structural integrity. The intracellular force for the i-th bead can be explicitly written as follows:

Here, l_i is the bond length between i-th and (i − 1)-th bead pair, and l_i+1 is the bond length between (i + 1)-th and i-th beads, respectively. Also, and are the tangential unit vectors, and and are the corresponding outward normal unit vectors perpendicular to and , respectively (see Fig. S1). The first two terms represent the total spring-tension force, and the third term represents the outward pressure force on the i-th bead.

B. Effect of cell-cell interactions

The intercellular force arising due to cell-cell interactions consists of two parts: (i) the spring-like attractive force, which resembles the cell-cell adhesion through adhesive proteins (like E-cadherins), and (ii) the spring-like repulsion that prevents the interpenetration among neighboring cells. Two beads of two different cells would feel the force of adhesion or repulsion only when the Euclidean distance between them is below the cut-off range of the attractive or repulsive forces. Thus the interaction force on the i-th bead of α-th cell due to the j-th bead of β-th cell is as follows:

Here, is the unit vector of the Euclidean distance between the two interacting beads. K_adh and are the strength and cut-off range, respectively, for the attractive force, and K_rep and are that of the repulsive force.

Now, considering all the contributions from other cells within the interaction ranges, the total force of interaction on the i-th bead of the α-th cell would be :

C. Active force term via self-propulsion and Equations of motion

In addition to the intracellular and intercellular forces, cells can also move via self-propulsion. As in the selfpropelled particle (SPP) model [19] and self-propelled Voronoi (SPV) model [32], we also incorporate a polarity vector for every cell. In the direction of , each bead of a cell exerts a self-generated motility force of magnitude cυ₀ , where c is the coefficient of viscous drag, and υ₀ is the self-propulsion speed. Thus, in each cell, all the beads have the same polarity direction, but the polarity direction differs from cell to cell. All together, the overdamped equation of motion for the i-th bead of the α-th cell is:

And, the unit polarity vector undergoes a rotational diffusion described as:

where, ζ_α(t) is a Gaussian white noise with zero mean and variance 2D_r, which specifies the timescale for reorientation of the polarity vectors as 1/D_r.

Thus Eq. (4) and Eq. (5) are the governing equations of motion in our model.

D. Non-dimensional Parameters

To non-dimensionalize our governing equations, we take the unit of length as l₀ and the unit of time as (where c is the coefficient of viscous drag and K_s is the spring constant for a single cell). Thus, the non-dimensional position vector () would be :

And, non-dimensional time () would be :

Therefore, replacing the position vectors and time with the non-dimensionalized values in Eq. (4), we obtain

Here, K_int denotes the strength of adhesion or repulsion K_adh or K_rep in a compact notation. Also, , and are dimensionless pressure, interaction strength, and self-propulsion speed, respectively. Note that Eq. (8) is a non-dimensional form of Eq. (4). Treating Eq. (5) in a similar way, we get

and

Where, is the non-dimensional rotational noise strength. Thus, the main non-dimensional parameters of our model are , , and (ignoring , since we kept .

II. Simulation Details

We consider a system of N = 256 cells in a square box of size L = 460l₀ . Each cell contains M = 50 beads. The bead positions are updated by numerically integrating Eq. (4) and Eq. (5) using the Euler method. The time step is 10⁻³ and the total number of iterations is 10⁷. The simulations are carried out using a periodic boundary condition.

A. Initialization procedure

Generating nearly confluent tissue configuration

To populate the desired number of cells in a given square box, we choose a smaller initial radius for each cell. In our simulation, the initial radius of each cell is taken as a = 10l₀, and the box length is chosen to be L = 460l₀. To avoid cell-cell overlap during the random generation of cells, the minimum separation between two cell centers is chosen to be 2a + . In a square box of fixed length, L, we generate two random numbers from a uniform distribution as the X and Y coordinates of each cell center. The cell centers are then checked to meet the criterion of lying inside the box and also have no overlap with already generated cells. If there is an overlap, the cell center coordinates are generated repeatedly till the criteria are met. The 50 beads that make up the cell boundaries are then positioned in an equidistant manner along the perimeter of each cell such that the spring length between successive beads is 2na/50.

Equilibrium dynamics of the tissue configuration

The initially generated cell-system is then evolved by integrating Eq. (4) with , i.e., we evolve the passive system without Eq. (5). During this time window, the cells grow, and by interacting with each other, they approach a nearly confluent equilibrium configuration (with an area fraction above 0.8). This equilibrium configuration then serves as an initial configuration for simulations of the active system. In all our simulations, we equilibrate the randomly generated configurations up to the total iteration of 5 × 10⁴. Also note, to avoid any global drift of the system, we integrate Eq. (4) in the center of mass reference frame, i.e., we use the transformed position vectors for each bead, r^′i,α = r^i,α − R^cm, where R^cm is the position of the center of mass of the entire tissue monolayer.

B. Dynamics of the active system

Simulation of the active system

After the passive system attained an equilibrium configuration (described above), we reset the time to zero. The final equilibrated configuration of the passive system is then used as the initial configuration for the active dynamics of the tissues. The system is then evolved for a further 10⁷ iteration steps with a non-zero self-propulsion speed , and the dynamics is governed by both Eq. (4) and Eq. (5) simultaneously.

Measurement of observable quantities

We found that after 10⁵ iterations, the active system entered into a steady state, which was determined by observing the time evolution of the dynamical or structural observables like instantaneous area fraction and the shape index (as defined in the next section). The steady-state was confirmed when the observables fluctuated around a well-defined average value over time (as in Fig. 4(C-E), main text). We took this approximate steady-state reaching time as the reference time (t⁰ = 10⁵ iterations) for measuring the dynamical quantities like MSD and the non-Gaussian parameter (as defined in the next section).

Code details

The main simulation as well as codes for analysing the data are all written in FORTRAN 90. We implemented the ‘cell linked-list’ algorithm [76] (which provides a significant speed-up) in the neighbor search process to compute the intercellular forces. The codes are available in the GitHub repository: https://github.com/PhyBi/Collective-Cell-Dynamics/tree/legacy.

TABLE S1.

III. Comparison of Our Model with the Topological Energy Function of Vertex and Voronoi Models

The free energy of our model fundamentally differs from the topological energy function of Vertex/Voronoi-type models. For instance, consider a single isolated cell at equilibrium. Under an isotropic pressure, the cell will attain the equilibrium shape of an n-sided regular polygon having n beads and n springs (with spring constant K_s and rest length l₀). Following Fig. S2, the free energy is:

FIG. S2.

where the first term is the work done by the pressure, the second term is the spring energy, and we invoke the geometric relation, Δl = (r − r₀)sin(π/n) (see Fig. S2). Differentiating Eq. (11) with respect to r, we obtain at equilibrium:

Using the geometric relation, l₀ = 2r₀sin(π/n) (see Fig. S2), Eq. (12) can be reduced to

where T = K_s (2Δl) can be considered as the local surface tension at each edge. The Eq. (13) is basically the YoungLaplace (YL) equation, which is thus satisfied in our model. In contrast, the vertex energy function for a single cell is given by

where p and a are perimeter and area of a cell, while p₀ and a₀ are target perimeter and target area parameters, respectively. The first term in Eq. (14) is similar to the surface spring energy in our model. However, the 2nd term is the ‘area elasticity,’ which represents a pressure dependent on the cell’s instantaneous area (pressure, P = K_A(a − a₀)). This is inconsistent with the YL equation, which is arguably responsible for the tendency of rounding of a single cell during detachment [77, 78]. The YL equation was also invoked to quantify adhesion [78, 79], and it formed the basis of inferring tension at cellular edges from analysis of tissue images [80, 81]. Since the vertex energy function does not obey the YL equation, it can lead to skewed non-rounded shapes for a single cell at equilibrium, depending on independent choices of ao and po . Thus, a single isolated cell is poorly defined in this framework. An explicit ‘pressure’ term must be incorporated into the energy function to make it consistent with the YL equation (see [82]).

IV. Measured Observables

A. Mean Square Displacement (MSD)

To analyze the cell motions in the monolayer, we measured the mean square displacement (MSD) of the cell centers. The MSD is defined as

where

is the position vector of the center of mass of the α-th cell at time instant t. The angular bracket 〈...〉 denotes the average over all the cells and over different ensemble runs.

B. Effective Diffusivity

From the MSD, we determined the effective diffusivity D_eff as a dynamical order parameter that differentiates the solid and fluid phases and is defined as

where D₀ = υ₀²/2D_r is the self-diffusivity for an isolated single cell [32]. D_eff is estimated from the slope of the linear fits to the MSD curves at long time (from 5 × 106 to 10⁷ iterations).

C. Non-Gaussian Parameter

The non-Gaussian parameter (α₂(Δt)), in 2D, is defined as [37] :

where Δ_r is as defined in Eq. (16). This parameter includes the 2^nd and 4^th order moments of the distribution of Δ_r. The parameter value is zero for a Gaussian distribution and becomes non-zero when the distribution deviates from the Gaussian. The higher-order moments capture the deviation of the distribution, especially in the tail part, which cannot be accounted for by MSD.

D. Area Fraction

The area fraction in our model can be calculated in the following two ways. First, the area fraction (ϕ) of the tissue monolayer at a particular time is given by

Here A_α is the area of the α-th cell, N is the cell number, and A_box = L², is the area of the simulation box. This definition considers the beads as point particles. However, since we have a cut-off range of the repulsion interaction , each bead can be treated as a disc of diameter , contributing a factor of (π/2)(/2)² to the single cell area. Thus, according to this factor, we can have a higher value of the area fraction given by

Here M is the bead number per cell. We calculated that ϕ ranges 0.85 − 0.88 across the extreme fluid and solid phases (corresponding to Fig. 2 of our manuscript), while ϕ^′ ~ 0.93 − 0.95 during the transition. Both these estimates show high confluency, and the area fraction changes negligibly across the transition, as found during the liquid-to-solid transition in dense zebrafish tissues [1].

E. Shape Index

Following previous studies [26, 32], the cell shape index (S) is defined as:

Where P_α is the perimeter and A_α is the area of the a-th cell at a particular time instant. S denotes the averaged shape index over all the cells at a particular instant. And, 〈S〉 indicates the average of S over time and over different simulations (ensemble average). While calculating 〈S〉, we averaged over the time window from 5 × 10⁶ to 10⁷ iterations when the system is deep in the steady state, as confirmed from the temporal evolution of shape index (S) (see Fig. 4E, main text).

We also measured 〈S〉_{V or}, the averaged shape index obtained from the Voronoi tessellated polygons. Following standard protocol [35], we performed the Voronoi space tessellation from the actual cell center data in a simulated tissue, to obtain the corresponding Voronoi polygons (see Fig. S9A). Then we calculated the shape index, S_{V or}, for each of the Voronoi polygons using Eq. 21. Note that, 〈S〉_{V or} represents the ensemble average as well as the time average in the steady-state (as stated above for measuring 〈S〉).

F. Edge Tension

For the i-th bead of a particular cell at an instant, the edge tension is estimated as the extension of the corresponding spring length (apart from a multiplicative constant):

where, l_i is the instantaneous length of the spring between i-th and (i − 1)-th beads of the cell and l₀ is the rest length. Then we calculated the average tension over all the beads of a cell and over all the cells as below:

G. Coordination Number

By definition, the number of cells with which a particular cell is ‘in contact’, is the coordination number (Z) of that particular cell [83]. Now, in our model, a particular cell, say the a-th cell, would be ‘in contact’ with the β-th cell and vice versa, if any two beads of the two cells come within the interaction ranges, either or . Thus, Z of a particular cell at a particular time instant is the total number of cells that are interacting with it through either repulsion or adhesion.

V. Jamming Onset in the ‘Passive Limit’ of Our Model

The ‘passive-tissue’ limit of our model essentially means the absence of any active self-propulsion vector (υ₀ = 0). In this limit, we studied the system without cell-cell adhesion (K_adh = 0) in order to compare with previous studies on the packing of soft grains [35, 43]. The intercellular interactions are purely repulsive in this case. We considered a system of N = 256 cells in a square box of size L = 460l₀, each cell containing M = 50 beads (same system size as described in the ‘Simulation Details’ section). We followed the same initialization protocol stated in the ‘Simulation Details’ section. Here, we performed the simulations using a ‘stiff’ or ‘rigid’ boundary condition. The stiff/rigid boundary condition was implemented as follows: At a particular time instant t, if any bead of a cell touches the boundary lines of the square box, i.e., if the Euclidean distance between the bead and any of the boundary lines becomes smaller than the range of repulsion , then the bead position was not updated, and the bead was kept fixed to its previous position at the time instant t − dt (where dt is the integration timestep).

The system was then evolved by integrating Eq. 4 with = 0 and (i-e., without Eq. 5) using the Euler method with integration time step 10⁻³. We equilibrated the system by iterating up to 10⁵ steps. The system then reached an equilibrium configuration with a particular level of confluency, depending on the value of the intracellular pressure . With increasing intracellular pressure, the equilibrium configurations slowly went from sub-confluent to highly confluent. Thus, we took a series of values to generate configurations with increasing packing fractions (see Fig S8(A-C)), similar to a protocol used in prior studies [43].

Finally, for each value of , we calculated the observable quantities like cell shape index, area fraction, coordination number, etc, by averaging over 50 independent simulation runs. During the calculations of the quantities, we considered only the cells in the bulk, not those in contact with the box boundaries. Notably, we observed the onset of crowding-dependent jamming in our passive model similar to the Deformable Particle (DP) model [35] published before (see Fig. S8(D-E)). As shown in Fig. S8D, the tissue reaches the confluence when the normalized area fraction reaches unity (ϕ/ϕ_max = 1.0 with ϕ_max = 0.883), at the measured shape index 〈S〉 ≈ 3.81. The average coordination number Z ~ 5.7 indicates a jammed confluent tissue at 〈S〉 = 3.81 (Fig. S8E). Note that this value of Z at the confluence is higher than the isostatic contact number for the Vertex model (Z_iso = 5), as argued from the Voronoi tesselation of the 2D space [26]. Interestingly, the Self-propelled Voronoi (SPV) model [32] also shows a liquid-to-solid transition at 〈S〉 ≈ 3.81 in the zero motility limit. Also note that our model, despite having a distinct single-cell energy function (in equilibrium) compared to the DP and Vertex models, still shows the jamming onset at 〈S〉 ≈ 3.81. This indicates that the jamming onset can be explained from the argument of isostaticity based on the polygonal tiling of the 2D space (since Z_iso = 5 corresponds to the 〈S〉 = 3.81 for a pentagon), as described in [26].

VI. The Confluent Limit of Our Model

As discussed in the previous section (Section V of the SI), the level of confluency may be described by the average area fraction normalized by the maximum area fraction in a parameter regime (ϕ/ϕ_max). The tissue becomes confluent as ϕ/ϕ_max → 1. Note that a particle-based description like our model cannot have ϕ_max = 1 since two adjacent cells have distinct, non-overlapping adjoining edges unlike the Voronoi-tessellation-based models, where adjacent cells share a common edge with ‘zero-thickness’. In our case, ϕ_max → 1 only when the bead number per cell approaches infinity. Nevertheless, 100% confluency is an idealization since in vivo tissues still show 80%-90% confluency (see Fig. 2e, 2h in [1]).

In the main paper, we kept the tissue at a near-confluent regime (where (ϕ/ϕ_max) ~ 0.94 − 1 across the transition). We calculated that ϕ ranges 0.83 − 0.88 across the fluid and solid phases (corresponding to Fig. 2–4, main text), while ϕ^′ ~ 0.91 − 0.96 during the transition (note the slightly different definition of packing fractions, ϕ, and ϕ^′, as given in Eq. 19–20 in SI). These estimates show near-confluency, and the area fraction hardly changes across the transition.

As mentioned above, since our model does not attain ϕ = 1, we cannot directly compare the average shape indices of cells with the shape indices of the Voronoi cells in the Vertex and SPV models. Thus, for comparison, we performed the Voronoi tessellation of space (using the cell centers) to obtain the Voronoi polygons corresponding to the cells (Fig. S9A), and we then measured the shape indices of the Voronoi cells (〈S〉_{V or}, as defined in Sec. IV E, SI). Interestingly, in this near-confluent regime, 〈S〉_{V or} measured from the Voronoi polygons shows that the liquidsolid transition takes place at 〈S〉_{V or} ≈ 3.81 for the low values of rotational noise strength (D_r), similar to the behavior of Vertex and SPV models [26, 32] (see Fig. S9B,C). However, for higher D_r, the phase boundary shifts from 〈S〉_{V or} = 3.81 (Fig. S9C). Thus, our model replicates the behavior of the Vertex model in certain regimes if we treat 〈S〉_{V or} as a structural order parameter. But, note that the actual shape index 〈S〉 measured from cell contours are nonmonotonic near the phase boundary (Fig. 4F, main text), contrasting the behavior of the Vertex model.

Furthermore, as our model has tunable confluency depending on parameter variations, we could achieve high confluency by increasing the natural spring length (l0 from 0.1 to 0.14) and explored the effect of increasing adhesion from zero to a much higher value (see Fig. S10(A-B)). Though the packing fraction did not alter up to the second decimal place (ϕ = 0.891 − 0.893), the MSD of cell centers still showed a liquid-solid transition, where the cell dynamics changed from sub-diffusive to caged (plateau) behavior with adhesion (Fig. S10C, D). The transition points are located based on a cut-off value of the MSD exponent (see Fig. S10D, where the cut-off exponent was set to 0.2). Notably, the average packing fraction (both ϕ and ϕ^′) hardly changed across the liquid-solid transition, suggesting a rigidity-driven glass transition, not a crowding-dependent jamming (Fig. S10E). Also, note that the values of the area fraction in Fig. S10E are well above the critical packing fraction for jamming transition in passive foams — in previous simulations of the ‘deformable polygon model’ and ‘dynamic vertex model,’ the critical packing fraction for jamming was reported to be ϕ_c ~ 0.83 in 2D (see Fig. 2f in [7] and Fig. 4a in [84]). Thus, for the chosen parameters of our model, the tissue is already ‘jammed’ without adhesion, showing sub-diffusive dynamics with active cellular motility (Fig. S10C). However, with increasing adhesion, solidification arises from dampening the active motility of cells, which is caused by the increased resistance to shear forces between cells due to higher adhesion.

Moreover, the solidification in the highly confluent state accompanied a negligible change in the average contact number per cell (Z changes from 5.90 to 5.98 in Fig. S10F corresponding to the red and blue curves, representing the confluent limit). In this regard, a simple counting of constraints for passive non-adhesive spheres shows that the critical isotactic contact number for jamming is Z_c ~ 4 in 2D. Note, the degrees of freedom for N cells is n_DOF = 2N, while the number of constraints is n_c = ZN/2 as each contact is shared between two cells. This leads to Z_c = 4 when n_DOF = n_c (see [85]). However, another constraint counting argument based on the Voronoi tesselation of the 2D space suggests that the isostatic value of Z at the confluence should be , as explained for the Vertex model [26]. Importantly, as shown in Fig. S10F, our observed contact number in the confluent limit is above both the estimated isostatic values (Z_c = 4 and ). Thus, the solidification is essentially an adhesion-dependent active glass transition, not crowding-dependent jamming.

Additionally, in the confluent limit, we measured both the average shape indices from the Voronoi cells (〈S〉_{V or}) and from the actual cell contours (〈S〉) (Fig. S10G, H). Notably, the liquid-solid transition boundaries (pinpointed by arrowheads in Fig. S10D) closely coincide with 〈S〉_{V or} ≈ 3.81 (Fig. S10G), again exhibiting the Vertex and SPV model-like behavior [26, 32] for the given parameter choices. Nevertheless, the actual shape indices still show nonmonotonic behavior across the transition (Fig. S10H), similar to what was observed before (as in Fig. 4F, main text).

VII. Supplemental Figures

FIG. S3.

FIG. S4.

FIG. S5.

FIG. S6.

FIG. S7.

FIG. S8.

FIG. S9.

FIG. S10.

Additional files

Movie S1. Fluid tissue phase (low adhesion regime). The steady-state dynamics of the tissue exhibiting a liquid-like phase at (low) adhesion strength = 8.3 × 10⁻⁵. A group of neighboring cells are marked in red and green to highlight a T1 transition. Also, note the freely diffusive dynamics of the cells. Parameter values: = 0.2, = 16.6 × 10⁻³, = 5.2 × 10⁻⁴. Other parameters are from Table S1.

Movie S2. Solid tissue phase (high adhesion regime). The steady-state dynamics of the tissue exhibiting a solid-like phase at (high) adhesion strength, = 0.25. A randomly chosen cell is marked red to highlight that its motion is completely caged by its neighbors. Rest of the parameter values are the same as in Movie S1.

Acknowledgements

DD and SR thank DBT (Government of India, Project No. BT/RLF/Re-entry/51/2018) and IISER-Kolkata for financial support. We also thank Dapeng Bi (Northeastern University), Dibyendu Das (IIT Bombay), Anirban Sain (IIT Bombay) and Somajit Dey (IISER Kolkata) for their useful discussions.