The cellular Potts model on disordered lattices

Abstract

The cellular Potts model, also known as the Glazier–Graner–Hogeweg model, is a lattice-based approach by which biological tissues at the level of individual cells can be numerically studied. Traditionally, a square or hexagonal underlying lattice structure is assumed for two-dimensional systems, and this is known to introduce artifacts in the structure and dynamics of the model tissues. That is, on regular lattices, cells can assume shapes that are dictated by the symmetries of the underlying lattice. Here, we developed a variant of this method that can be applied to a broad class of (ir)regular lattices. We show that on an irregular lattice deriving from a fluid-like configuration, two types of artifacts can be removed. We further report on the transition between a fluid-like disordered and a solid-like hexagonally ordered phase present for monodisperse confluent cells as a function of their surface tension. This transition shows the hallmarks of a first-order phase transition and is different from the glass/jamming transitions commonly reported for the vertex and active Voronoi models. We emphasize this by analyzing the distribution of shape parameters found in our state space. Our analysis provides a useful reference for the future study of epithelia using the (ir)regular cellular Potts model.

---

1 Introduction

Collective cell migration plays an essential role in many biological settings, ranging from morphogenesis,^1–5 to wound healing,^6–9 to cancer metastasis.^4,10–15 Confluent cell monolayers, such as those found in epithelial tissues, have attracted considerable attention from the modeling community. Firstly, their effective two-dimensional (2D) nature makes them relatively simple to study, yet these tissues can exhibit a wide variety of collective^9,10,16,17 and rheological^18–23 behaviors. Secondly, the fact that epithelia form the outer surfaces of our organs makes them susceptible to cancer^24,25 and play a role in a wide variety of diseases,^26–28 including pulmonary fibrosis²⁹ and asthma.¹⁹ This gives the study of epithelia a direct biomedical relevance.

Over the past decades, various quantitative (computational) models have been developed that capture the principal features of cell populations, without being overburdened by complexity.^30,31 Many are agent-based in nature and have seen concurrent use in the study of active matter.³² For tissues and cells, these models include lattice-based methods like cellular Potts model (CPM)^33–40 and cellular automata (CA),^41–44 particle-based,^3,45–48 vertex,^18,49–52 Voronoi,^53–56 phase-field,^57–63 and Fourier-contour models.⁶⁴ We should emphasize that the constraint of confluency is an important difference between particle-based models and models that have this feature built into their description, e.g., the vertex, Voronoi, and CPM. Evolutionary dynamics can also be included in models to account for (cancerous) mutations and invasion.^65–69 We refer to ref. 11, 30, 57 and 70 for overviews of the different models available and their respective ranges of application. Ref. 11 provides a detailed list of open-source implementations of these models, should one wish to experiment. Lastly, ref. 71 and 72 have made comparisons between different models, which included their ability to describe cell shapes and the way cells exchange neighbors.

In this paper, we will limit ourselves to the CPM, which has good performance on these two points, i.e., cell-shape descriptiveness, and natural neighbor exchange. These abilities of the CPM, naturally lead to a greater computational cost when compared to the models in which cells are described as single particles.^71,72 However, the method is generally considered efficient and has seen widespread adoption. For example, the CPM has been successfully used to study cell sorting and rearrangements,^33,73,74 chemotaxis,^75–77 topotaxis,⁷⁸ cell migration patterns,^79–83 wound closure,^84,85 tumor growth and invasion,³⁷ and the interactions of cells with extracellular matrix.^86–88 Several open-source packages are available for using the CPM^89–94 and the approach has been extensively reviewed, e.g., see ref. 34–36 and 95. However, despite its popularity and qualities, the method is known to suffer from lattice artifacts.³⁶ That is, in certain regimes of model parameters, the shape of cells and macroscopic properties of the system are strongly affected by the symmetries of the underlying lattice. Some of the known issues have been addressed. For example, cells can be fragmented at a sufficiently high rate of cell membrane fluctuations (low surface tension). This was resolved by Durand and Guesnet⁹⁶ by adding an efficient connectivity check. In addition, a node-based version of CPM has been recently proposed with the goal of reducing lattice artifacts.⁹⁷ This model describes the cells as polygons and tracks their vertices.

Here, we will take a different route toward removing/reducing lattice artifacts. We employ an irregular and on-average homogeneous lattice to support our CPM, which derives from a separate simulation of a fluid-like state via Voronoi tessellation. Using this supporting lattice, we study how the absence of long-range order in the lattice affects the transition between (disordered) fluid-like and solid-like states that can occur in model epithelia. Such changes of state are experimentally reported to include‡ jamming, rigidity, and glass transitions,^{14,19,98–104} and reproduced in a variety of vertex- and Voronoi-based computational studies,^{18,53,105–107} which also include polydisperse systems.¹⁰⁸ Here, we should mention that in Voronoi-based models, although dynamical transitions are observed,⁵³ they are known to be absent in the athermal version of the model.¹⁰⁹ These transitions have attracted attention as key players in the development of the aforementioned tissue-related diseases. The change from fluid-like to solid-like is commonly referred to as a phase transition in the literature^{14,19,103,110,111} and we adopt the nomenclature. Note, however, that biological tissues are intrinsically out of equilibrium.

Disordered arrested dynamics and fluid-to-solid transitions have been reported for the CPM.^112–114 However, Durand and Heu¹¹⁵ used the CPM to study soft cellular systems and instead found an order-to-disorder phase transition. We revisit the work by Chiang and Marenduzzo¹¹² and show that lattice artifacts present in their systems are removed by our irregular-lattice CPM. Our results further demonstrate that for the CPM with their (simple) Hamiltonian, there is an order-to-disorder transition rather than a disordered solidification. This transition has all the hallmarks of a first-order phase transition from a fluid to a hexagonal solid. We verified the nature of the transition by examining in detail the geometric features of the cells in the tissue and their neighborhoods. Such features include the (distribution of) the isoperimetric quotient and the circularity. We have also studied how our irregular lattice compares to the use of a hexagonal one and find that the latter has spurious dynamics in the hexagonal crystal state.

The benefits of using a CPM on an irregular lattice come at only a small computational overhead compared to using the regular CPM. This makes it a suitable alternative for the study of real biological tissues, which we aim to pursue in future work.

2 Methodology

In this section, we cover the main features of our variant of CPM: the creation and characterization of irregular lattices, and the means by which we have modified the traditional CPM to work with these lattices. We also discuss the various means, including the mean-squared displacement (MSD) and isoperimetric quotient, by which we characterize the outcomes of our simulations. We further provide our standard choices for the system parameters and simulation ranges that we have considered.

2.1 The cellular Potts model

We introduce the basic CPM algorithm here so that the background for our extension is set. Assume that we have a lattice with N_s sites and N_c cells, we can define a function σ: {1,…,N_s} ↦ {1,…,N_c} that describes the configuration of the cells on the lattice, uniquely. That is, σ(i) indicates the cell index, to which the site i of the lattice belongs. As the development of CPM was inspired by the Potts model¹¹⁶ – originally used to study spin systems – let us call σ(i) the spin of site i. Note that the real biological system of interest for the CPM has nothing to do with magnetism and the spin indices.

Now that we can describe a configuration using σ, we can specify the Hamiltonian that gives the total energy of the system, H = H(σ,[scr P, script letter P]), where [scr P, script letter P] represents the set of all physical parameters of the system, e.g., surface tension and target cell area. In general, the lattice can be of any dimension, but we will restrict ourselves to 2D CPMs here. These are suited to describe cell monolayers, which form a class of epithelia, found in the cells lining blood vessels and alveolar sacs of the lung.¹¹⁷ The simplest form of Hamiltonian that is used on 2D regular lattices is given by

The first term gives the total interaction energy between the cells that comes from the surface tension between the cell membranes. The indices i and j indicate the lattice sites and the summation is carried out over the site pairs within each other's interaction neighborhood. The Kronecker delta δ_ij (1 if i = j and 0 if i ≠ j for any indices) is used to indicate that only adjacent sites that belong to cells with different spins contribute to the surface tension. The factor α indicates the surface tension and is typically considered uniform between cells, though the method can be used to study mixtures of different cell types as well.^33,73,118

The second term indicates that each cell has a preferred (or target) area A₀. Departures of the instantaneous area a_σ away from A₀ are penalized using a Hookean potential with spring constant 2λ. This choice models the tendency of cells to have a constant volume, as well as their connection to their neighbors within an epithelium. That is, any change in shape (elongation/shrinking out of the plane) will lead to a change in the in-plane area (volume conservation). It is assumed that any deviation from the natural shape is associated with an energy cost, which at the lowest order would be quadratic.

The ‘dynamics’ of the CPM is propagated using a Monte-Carlo (MC) approach. This consists of trial moves weighted by the change in the energy. A trial move involves changing the spin of a randomly chosen site to that of its neighbors. In brief, the basic implementation is

(1) Store the energy of the current configuration E_old = H(σ_old,[scr P, script letter P]) and choose a lattice site randomly, say, site i. We call it the candidate site.

(2) Choose a site from its surroundings§, [scr N, script letter N]_c(i), say, site j. We call it the invading site.

(3) If σ(i) = σ(j), go to step 1. Otherwise, temporarily change σ(i) into σ(j), and calculate the new value of energy, E_new = H(σ_new,[scr P, script letter P]). We call this configurational change an attempt.

(4) Call the change in energy ΔE = E_new − E_old. When ΔE ≤ 0, accept the attempt. When ΔE > 0, accept the attempt with the probability

or reject it with probability 1 − P_acc. Here, the energy difference is normalized by the thermal energy, where k_B is the Boltzmann constant and T is the temperature. These steps will repeat until a user-defined maximum number of iterations is reached.

It should be noted that the use of thermal energy in the algorithm, reflects its origin as a tool to study spin systems. The interpretation of the temperature in the context of a cell membrane is to set a rate, at which the various (internal) cell activities change the boundaries. This rate can be given the interpretation of a time scale for membrane fluctuations,^36,74,96 provided only local trial moves are used. Nonetheless, the CPM's dynamics do not represent the evolution of the system, as would follow from say a Langevin description. However, we should note that there have been attempts to reconcile the ‘time’ in the CPM with a physically meaningful time through the introduction of Poissonian statistics.¹¹⁹

According to the algorithm that we have described above, the site pairs which belong to the same cell, are also allowed to be chosen. Such picks are always disregarded for updates after checking that the spins are the same. However, picking these pairs comes at a computational cost. It is possible to identify all ‘allowed site pairs’ in linked list data structures¹²⁰ and only choose from the elements of these lists. However, in practice, for our typical parameter choices, it turned out that searching and updating the linked lists was computationally disadvantageous. Hence, we utilized the above algorithm. Note that linked lists should become more efficient when the number of border sites is considerably less than the total number of lattice sites. Examples of this include the simulation of single cells¹²¹ and non-confluent cell populations.^122,123

2.2 Using (Ir)regular lattices

The above algorithm for site updates holds for any underlying lattice, provided H(σ,[scr P, script letter P]) and [scr N, script letter N]_c(i) for 1 ≤ i ≤ N_s are well-defined. As mentioned in the introduction, CPM has predominantly been applied to regular lattices, i.e., square and hexagonal lattices. We are aware of only one instance of a study of an irregular lattice mentioned in the literature;¹²⁴ in ref. 125 a graph is shown of what appears to be a study performed on an irregular lattice, but the original text could not be obtained. For the sake of completeness, we should also mention the node-based version of CPM,⁹⁷ which describes cells through surfaces rather than through volumes, making it principally different from the standard CPM.

When moving to irregular lattices, we will assume that for any lattice site i, we have [scr N, script letter N]_c(i) = [scr N, script letter N](i). In addition, we have assumed that the sites i and j are neighbors to each other, if and only if the polygons that contain these sites have at least one vertex in common. This simple rule on square lattice leads to Moore neighborhood which considers 8 neighbors for each site, see Fig. 1(a). On a hexagonal lattice, this leads to six neighbors, while the number of neighbors will vary per site on an irregular lattice, see Fig. 1(b) and (c), respectively. All of the neighborhoods considered thus far are in contact with the central cell. This is intuitive, as the interaction term considered in the CPM Hamiltonian describes surface tension. For discussion on other definitions of a neighborhood (for regular lattices) we refer to, for example, ref. 96 and 126.

Fig. 1 The definition of a neighborhood on different lattice types. (a) The Moore neighborhood on a square lattice, where the central site is blue and its 8 neighbors are indicated in yellow. (b) On the hexagonal lattice, there are 6 nearest neighbors to a central site. (c) On the irregular lattice, the sites having at least one shared vertex are neighbors.

To work with an irregular lattice, the Hamiltonian of eqn (1) needs to be modified. As the surface energy is proportional to the contact length of cells, we introduce a weight factor in the interaction term of the Hamiltonian. This weight ensures that the surface-energy penalty is proportional to the contact length of neighboring sites. The Hamiltonian then reads

in which the notation is mostly the same as in eqn (1). Here, the first summation on the right-hand side runs over all neighboring sites based on the aforementioned neighborhood, and w_ij is the weight factor shared between sites i and j. This weight is defined

w_ij = l_ij/l̄, (4)

where l_ij is the length of the contact edge between the sites i and j, and l̄ is the average length of the edges taken over the entire lattice. Taking w_ij = 1, which is appropriate for regular lattices, the Hamiltonian of eqn (3) reduces to that of eqn (1).

We generate our irregular lattices from a set of points using Voronoi tessellation. There are many ways in which to choose the generating points, e.g., by choosing these randomly on the plane using a uniform distribution. However, this leads to the presence of many small Voronoi cells and several large ones.^127,128 Here, we want our lattice sites to have roughly the same size and number of neighbors, whilst maintaining an isotropic character to the distribution of points. Therefore, we choose to base our lattice on the center of mass (CMS) of uniformly sized particles in a fluid phase. A regular CPM can be easily implemented and can serve to generate such a configuration for a suitably chosen value of α = 0.8, which places the configuration in the fluid phase. We performed a large-scale simulation to obtain an equilibrated fluid, see Table 1 for our choices, by which we obtained approximately 200² lattice centers. This is illustrated in Fig. 2(a)¶.

Table 1 Parameter choices and details of the simulations carried out for this study. These include the size of the simulation boxes (L_x × L_y), the number of lattice sites (N_s) and number of cells (N_c), the preferred cell area (A₀) and the values of surface tension (α). Waiting and sampling times (t_w and t_s), and the number of independent simulations for each value of α are presented in the final three rows, specified by α. The right-most column presents the values for the square-lattice simulation, by which we generated the fluid state from which the irregular lattice was obtained

	Square	Hexagonal	Irregular	Irr. gen.
L x × Ly	200^2	199.15 × 199.87	200.11 × 197.08	1280^2
N s	4 × 10^4	39 804	40 960	1280^2
N c	1000	1000	986	40 960
A 0	40	39.80	40	40
α	[1,4]	[1,4]	[1,4]	0.8
t w (MCS)	2 × 10^5, α ∈ [1,1.4]	2 × 10^5, α ∈ [1,1.66]	2 × 10^5–7 × 10^5, α ∈ [1,2.12]	2 × 10^4
	2 × 10^5, α ∈ [1.6,2.4]	2 × 10^5–1 × 10^6, α ∈ [1.68,2.4]	3 × 10^5–1 × 10^6, α ∈ [2.16,2.6]
	3 × 10^6∼ 8 × 10^6, α ∈ [2.6,4]	1.2 × 10^6, α ∈ [2.6,4]	1 × 10^6, α ∈ [2.8,4]
t s (MCS)	2 × 10^5–1 × 10^6, α ∈ [1,1.4]	2 × 10^5, α ∈ [1,1.66]	2 × 10^5–1 × 10^6, α ∈ [1,2.12]	—
	2 × 10^6–8 × 10^6, α ∈ [1.6,2.4]	2 × 10^6–6 × 10^6, α ∈ [1.68,2.4]	1 × 10^6–5 × 10^6, α ∈ [2.16,2.6]
	2 × 10^6, α ∈ [2.6,4]	2 × 10^6, α ∈ [2.6,4]	1 × 10^6, α ∈ [2.8,4]
Number of simulations per each α	50	20	20	1

---

Fig. 2 Creation and characterization of the irregular lattice. (a) Part of a snapshot of the simulation, which was used to create a fluid state of cells. Nearly 0.1% of all simulated cells are shown, outlined in black. The blue dots indicate the cell's centers of mass. (b) The Voronoi tessellation based on the snapshot in panel (a). Each polygon will serve as a lattice site for our irregular lattice. (c) Fast Fourier transform of the picture of the center of mass of all cells comprising the fluid-like state – a portion of which is shown in panel (a). We express the inverse wavelength in terms of a length unit u, which represents the average cell spacing.

Next, we applied Voronoi tessellation to these centers, imposing periodic boundary conditions, using the package Voro++,¹²⁹ see Fig. 2(b). Finally, we verified that the newly formed irregular lattice does not have any long-range orientational structure. Using the image analysis software ImageJ,¹³⁰ we performed a fast Fourier transform (FFT) on a snapshot of the center of mass of the lattice sites. The result is shown in Fig. 2(c), which reveals the uniform rings, that are indicative of structural homogeneity. Further simulation details will be provided in Section 2.3.

2.3 Simulation parameters

The dynamics of the system is generated through MC attempts, as described in Section 2.1. We refer to N_s MC attempts as a sweep and we measure ‘time’ in terms of MC sweeps (MCS). As the temperature here is merely a scaling of the cell membrane fluctuations, we set k_BT = 1 throughout. The number of MCSs used depends strongly on the state point under consideration since the system features slow dynamics. Table 1 provides the relevant choices for both the equilibration (or waiting) time t_w and the time over which we sampled t_s. We examined the evolution of total energy, the shape parameters, and the order parameter, both of which will be introduced shortly, to establish the appropriate value of t_w.

We wanted to evaluate the phase transition induced by the variation in the surface tension of the cells. Therefore, we assumed λ to be constant and equal to 1.0, while we changed α as the control parameter in the range of α ∈ [1.0,4.0] on all lattices to go from a disordered diffusive dynamics of the cells to an ordered arrested one. The value of α was changed in steps of 0.2 and smaller steps of (0.02–0.04) near the transition point, as appropriate. In all cases, we used periodic boundary conditions and modeled nearly 1000 cells, see Table 1 for the details. For the irregular lattice, we used a simulation box of size L_x = 200.11 and L_y = 197.08, such that it is perfectly tileable by hexagons having a target area of A₀ = 40. Every configuration studied was set up (and remained) confluent. We also performed several independent simulations to generate statistics.

As the initial condition on square lattice, we considered a rectangular arrangement of cells, each of which has dimensions 5 × 8. On the hexagonal and irregular lattices, we did the same for α ≤ 2.2, while taking an equilibrated snapshot (prepared at α = 2.2) for higher values of α. The reason for this was the slow equilibration of the system in the high-α range. Again, we emphasize that we started sampling after the system was equilibrated and all the transient effects were decayed.

2.4 Characterization of the results

We characterize the outcomes of our simulations in several ways. First, we compute the MSD, 〈r²(t)〉, by extracting the CMS of each cell after equilibrating the system, based on eqn (5).

〈r²(t)〉 = 〈(r(t + t_w) − r(t_w))²〉. (5)

Here, r is the position vector for the center-of-mass of each cell and t denotes ‘time’. Henceforth, we identify t as the number of MCS||. The angle brackets indicate averaging over the cells and taking an ensemble average. We also fitted power-laws to the extracted MSDs to evaluate the diffusive behavior of the cells. Here, we focused on the long-time behavior only, which we found to follow

〈r²(t)〉 ∝ t^β, (6)

where β is the scaling coefficient. In practice, we determined β from the slope of the MSD after taking the log of both the time and MSD. To calculate the effective diffusion coefficient D_eff, we use

〈r²(t)〉 = 4D_efft, (7)

for the long-time behavior, whenever β ≳ 0.95. By using D_eff, we identified solid-like and fluid-like states of the tissue. We did this by fitting a polynomial function to the diffusion coefficient, and pinpointing where the curvature of the function is maximally negative. Diffusion coefficient has been used in similar studies for the identification of phase transition in tissues.⁵³

Second, we consider the local bond-order parameters to establish the degree of hexatic order. This is computed as follows for the cell indexed k

In this equation, [scr N, script letter N](k) is the set of nearest neighbors to the k-th cell and N_n is the number of neighbors. Here, [scr N, script letter N](k) is considered simply six nearest neighbors** of the cell k. This means that N_n = 6 for all the cells. The angles θ(j,k) represent the counterclockwise angle between the x-axis and the vector connecting the center-of-mass of the cell k to that of j. The ι represents the complex identity, ι² = −1. For a perfect hexagonal arrangement, the absolute |ψ₆(k)| = 1, and the value of the hexatic order decreases with disorder. Typically, we average |ψ₆| over all cells and the production time.

Third, we complement the hexatic-order analysis by considering the dimensionless quantity called the shape index. For any 2D shape, this quantity is defined as , where P and A are the perimeter and the area, respectively. This parameter has been studied extensively in tissue mechanics, e.g., see ref. 18, 19, 23, 53, 106, 112, 133 and 134. A slightly modified version of this quantity is the isoperimetric quotient, which is defined as . It is dimensionless as well, and equal to 1.0 for the circle, while being close to 0 for highly elongated shapes.

Here, we should note that the cell perimeter derived from lattice-based models is not readily comparable to that in off-lattice models. This is because the borders of the cells are defined by the edges of the lattice sites, which enforces excess jaggedness to the cells – in mathematical language, a natural metric to distance calculations. To overcome this issue, we applied Voronoi tessellation to the center-of-mass (CMS) positions of the cells. This enabled us to study the effective areas and perimeters of the generated polygons and compute the associated isoperimetric quotient. Note that the effective cell shape obtained by applying Voronoi tessellation can be different from the original. There are alternative ways to smoothen the jagged borders of cells on a lattice, e.g., using elliptic Fourier analysis.^135–137 However, since one of the goals of this study is to compare the shape characteristics in CPM with those in vertex and Voronoi models, we decided to construct our effective polygons in the same manner. We use the subscript ‘V’, for example, q_V, to indicate the use of our Voronoi procedure. Fig. 3(a)–(c) shows a zoomed-in view of a snapshot that illustrates the process, as well as the obtained value of q_V. We should note that a correction was also introduced to evaluate the perimeter of the cells on the CPM.^123,126 However, we chose to use Voronoi tessellation to make a clear comparison between this study and the vertex and Voronoi models that are being used in this context of epithelia.

Fig. 3 Using Voronoi tessellation to determine the shape index of lattice-based cells. (a) A snapshot of cells on the disordered lattice. The colors indicate individual cells and the black circles their centers of mass. (b) Voronoi tessellated version of the cells based on the position of the centers of mass. (c) Voronoi tessellated cells are colored by their isoperimetric quotient q_V, as indicated by the color bar shown on the right-hand side. The data derives from a part of a simulation with α = 1.8. (d) Isoperimetric quotient of regular polygons, q_reg(n), as well as its generalization q_V(n*).

The isoperimetric quotient of regular polygons with n edges, q_reg(n), can be readily computed and reads

where A_reg(r,n) and P_reg(r,n) are the area and the perimeter of a regular polygon having n edges and a circumscribing circle with radius r. This function maps the number of the edges of regular polygons to their isoperimetric quotient, and vice versa. We can now straightforwardly apply the right-hand side of eqn (9) for non-integer values of n, which we call the generalized edge number and which we will identify using n*. Numerically inverting this function, we uniquely obtain n* for any given value of q_V. By doing so, we can assess the effective number of edges for a given (convex) cell shape. The dependence is shown in Fig. 3(d).

Lastly, we also considered the circularity C of our model cells. This parameter was introduced by Zunic and Hirota,^138,139 and is calculated as follows. Given the area A of a 2D connected object, and the moment of inertia tensor Ī, we write

where in the denominator, the trace of the tensor is taken. Since the trace is invariant under translation and rotation, C is independent of the coordinate system. One should be careful in connection to ref. 138 and 139, to note that people also use the term ‘circularity’ to refer to the quantity 4πA²/P, which is q² in our characterization. However, throughout this study, by circularity, we mean the quantity calculated using eqn (10). The reason we study this parameter is that, unlike the isoperimetric quotient, it is not directly dependent on the perimeter. As we have indicated above, there are issues in defining a perimeter length in lattice-based models. Circularity bypasses this issue and thus helps us learn how different ways of measuring cell roundness compare to each other. Circularity is in fact, one of the Hu moment invariants¹⁴⁰ that are widely used in image processing and pattern recognition.¹⁴¹ This quantity has also been used in studying cell morphology,¹⁴² branching patterns in organogenesis simulations,⁴¹ and pathological cell nucleus shape analysis.¹⁴³ However, to the best of our knowledge, thus far it has not been studied in the context of confluent tissue transitions. Similar to the isoperimetric quotient, we rely on a Voronoi tessellation to determine C_V. The circularity is much more strongly nonlinear than q in the number of edges of a regular polygon. Therefore, we do not invert it to establish an equivalent circularity-derived generalized edge number n°.

3 Results

In this section, we introduce the main results of our CPM simulations. As explained in Section 2.1, the control parameter in our simulations is α, which models the surface tension between cells. We start by showing that increasing α leads to a first-order phase transition from a disordered fluid to an ordered hexagonal phase on all lattices.

3.1 Hexatic order, artifacts, and diffusion

Fig. 4 shows the (average) hexatic orientational order parameter 〈|ψ₆|〉 as a function of α and several representative snapshots of a part of the simulation area. The average hexatic order is generally an increasing function of α, except for high values of α on the square lattice, which we will return to shortly. For all our lattices there is a substantial increase in 〈|ψ₆|〉 in the range 0.7 to 0.8, which appears to connect a disordered and an ordered ‘branch’ in the state space. This is indicative of the presence of a first-order phase transition from a disordered fluid to a hexagonal solid. The nature of the structural change is further confirmed by contrasting the first two columns of snapshots in Fig. 4(b) with each other. The centers of mass of the cells assume a hexagonal arrangement in the three subpanels.

Fig. 4 The formation of hexatic order on the three lattice types considered by increasing the surface tension. (a) The system-averaged hexatic order 〈|ψ₆|〉 as a function of α for the CPM on a square (blue), irregular (green), and hexagonal (red) lattice. The error bars show the standard error of the mean. The black dashed line indicates the transition value reported by ref. 115. (b) Snapshots of the simulations specified by their labels in panel (a); only a small part of the simulation box is shown (∼6%). The color of the cells indicates their local hexatic order.

In the case of a square lattice, the transition appears sharp. For the other two lattices, there is a smoother transition, which can be attributed to the fact that we work in the NVT ensemble. That is, when we average over several realizations of the system, these likely contain both fluid and hexagonal configurations, beyond the transition α. In support of this, we found patterns of coexistence between the disordered and ordered phases in some of the simulations, see Fig. S1 (ESI†) for an example snapshot. This provides further (tentative) evidence for the presence of a first-order transition. Note that unlike in fluid–solid phase coexistence for particle-based systems,^3,45–48 we have an area-constraint term in our Hamiltonians (1) and (3). This presumably narrows any density gap that could be present between the disordered and ordered phases, and makes it difficult to find strong evidence of coexistence. We will also not concern ourselves with the existence of a potential intermediate hexatic phase here, but referencing the literature,^115,144,145 it should be present.

Increasing α beyond the transition value, 〈|ψ₆|〉 ≈ 0.76, we find that the hexatic order increases further. Note that our transition value is comparable to that reported in ref. 115, providing additional confidence in our results. On a square lattice, 〈|ψ₆|〉 eventually assumes a maximum. We can appreciate the underlying cause of the high-α reduction by examining snapshot 7 in Fig. 4(b). This reveals that for the square lattice, the cells become distorted into a configuration with staircase-like borders. That is, the borders of cells in large parts of the simulation locally follow the lattice (are at 90° angles). This leads to the cell edges overall being directed at angles ±45° with respect to the horizontal (or vertical) and the cells assuming a rhombus-like shape. This, in turn, gives rise to the strong, unphysical peak in the distribution of hexatic order, as can be seen in Fig. S3 (ESI†). For the disordered lattice, there is no such a peak, and the distribution of hexatic order is smooth within the error, as can be appreciated from Fig. S3b (ESI†). Together these results underpin that on the disordered lattice, the artificial (rhomboid) cell shapes are not present.

The orientation of the cell borders on the hexagonal lattice is also a lattice artifact. Fig. 4(b), snapshots 5 and 8, show that the borders of neighbor cells at intermediate-to-high surface tension, follow specific directions. These preferential orientations are dictated by the hexagonal symmetry of the underlying lattice. The effect appears less ‘severe’ than for the cells on the square lattice, where artifacts lead to significant cell-shape distortions. However, we caution against drawing this conclusion, as the natural high-surface-tension shape of the cells is hexagonal, which masks the extent to which cells are constrained by the underlying lattice.

We further characterized the behavior of our systems by examining the MSD scaling exponent β, and the diffusion coefficient, D_eff, as defined in eqn (6) and (7), respectively. A plot of MSD on the disordered lattice is shown in Fig. S2 (ESI†). First, we evaluated β and D_eff by fitting lines to log–log plots of the late-time MSD. Fig. 5 shows both quantities for different lattices as a function of α (β is shown in the inset). For low values of α, we readily obtain diffusive scaling β = 1 within the error, and the system is in a fluid-like state. Since the definition (7) holds for β = 1, we only show D_eff for the cases where β departs from 1 by less than one standard error of the mean. Note that the decrease in D_eff can be in orders of magnitude over the entire range of considered α, though we do not consider this indicative of glassy behavior, as reported in other sources;^53,112,113 we will return to this point in Section 4. For high values of α we observed sub-diffusive behavior. It is likely that for some of these values, the dynamics eventually becomes diffusive when they are evaluated for a sufficiently long time. However, we did not test this, as these values of α are far enough away from the transition that they do not affect our analysis and conclusions.

Fig. 5 Particle diffusion indicates where the system transitions from fluid-like to solid-like behavior. The effective diffusion coefficient D_eff of cells as obtained from their long-time diffusive dynamics as a function of the interfacial energy α for the three different lattices as labeled. We only show data for which the scaling exponent of the MSD, β ≈ 1.0, and the (long-time) dynamics is indeed diffusive. The inset shows β as a function of α. The dashed lines in the main panel and the inset serve as guides to the eye. The star symbols localize the inflection points to fitted data, and the square symbols indicate the transition points which were determined from the second derivative of the fitted data. See the main text for the procedure.

Fig. 5 reveals that D_eff generally decreases with α. It is difficult to identify where the phase transition happens without performing an exhaustive analysis. Here, we therefore computed two properties. First, the decrease has an inflection point at a given value of α that is lattice-dependent. This inflection point is posited to be indicative of coexistence, i.e., it is caused by blending fluid-like and solid-like behavior in equal parts. Thus, we expect the inflection point to overestimate the value of α for which there is a transition. To locate the inflection α, we fitted polynomial curves to the semi-log plot of D_eff. The obtained points are indicated using stars in Fig. 5. We will refer to the associated values using and , respectively. Second, we identify the point where D_eff starts to drop rapidly by examining the second derivative of the fitted D_eff(α). To do so, we considered the point at which, the curvature of the fitted curve is maximally negative. This point is posited to be closer to the transition, i.e., at the end of the purely disordered branch. We denote this point using and which are showed by square symbol in Fig. 5.

3.2 Isoperimetric quotient and circularity

Fig. 6(a) shows the equilibrium distribution of the Voronoi-based isoperimetric quotient, q_V, for two different values of α, obtained by using an irregular lattice. Because our model tissue is confluent, the q_V for the majority of the cells assumes values in the square-to-hexagon range. Note that for the fluid-like state, there is a dip in the value q_V around that of a regular pentagon. This is a consequence of the definition of q_V, rather than a profound result. Voronoi cells with 5 vertices always have values of q_V ≤ q_reg(5) = 0.930, and cells with 6 vertices always have q_V ≤ q_reg(6) = 0.952. Whether the role of an arbitrary polygon in tiling the plain is similar to that of a pentagon or a hexagon, depends mostly on whether its q_V is closer to q_reg(5) or q_reg(6), rather than the number of its edges. It is clear that upon undergoing a phase transition, the distribution shifts directly from being double-peaked to being strongly peaked close to a hexagonal value – another indicator of a first-order phase transition. The fact that the model tissue has defects in combination with the properties of q_V, makes it that the peak is not exactly centered about q_V ≈ 0.952. The defects can be seen in Fig. 4(b), snapshots 6 and 9.

Fig. 6 The distribution of Voronoi-based isoperimetric quotient q_V for a liquid- and solid-like state. (a) For a disordered lattice, we computed the probability density function (PDF) of q_V for α = 1.8 (blue) and α = 2.2 (orange), which are in the fluid and crystal states, respectively. (b) The same data, but now converted to an effective number of edges n*. In both graphs, the error bars indicate the standard error of the mean, and values of q_V and n* that correspond to regular polygons are indicated using dotted vertical lines and the use of symbols/numbers.

The distribution of q_V is insightful, but to help understand the properties of the underlying system, it is beneficial to convert it to an effective number of edges n*, see Fig. 6(b) and the definition in Section 2.4. The advantage of working directly with n*, rather than q_V, is that it allows us to examine the partition of the distribution. We identify the nearest integer number to a given n* by ⌊n*⌉, which factors our range into distinct subsets of triangular, square, pentagonal, hexagonal, heptagonal, etc. neighborhoods. For example, all the cells having ⌊n*⌉ = 5 (i.e., 4.5 < n* ≤ 5.5) can be referred to as pseudo-pentagons. Integrating the PDF belonging to the pseudo-pentagons provides insight into the fraction of pentagonal cells in the system. This includes cells with 5 or more edges, but for which some of the edges are very small.

Fig. 7(a) shows fractions of pseudo-polygons f_n* present in the system as a function of α. Interestingly, the curves for f₅ and f₆ show a crossover that matches well where we locate the phase transition based on our analysis of the diffusion coefficient drop. The physical intuition in this representation is clear: to crystallize, the system must have a majority of pseudo-hexagons. Fig. 7(b) shows the relation between the crossover and the transition on different lattices. Here, we have plotted D_eff as a function of f₆/f₅. The transition value of diffusion coefficient, i.e., , is indicated by the horizontal dashed lines. In Fig. 7(b), it can be seen that when the ratio exceeds 1, the diffusion coefficient drops sharply, which is a precursor to full crystallization, on both the square and irregular lattices. However, on the hexagonal lattice, there is a slight mismatch between the crossover and the transition. The ESI† goes into a longer discussion of this crossover for regular lattices and, also considers the mean and median of the distribution. The latter quantity also appears to be a good quantifier of the transition for these systems, see Fig. S4 and S5 (ESI†).

Fig. 7 The trend in the fractions of pseudo-polygons can be used to determine the transition point. (a) The fractions of pseudo-polygons f_n* as a function of the surface tension α obtained for a disordered lattice. The number of edges for the pseudo-polygons is indicated in the legend and the colored area around the respective curves indicates the standard error of the mean. The vertical gray bar shows the range for which the system transitions from a fluid to a hexagonal solid; i.e., where the diffusion coefficient drops steeply. (b) The diffusion coefficient D_eff as a function of the ratio f₆/f₅. The horizontal dashed lines show for different lattices. The vertical dashed line highlights the crossover at f₆/f₅ = 1.

Returning to the data presented in Fig. 6, we see that the distributions for the fluid-like state (e.g., α = 1.8) have two peaks, while the solid-like state (e.g., α = 2.2) is characterized by a single peak. This statement holds in general and we can extract the positions of n* for the maximum; or the two maxima and the local minimum, see Fig. 8. Similar data for the square and hexagonal lattices are provided in Fig. S6 (ESI†). We note that the behavior of the peak bears the hallmarks of a first-order phase transition, with a jump in the value of the ‘effective order parameter’ at α ≈ 2.1, as we alluded to before. The minimum is effectively at n* = 5, as is a feature of the definition of q_V. When the system transitions into the solid state, the peak reaches a nearly constant value of n* ≈ 5.7, which is shared among the studied underlying lattices. This is commensurate with the model cells transitioning to a state that is geometrically similar between the square, irregular, and hexagonal lattices. However, we emphasize that this does not imply that the transition is unaffected by the properties of the underlying lattice, as we have provided evidence for above.

Fig. 8 The fluid-to-solid transition in the model cell system as characterized using the extrema in the distribution of the Voronoi-based isoperimetric quotients. The n* values for the local maxima are indicated in red and the local minimum in blue when it is present. The dashed lines are guides to the eye and the grey vertical bar indicates α interval in which we locate the transition, similar to Fig. 7(a). The dotted horizontal lines indicate the values for a regular pentagon and hexagon, respectively, as labeled. The error bars are smaller than the point size. All data was obtained for a disordered lattice underlying the CPM dynamics.

Lastly, we calculated the circularity of the Voronoi-generated polygons for our model tissues. Fig. 9(a) shows the distribution of C_V of the cells on the disordered lattice, for two values of α already considered in Fig. 6. It turns out that the distributions of C_V for all the values of α studied here are unimodal, unlike those of q_V, e.g., see Fig. 9 and Fig. S7 (ESI†). We also observe from our data that at the transition point, the mode of the distribution of C_V almost matches with the circularity of a regular pentagon, which we identify with C(5). If we apply the same strategy as we did to calculate n* from q_V, to calculate n° from C_V, the mode of the distribution lies around n° ≈ 5.1. This relation is shown in Fig. 9(b) for the disordered lattice, where we calculated D_eff for the long-time diffusive regime, as a function of the departure of the mode from the pentagonal value, i.e., mode (C_V) − C(5). Clearly, the value C(5) is an excellent estimator for the phase transition. Lastly, it should be reemphasized that circularity compares the area with moment of inertia. This gives it an advantage over the isoperimetric quotient, which depends on perimeter length.

Fig. 9 The circularity C_V of the Voronoi tessellation as an indicator of the phase transition. (a) The PDF of C_V for two different values of α on the disordered lattice. One is representative of the fluid-like regime (α = 1.8, blue), and the other, of the solid-like regime (α = 2.2, orange). The position of the peak of the distribution (green circles) is extracted from the red (polynomial) fits. The two dotted vertical lines and symbols indicate the value of circularity for a regular pentagon and hexagon, respectively. (b) The diffusion coefficient, D_eff as a function of the departure between the mode of the distribution, mode (C_V), and the value of the circularity for a regular pentagon C(5). The horizontal dashed line, and the gray box, show the average and standard error of . The data points are colored by the value of α, as indicated in the color bar to the right. In both panels, the error bars show the standard error of the mean.

4 Discussion

In this section, we provide additional context to our research. We have split this into five parts for convenience. The first concerns itself with the quality of our phase diagram and the lack of a glassy phase. The second makes the connection to experimental observations to further justify the use of the CPM. The third discusses resolving the lattice artifacts present on regular lattices. The fourth provides insight into the efficiency of our algorithm and the challenges in working with an irregular-lattice CPM. We close with a larger discussion on the use of shape parameters.

4.1 The phase transition

On three different lattices studied here, we observed an order–disorder transition. This transition was reported previously for the square-lattice version of CPM,¹¹⁵ where it was studied as a function of temperature. Here, we show the transition as a function of the surface tension and show that it is not limited to the classical square-lattice version of CPM. By the Mermin–Wagner theorem, there should be a narrow, intermediate band of hexatic phase in the phase diagram, as also reported in the literature.^115,144,145 However, this phase is notoriously difficult to identify and we did not carry out the detailed analysis required to establish its presence. Our focus is instead on the effect of lattice artifacts on the overall trends.

Bearing this caveat in mind, our transitions have the hallmarks of a first-order phase transition. These include a sharp change in the value of the hexatic order and the diffusion coefficient, and the presence of what appears to be coexistence between ordered and disordered cell populations, as shown in Fig. S1 (ESI†). One should be aware that the simulations are in NVT ensemble, thus, the presence of a coexistence region tends to smoothen averaged curves. To establish the coexistence behavior proper free-energy calculations would have to be performed, but this was not the main goal of our present work.

The value of the isoperimetric quotient at the transition, see Fig. S5 (ESI†), lies between that of a regular pentagon and a regular hexagon, which is commensurate with this observation. Additionally, as is plotted in Fig. 4(a), the transition value of hexatic order matches with the value reported in ref. 115, providing additional support for the similarity between our findings and theirs. Other studies in the literature that used ‘deformable’ particles, e.g., phase-field^60,63 and Fourier contour⁶⁴ approaches, have also reported the emergence of hexagonal solid phases in model tissues, lending further credence to our result.

This makes the transition in the CPM markedly different from the jamming transitions reported for the vertex^18,106 and active Voronoi models.^53,146 These appear to leave the structure amorphous and occur at an average (target) value of q = q_reg(5). The regular (perfect) pentagon is a shape that cannot tile the plane confluently. Thus, imposing this on all cells in the model tissue leads to frustration between local and global constraints, resulting in a disordered arrest. Considering the similarities between cellular-Potts, phase-field, and Fourier-contour models, we surmise that a lack of border flexibility may be partly responsible for the difference. That is, in the vertex and Voronoi descriptions, the cells are described as polygons whose dynamics is ruled by their vertices and centroids, respectively. This means that they have far fewer degrees of freedom to their dynamics when compared to CPM at an equal cell number.

In this context, the work by Sadhukhan and Nandi¹¹³ should be mentioned, who have reported glassy behavior in the CPM with a Hamiltonian restricting the perimeter of the cells, as well as their area. This choice brings their version of CPM closer to the vertex and Voronoi models. However, we believe that the observation of square cells by these authors may be attributed to lattice artifacts.³⁶ For their large preferred perimeter regime, the interlocking of cells can also be attributed to the underlying lattice. As a lattice-based model, CPM has the inherent weakness of dealing with target perimeter lengths. Our work most closely follows that of Chiang and Marenduzzo¹¹² who identify a glass transition for the square-lattice CPM. However, in reproducing their results for the case of no self-propulsion, we have conclusively shown that this glassy behavior is not present. We did this by studying the hexatic orientational order parameter across the transition. Lastly, we turn to the results of ref. 114. They report a similar value of the shape parameter at the transition. However, they identify this transition as a jamming one. This is likely a misinterpretation of their findings, and this issue could be resolved by examining the hexatic order in their ‘jammed’ phase.

4.2 Connection to experiment

There have been many studies focused on jamming transitions and glassy behavior in living tissues, e.g., reviewed by ref. 15 and 103. However, these are not the only arrested patterns that can exist in these systems. Hexagonally packed cell arrangements are also frequently found in epithelia across a range of species and developmental stages.^147–163 Although more complicated mechanisms such as adherens junctions¹⁶⁴ can play a role in maintaining cellular order in hexagonal cell packings, it is informative to distinguish models that can capture this type of packing by their nature. The importance of this study is that we showed that regular and irregular versions of CPM, show an order/disorder transition, above which they were able to readily achieve the crystallized configuration. This makes them suitable for studying tissues with hexagonal cell packing.

4.3 Resolving lattice artifacts

The main advantage we derive from our irregular underlying lattice is to rid our simulations of two types of (obvious) lattice artifacts, without having to resort to higher-order neighbor coupling, as has been discussed in, e.g., ref. 36. Fig. 4 and Fig. S3 (ESI†) both clearly provide evidence that a square lattice CPM with a Moore neighborhood rule can lead to an abnormal distribution of hexatic order at high surface tension. This happens because on the square lattice, the cells are essentially forced to assume the rhombus-like shape that is shown in Fig. 4(b), snapshot 7. There are lattice artifacts on the hexagonal lattice too, although they appear not to be as severe as they are on the square lattice. The cells become hexagonal, but these hexagonal cells are forced to align with the underlying lattice. The disordered lattice is free of both of these artifacts within the error. Firstly, the distribution of hexatic order for the simulations on the disordered lattice is completely smooth and free of any unphysical peak within error. Secondly, there is no obvious preferred direction for cell borders on the irregular lattice.

By resolving these two lattice artifacts, we decreased the influence of the lattice geometrical symmetries on the shape of cells. We should mention that artifacts are an inevitable feature of lattice-based models. The shape of cells is always to some extent dependent on the geometry of the underlying lattice. However, as we have shown, by breaking the orientational symmetries of the underlying lattice, one can significantly weaken this dependency. This is an important advantage because as discussed in several studies,^{19,146,165–167} the shape of the cells can be strongly connected to the dynamics of the tissue.

4.4 Efficiency and considerations

Our irregular-lattice version of CPM not only is free of demonstrable artifacts but it can also be applied to any lattice that derives from a Voronoi construction, including regular ones. This provides our description with greater flexibility. In fact, for simulations on the hexagonal lattice, we used the same code, but with a perfect hexagonal Voronoi lattice as an input. When using an irregular lattice, it may be prudent to generate several independent realizations of the lattice. This should further reduce the correlations when taking statistical averages; note that we used only one disordered lattice in the current study.

In terms of computational efficiency, an extra run time is incurred using our approach compared to the square-lattice CPM. This was, however, very reasonable, as the average simulation time per 10⁶ MC sweeps for the square-lattice CPM was generally in the range of 7 to 11 hours, while it was in the range of 8 to 15 hours for our generalized CPM. This timing data was obtained using cluster nodes equipped with, on average, 20 CPU cores and 64 GB of RAM. Even in the worst case, the run time was less than twice that of the regular CPM.

Note that to properly evaluate the perimeter of our cells, we used Voronoi tessellation on their centers of mass. Even on the irregular lattice, the contour length of the cells is artificially high, as we established by examining a very large circular cell in isolation. This is a limitation of the model, which could be overcome by using alternative perimeter evaluation algorithms such as Voronoi tessellation.

4.5 Shape-parameter-based characterization

Several studies have looked at distributions of properties of the cell neighborhood. Some considered the aspect ratio of the cells,^146,168 while others extracted the number of neighbors.^{49,161,169,170} Here, we calculated the distribution underlying the mean isoperimetric quotient, and from these, we extracted fractions of pseudo-n*-gons to gain a deeper insight into the system. Fig. 7 for the irregular and Fig. S4 (ESI†) for the regular lattices, respectively, show that the majority of the polygons in the subdiffusive state are pseudo-hexagons, while in the sufficiently fluid-like regime, they are pseudo-pentagons. Close to the transition, we have observed that their relative abundance in the sample crosses over.

We found a single study by Saito and Ishihara,⁶⁴ who examined isoperimetric-quotient distributions across the transition in their Fourier-contour model. They also reported a unimodal-to-bimodal transition in the distributions, as we have observed. The bimodality in the fluid-like regime is indicative of a coexistence between more-rounded and less-rounded cells, which is understood to lead to fluidity of the tissue. We note that in the current literature^{60,64,112–114,171} there is a tendency to focus on the mean value of the shape parameter as the indicator of transition. This is undoubtedly motivated by its relevance to the behavior of the vertex and Voronoi models.^{18,20,53–55,106,133,146,172–175} However, as we have argued above, these models fall into a different class. We, therefore, believe it to be informative in general to examine shape-parameter distributions rather than only means.

As a complementary quantity to describe the shape of the cells we studied the circularity parameter C,^138,139 which is in essence, one of Hu moment invariants.¹⁴⁰ This compares the area of shapes to their moment of inertia, rather than to their perimeter, which is the case for q. This makes C less sensitive to noise at the boundary or definition issues with perimeter length. That is, we could have performed the circularity analysis on the original cell shapes rather than their Voronoi-based counterparts. However, we found that applying Voronoi tessellation to the cells led to a better comparison between the two roundness measures.

5 Summary and outlook

In this study, we extended the regular (square-lattice) cellular Potts model to work on arbitrary (ir)regular lattices. The main motivation of this study was to eliminate lattice artifacts observed for square-lattice CPMs in the literature. We demonstrated that such artifacts could indeed be eliminated by using an irregular lattice generated from a fluid-like state using Voronoi tessellation. Our generalized CPM maintains many of the desirable features of a base CPM, at the price of a small computational overhead.

We gained the following insights using our algorithm on both irregular, square, and hexagonal lattices. First, there is an order–disorder transition in confluent cell monolayers when using CPM on all lattices that we considered. This sets the base CPM apart from the active Voronoi and vertex models, in which there is a rigidity transition that maintains the disorder of the fluid-like state for shape parameters close to those of a regular pentagon. This makes the CPM suited to describe epithelia in which ordered cell arrangements form, as these can be reached directly and straightforwardly from the disordered state. This, however, does not mean to imply that the CPM can never describe glassy dynamics. For example, glassy dynamics could be realized in the CPM by modifying the Hamiltonian or by introducing restrictions on rearrangements leading to crystalline order. Second, for the square lattice, the cell shape can be impacted by artifacts above the transition, while for the hexagonal lattice, the borders of cells are affected by the symmetry of the lattice. Our irregular-lattice CPM did not exhibit these undesirable features.

In addition, we gained deeper insight into the transition by studying the distributions of different quantities across the transition. These included the hexatic bond order, the isoperimetric quotient, and the circularity, which overlap in their descriptiveness of local neighborhoods. For the isoperimetric quotient, the transition is closely correlated to the crossover in the fraction of pentagonal and hexagonal neighborhoods found in our simulations. However, the transition does not lie at an average isoperimetric quotient of a pentagon, as is the case for the vertex and active Voronoi models. This is because, unlike these two systems, the shape is not a target parameter in the Hamiltonian describing the system. It is further important to realize that examining the mean isoperimetric quotient, or the related shape index, may give an incomplete picture of the behavior of the system. We exclude the athermal Voronoi model in this comparison, where the rigidity transition is known to be absent.¹⁰⁹

Looking forward, we note that artifacts can be straightforwardly eliminated through the introduction of an irregular grid. This may have additional advantages when coupling the dynamics of the cells to that of external fields, such as food or chemical fields in bacterial colonies.¹⁷⁶ We will explore these directions in modeling real biological tissues using our generalized CPM in the future and hope this will lead to wider adoption of our approach.

Author contributions

The authors contributed as follows to the study: conceptualization (JdG & HN), data curation (HN), formal analysis (HN), funding acquisition (JdG), investigation (HN), methodology (JdG & HN), project administration (JdG), resources (JdG), software (HN), supervision (JdG), validation (HN & JdG), visualization (HN), writing – original draft (HN), writing – review & editing (JdG).

Data availability

Open data package containing the means to reproduce the results of the simulations available at: [https://doi.org/10.24416/UU01-IBVNT1].

Conflicts of interest

The authors declare that there are no conflicts of interest.

Acknowledgements

The authors acknowledge financial support from the Netherlands Organization for Scientific Research (NWO) through Start-Up Grant 740.018.013. We are grateful to Bryan Verhoef, Kim William Torre, Meike Bos, and Davide Marenduzzo for useful discussions.

Notes and references

1. M. Merkel and M. L. Manning, Semin. Cell Dev. Biol., 2017, 67, 161–169.
2. S.-Z. Lin, S. Ye, G.-K. Xu, B. Li and X.-Q. Feng, Biophys. J., 2018, 115, 1826–1835.
3. M. Akiyama, T. Sushida, S. Ishida and H. Haga, Dev., Growth Differ., 2017, 59, 471–490.
4. P. Friedl and D. Gilmour, Nat. Rev. Mol. Cell Biol., 2009, 10, 445–457.
5. P. Rørth, Annu. Rev. Cell Dev. Biol., 2009, 25, 407–429.
6. M. P. De Leon, F.-L. Wen, G. J. Paylaga, Y.-T. Wang, H.-Y. Roan, C.-H. Wang, C.-D. Hsiao, K.-H. Lin and C.-H. Chen, Nat. Phys., 2023, 1–9.
7. H. Xu, Y. Huo, Q. Zhou, L. A. Wang, P. Cai, B. Doss, C. Huang and K. J. Hsia, Proc. Natl. Acad. Sci. U. S. A., 2023, 120, e2221040120.
8. A. Brugués, E. Anon, V. Conte, J. H. Veldhuis, M. Gupta, J. Colombelli, J. J. Muñoz, G. W. Brodland, B. Ladoux and X. Trepat, Nat. Phys., 2014, 10, 683–690.
9. M. Poujade, E. Grasland-Mongrain, A. Hertzog, J. Jouanneau, P. Chavrier, B. Ladoux, A. Buguin and P. Silberzan, Proc. Natl. Acad. Sci. U. S. A., 2007, 104, 15988–15993.
10. K. J. Cheung, E. Gabrielson, Z. Werb and A. J. Ewald, Cell, 2013, 155, 1639–1651.
11. J. Metzcar, Y. Wang, R. Heiland and P. Macklin, JCO Clin. Cancer Inform., 2019, 1–13.
12. O. Ilina, P. G. Gritsenko, S. Syga, J. Lippoldt, C. A. M. La Porta, O. Chepizhko, S. Grosser, M. Vullings, G.-J. Bakker, J. Starruß, P. Bult, S. Zapperi, J. A. Käs, A. Deutsch and P. Friedl, Nat. Cell Biol., 2020, 22, 1103–1115.
13. S. Grosser, J. Lippoldt, L. Oswald, M. Merkel, D. M. Sussman, F. Renner, P. Gottheil, E. W. Morawetz, T. Fuhs, X. Xie, S. Pawlizak, A. W. Fritsch, B. Wolf, L.-C. Horn, S. Briest, B. Aktas, M. L. Manning and J. A. Käs, Phys. Rev. X, 2021, 11, 011033.
14. L. Oswald, S. Grosser, D. M. Smith and J. A. Käs, J. Phys. D: Appl. Phys., 2017, 50, 483001.
15. E. Blauth, H. Kubitschke, P. Gottheil, S. Grosser and J. A. Käs, Front. Phys., 2021, 9, 666709.
16. X. Yin, Y.-Q. Liu, L.-Y. Zhang, D. Liang and G.-K. Xu, Nano Lett., 2024, 24, 3631–3637.
17. J.-M. Armengol-Collado, L. N. Carenza, J. Eckert, D. Krommydas and L. Giomi, Nat. Phys., 2023, 1–7.
18. D. Bi, J. H. Lopez, J. M. Schwarz and M. L. Manning, Nat. Phys., 2015, 11, 1074–1079.
19. J.-A. Park, J. H. Kim, D. Bi, J. A. Mitchel, N. T. Qazvini, K. Tantisira, C. Y. Park, M. McGill, S.-H. Kim, B. Gweon, J. Notbohm, R. Steward Jr, S. Burger, S. H. Randell, A. T. Kho, D. T. Tambe, C. Hardin, S. A. Shore, E. Israel, D. A. Weitz, D. J. Tschumperlin, E. P. Henske, S. T. Weiss, M. L. Manning, J. P. Butler, J. M. Drazen and J. J. Fredberg, Nat. Mater., 2015, 14, 1040–1048.
20. O. K. Damavandi, E. Lawson-Keister and M. L. Manning, bioRxiv, 2022,  DOI:10.1101/2022.06.01.494406.
21. A. Hopkins, M. Chiang, B. Loewe, D. Marenduzzo and M. C. Marchetti, Phys. Rev. Lett., 2022, 129, 148101.
22. S. M. Fielding, J. O. Cochran, J. Huang, D. Bi and M. C. Marchetti, Phys. Rev. E, 2023, 108, L042602.
23. M. J. Hertaeg, S. M. Fielding and D. Bi, Phys. Rev. X, 2024, 14, 011027.
24. A. Chamoli, A. S. Gosavi, U. P. Shirwadkar, K. V. Wangdale, S. K. Behera, N. K. Kurrey, K. Kalia and A. Mandoli, Oral Oncol., 2021, 121, 105451.
25. C. Levra Levron, L. Elettrico, C. Duval, G. Piacenti, V. Proserpio and G. Donati, Cell Death Differ., 2024, 1–12.
26. R. Okamoto and M. Watanabe, J. Gastroenterol., 2016, 51, 11–21.
27. M. J. Holtzman, D. E. Byers, J. Alexander-Brett and X. Wang, Nat. Rev. Immunol., 2014, 14, 686–698.
28. Retinal Pigment Epithelium in Health and Disease, ed. A. K. Klettner and S. Dithmar, Springer International Publishing, Cham, 2020, pp. 117–185.
29. I. T. Stancil, J. E. Michalski, D. Davis-Hall, H. W. Chu, J.-A. Park, C. M. Magin, I. V. Yang, B. J. Smith, E. Dobrinskikh and D. A. Schwartz, Nat. Commun., 2021, 12, 4566.
30. R. Alert and X. Trepat, Annu. Rev. Condens. Matter Phys., 2020, 11, 77–101.
31. B. A. Camley and W.-J. Rappel, J. Phys. D: Appl. Phys., 2017, 50, 113002.
32. M. R. Shaebani, A. Wysocki, R. G. Winkler, G. Gompper and H. Rieger, Nat. Rev. Phys., 2020, 2, 181–199.
33. F. Graner and J. A. Glazier, Phys. Rev. Lett., 1992, 69, 2013–2016.
34. M. Scianna and L. Preziosi, Cellular Potts models, CRC Press, Boca Raton, FL, 2013.
35. M. Scianna and L. Preziosi, Multiscale Model. Simul., 2012, 10, 342–382.
36. A. F. M. Marée, V. A. Grieneisen and P. Hogeweg, in The Cellular Potts Model and Biophysical Properties of Cells, Tissues and Morphogenesis, ed. A. R. A. Anderson, M. A. J. Chaplain and K. A. Rejniak, Birkhäuser Basel, Basel, 2007, pp. 107–136.
37. A. Szabó and R. M. H. Merks, Front. Oncol., 2013, 3, 87.
38. T. Hirashima, E. G. Rens and R. M. H. Merks, Dev., Growth Differ., 2017, 59, 329–339.
39. S. E. M. Boas, Y. Jiang, R. M. H. Merks, S. A. Prokopiou and E. G. Rens, Cellular Potts Model: Applications to Vasculogenesis and Angiogenesis, Springer International Publishing, 2018, pp. 279–310.
40. P. Hogeweg, J. Theor. Biol., 2000, 203, 317–333.
41. A. Poorkhanalikoudehi and K.-H. Zimmermann, WSEAS Trans. Biol. Biomed., 2021, 18, 170–182.
42. G. B. Ermentrout and L. Edelstein-Keshet, J. Theor. Biol., 1993, 160, 97–133.
43. H. Hatzikirou and A. Deutsch, Current Topics in Developmental Biology, Academic Press, 2008, vol. 81, pp. 401–434.
44. J. M. Nava-Sedeño, H. Hatzikirou, F. Peruani and A. Deutsch, J. Math. Biol., 2017, 75, 1075–1100.
45. Z. Kirchner, A. Geohagan and A. Truszkowska, Integr. Biol., 2024, 16, zyae005.
46. B. A. Camley, J. Zimmermann, H. Levine and W.-J. Rappel, PLoS Comput. Biol., 2016, 12, e1005008.
47. M. Basan, J. Elgeti, E. Hannezo, W.-J. Rappel and H. Levine, Proc. Natl. Acad. Sci. U. S. A., 2013, 110, 2452–2459.
48. M. Campo, S. K. Schnyder, J. J. Molina, T. Speck and R. Yamamoto, Soft Matter, 2019, 15, 4939–4946.
49. R. Farhadifar, J.-C. Röper, B. Aigouy, S. Eaton and F. Jülicher, Curr. Biol., 2007, 17, 2095–2104.
50. S. Alt, P. Ganguly and G. Salbreux, Philos. Trans. R. Soc., B, 2017, 372, 20150520.
51. A. G. Fletcher, M. Osterfield, R. E. Baker and S. Y. Shvartsman, Biophys. J., 2014, 106, 2291–2304.
52. K. Landsberg, R. Farhadifar, J. Ranft, D. Umetsu, T. Widmann, T. Bittig, A. Said, F. Jülicher and C. Dahmann, Curr. Biol., 2009, 19, 1950–1955.
53. D. Bi, X. Yang, M. C. Marchetti and M. L. Manning, Phys. Rev. X, 2016, 6, 021011.
54. F. Giavazzi, M. Paoluzzi, M. Macchi, D. Bi, G. Scita, M. L. Manning, R. Cerbino and M. C. Marchetti, Soft Matter, 2018, 14, 3471–3477.
55. X. Yang, D. Bi, M. Czajkowski, M. Merkel, M. L. Manning and M. C. Marchetti, Proc. Natl. Acad. Sci. U. S. A., 2017, 114, 12663–12668.
56. D. E. P. Pinto, M. M. T. da Gama and N. A. M. Araújo, Phys. Rev. Res., 2022, 4, 023186.
57. A. Moure and H. Gomez, Arch. Comput. Methods Eng., 2021, 28, 311–344.
58. G. Zhang and J. M. Yeomans, Phys. Rev. Lett., 2023, 130, 038202.
59. M. Nonomura, PLoS One, 2012, 7, e33501.
60. B. Loewe, M. Chiang, D. Marenduzzo and M. C. Marchetti, Phys. Rev. Lett., 2020, 125, 038003.
61. H. P. Jain, A. Voigt and L. Angheluta, Sci. Rep., 2023, 13, 10096.
62. D. Wenzel and A. Voigt, Phys. Rev. E, 2021, 104, 054410.
63. S. Monfared, G. Ravichandran, J. Andrade and A. Doostmohammadi, eLife, 2023, 12, e82435.
64. N. Saito and S. Ishihara, Active Deformable Cells Undergo Cell Shape Transition Associated with Percolation of Topological Defects, 2023.
65. J. Tkadlec, K. Kaveh, K. Chatterjee and M. A. Nowak, J. R.Soc., Interface, 2023, 20, 20230355.
66. H. Nemati, K. Kaveh and M. R. Ejtehadi, J. Theor. Biol., 2023, 564, 111436.
67. J. Renton and K. M. Page, J. Theor. Biol., 2021, 528, 110838.
68. J. Renton and K. M. Page, J. R.Soc., Interface, 2019, 16, 20180918.
69. V. S. K. Manem, K. Kaveh, M. Kohandel and S. Sivaloganathan, PLoS One, 2015, 10, e0140234.
70. A. Goriely, The Mathematics and Mechanics of Biological Growth, Springer New York, New York, NY, 2017, vol. 45, pp. 51–59.
71. C. Beatrici, C. Kirch, S. Henkes, F. Graner and L. Brunnet, Soft Matter, 2023, 19, 5583–5601.
72. J. M. Osborne, A. G. Fletcher, J. M. Pitt-Francis, P. K. Maini and D. J. Gavaghan, PLoS Comput. Biol., 2017, 13, e1005387.
73. J. A. Glazier and F. Graner, Phys. Rev. E: Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top., 1993, 47, 2128–2154.
74. J. C. M. Mombach, J. A. Glazier, R. C. Raphael and M. Zajac, Phys. Rev. Lett., 1995, 75, 2244–2247.
75. R. Z. Tan, J. Comput. Sci., 2024, 76, 102237.
76. G. Y. Ouaknin and P. Z. Bar-Yoseph, Biophys. J., 2009, 97, 1811–1821.
77. J. Käfer, P. Hogeweg and A. F. M. Marée, PLoS Comput. Biol., 2006, 2, e56.
78. L. V. Steijn, J. A. J. Wondergem, K. Schakenraad, D. Heinrich and R. M. H. Merks, Biophys. J., 2023, 122, 2791–2807.
79. L. Plazen, J. A. Rahbani, C. M. Brown and A. Khadra, Sci. Rep., 2023, 13, 4223.
80. I. M. N. Wortel, I. Niculescu, P. M. Kolijn, N. S. Gov, R. J. de Boer and J. Textor, Biophys. J., 2021, 120, 2609–2622.
81. I. Niculescu, J. Textor and R. J. de Boer, PLoS Comput. Biol., 2015, 11, e1004280.
82. K. Matsushita, H. Hashimura, H. Kuwayama and K. Fujimoto, J. Phys. Soc. Jpn., 2022, 91, 054802.
83. A. J. Kabla, J. R.Soc., Interface, 2012, 9, 3268–3278.
84. A. R. Noppe, A. P. Roberts, A. S. Yap, G. A. Gomez and Z. Neufeld, Integr. Biol., 2015, 7, 1253–1264.
85. M. Scianna, Comput. Biol. Med., 2015, 62, 33–54.
86. B. M. Rubenstein and L. J. Kaufman, Biophys. J., 2008, 95, 5661–5680.
87. E. Tsingos, B. H. Bakker, K. A. E. Keijzer, H. J. Hupkes and R. M. H. Merks, Biophys. J., 2023, 122, 2609–2622.
88. R. M. Crossley, S. Johnson, E. Tsingos, Z. Bell, M. Berardi, M. Botticelli, Q. J. S. Braat, J. Metzcar, M. Ruscone, Y. Yin and R. Shuttleworth, Front. Cell Dev. Biol., 2024, 12 DOI:10.3389/fcell.2024.1354132.
89. M. H. Swat, G. L. Thomas, J. M. Belmonte, A. Shirinifard, D. Hmeljak and J. A. Glazier, Methods in Cell Biology, Academic Press, 2012, vol. 110, pp. 325–366.
90. F. R. Cooper, R. E. Baker, M. O. Bernabeu, R. Bordas, L. Bowler, A. Bueno-Orovio, H. M. Byrne, V. Carapella, L. Cardone-Noott, J. Cooper, S. Dutta, B. D. Evans, A. G. Fletcher, J. A. Grogan, W. Guo, D. G. Harvey, M. Hendrix, D. Kay, J. Kursawe, P. K. Maini, B. McMillan, G. R. Mirams, J. M. Osborne, P. Pathmanathan, J. M. Pitt-Francis, M. Robinson, B. Rodriguez, R. J. Spiteri and D. J. Gavaghan, J. Open Source Software, 2020, 5, 1848.
91. M. Meier-Schellersheim and G. Mack, SIMMUNE, a tool for simulating and analyzing immune system behavior, arXiv, 1999, preprint, arXiv:cs/9903017 DOI:10.48550/arXiv.cs/9903017.
92. J. T. Daub and R. M. H. Merks, Vascular Morphogenesis: Methods and Protocols, Springer, New York, NY, 2015, pp. 67–127.
93. J. Starruß, W. de Back, L. Brusch and A. Deutsch, Bioinformatics, 2014, 30, 1331–1332.
94. I. M. Wortel and J. Textor, eLife, 2021, 10, e61288.
95. N. J. Savill and R. M. H. Merks, Single-Cell-Based Models in Biology and Medicine, Birkhäuser, Basel, 2007, pp. 137–150.
96. M. Durand and E. Guesnet, Comput. Phys. Commun., 2016, 208, 54–63.
97. M. Scianna and L. Preziosi, Comput. Biol. Med., 2016, 76, 94–112.
98. K. D. Nnetu, M. Knorr, J. Käs and M. Zink, New J. Phys., 2012, 14, 115012.
99. E.-M. Schötz, M. Lanio, J. A. Talbot and M. L. Manning, J. R.Soc., Interface, 2013, 10, 20130726.
100. L. Bocanegra-Moreno, A. Singh, E. Hannezo, M. Zagorski and A. Kicheva, Nat. Phys., 2023, 1–9.
101. E. Hannezo and C.-P. Heisenberg, Trends Cell Biol., 2022, 32, 433–444.
102. S. Garcia, E. Hannezo, J. Elgeti, J.-F. Joanny, P. Silberzan and N. S. Gov, Proc. Natl. Acad. Sci. U. S. A., 2015, 112, 15314–15319.
103. E. Lawson-Keister and M. L. Manning, Curr. Opin. Cell Biol., 2021, 72, 146–155.
104. T. E. Angelini, E. Hannezo, X. Trepat, M. Marquez, J. J. Fredberg and D. A. Weitz, Proc. Natl. Acad. Sci. U. S. A., 2011, 108, 4714–4719.
105. D. Bi, J. H. Lopez, J. M. Schwarz and M. L. Manning, Soft Matter, 2014, 10, 1885–1890.
106. E. C. Thomas and S. Hopyan, Biophys. J., 2023, 122, 4264–4273.
107. S. Pandey, S. Kolya, S. Sadhukhan and S. K. Nandi, The unusual glassy dynamics of confluent epithelial monolayer is nearly ideal for mode-coupling theory, arXiv, 2023, preprint, arXiv:2306.07250 [cond-mat, physics:physics] DOI:10.48550/arXiv.2306.07250.
108. D. M. Sussman, M. Paoluzzi, M. C. Marchetti and M. L. Manning, EPL, 2018, 121, 36001.
109. D. M. Sussman and M. Merkel, Soft Matter, 2018, 14, 3397–3403.
110. P.-F. Lenne and V. Trivedi, Nat. Commun., 2022, 13, 664.
111. N. I. Petridou, B. Corominas-Murtra, C.-P. Heisenberg and E. Hannezo, Cell, 2021, 184, 1914–1928.
112. M. Chiang and D. Marenduzzo, EPL, 2016, 116, 28009.
113. S. Sadhukhan and S. K. Nandi, Phys. Rev. E, 2021, 103, 062403.
114. A. J. Devanny, D. J. Lee, L. Kampman and L. J. Kaufman, bioRxiv, preprint, 2023,  DOI:10.1101/2023.07.10.548321.
115. M. Durand and J. Heu, Phys. Rev. Lett., 2019, 123, 188001.
116. R. B. Potts, Math. Proc. Cambridge Philos. Soc., 1952, 48, 106–109.
117. Color Atlas of Cytology, Histology, and Microscopic Anatomy, ed. W. Kühnel, Georg Thieme Verlag KG, 4th edn, 2003, p. 102.
118. M. S. Steinberg, Curr. Opin. Genet. Dev., 2007, 17, 281–286.
119. R. Belousov, S. Savino, P. Moghe, T. Hiiragi, L. Rondoni and A. Erzberger, Poissonian cellular Potts models reveal nonequilibrium kinetics of cell sorting, arXiv, 2024, preprint, arXiv:2306.04443 [cond-mat, physics:physics],  DOI:10.48550/arXiv.2306.04443.
120. J. L. Antonakos and K. C. Mansfield, Practical data structures using C/C++, 1999, pp. 165–190.
121. R. Allena, M. Scianna and L. Preziosi, Math. Biosci., 2016, 275, 57–70.
122. N. Guisoni, K. I. Mazzitello and L. Diambra, Front. Phys., 2018, 6, 61.
123. E. G. Rens and L. Edelstein-Keshet, PLoS Comput. Biol., 2019, 15, e1007459.
124. E. Leoncini, PhD thesis, INRIA, 2010.
125. P. Van Liedekerke, M. M. Palm, N. Jagiella and D. Drasdo, Comput. Part. Mech., 2015, 2, 401–444.
126. R. Magno, V. A. Grieneisen and A. F. Marée, BMC Biophys., 2015, 8, 8.
127. H. X. Zhu, S. M. Thorpe and A. H. Windle, Philos. Mag. A, 2001, 81, 2765–2783.
128. S. Kumar and S. K. Kurtz, Mater. Charact., 1993, 31, 55–68.
129. C. H. Rycroft, Chaos, 2009, 19, 041111.
130. C. A. Schneider, W. S. Rasband and K. W. Eliceiri, Nat. Methods, 2012, 9, 671–675.
131. P. Digregorio, D. Levis, A. Suma, L. F. Cugliandolo, G. Gonnella and I. Pagonabarraga, Phys. Rev. Lett., 2018, 121, 098003.
132. J. R. Borba, C. Brito, P. Migowski, T. B. Vale, D. A. Stariolo, S. R. Teixeira and A. F. Feil, J. Phys. Chem. C, 2013, 117, 246–251.
133. M. Krajnc, Soft Matter, 2020, 16, 3209–3215.
134. J.-j Li, Y.-l He and B.-q Ai, Soft Matter, 2023, 19, 3849–3858.
135. F. P. Kuhl and C. R. Giardina, Comput. Graphics Image Process., 1982, 18, 236–258.
136. L. Tweedy, B. Meier, J. Stephan, D. Heinrich and R. G. Endres, Sci. Rep., 2013, 3, 2606.
137. Y. E. Sánchez-Corrales, M. Hartley, J. van Rooij, A. F. Marée and V. A. Grieneisen, Development, 2018, 145, dev156778.
138. J. Žunić and K. Hirota, Measuring Shape Circularity, Springer, Berlin, Heidelberg, 2008, pp. 94–101.
139. J. Žunić, K. Hirota and P. L. Rosin, Pattern Recognit., 2010, 43, 47–57.
140. M.-K. Hu, IEEE Trans. Inf. Theory, 1962, 8, 179–187.
141. M. W. Nasrudin, N. S. Yaakob, N. A. Abdul Rahim, M. Z. Zahir Ahmad, N. Ramli and M. S. Aziz Rashid, J. Phys.: Conf. Ser., 2021, 1962, 012028.
142. L. Kopanja, Z. Kovacevic, M. Tadic, M. C. Žužek, M. Vrecl and R. Frangež, Histochem. Cell Biol., 2018, 150, 93–102.
143. M. Ammar, K. Baiche, S. Mahmoudi and M. Benazzouz, 2017 5th International Conference on Electrical Engineering – Boumerdes (ICEE-B), 2017.
144. A. Pasupalak, L. Yan-Wei, R. Ni and M. Pica Ciamarra, Soft Matter, 2020, 16, 3914–3920.
145. J.-j Li and B.-q Ai, New J. Phys., 2021, 23, 083044.
146. L. Atia, D. Bi, Y. Sharma, J. A. Mitchel, B. Gweon, S. A. Koehler, S. J. DeCamp, B. Lan, J. H. Kim, R. Hirsch, A. F. Pegoraro, K. H. Lee, J. R. Starr, D. A. Weitz, A. C. Martin, J.-A. Park, J. P. Butler and J. J. Fredberg, Nat. Phys., 2018, 14, 613–620.
147. D. J. Cislo, F. Yang, H. Qin, A. Pavlopoulos, M. J. Bowick and S. J. Streichan, Nat. Phys., 2023, 19, 1201–1210.
148. D. A. Sun and N. H. Patel, Wiley Interdiscip. Rev.: Dev. Biol., 2019, 8, e355.
149. A.-K. Classen, K. I. Anderson, E. Marois and S. Eaton, Dev. Cell, 2005, 9, 805–817.
150. F. Pilot and T. Lecuit, Dev. Dyn., 2005, 232, 685–694.
151. B. Aigouy, R. Farhadifar, D. B. Staple, A. Sagner, J.-C. Röper, F. Jülicher and S. Eaton, Cell, 2010, 142, 773–786.
152. T. Hayashi and R. W. Carthew, Nature, 2004, 431, 647–652.
153. H. Togashi, S. R. Davis and M. Sato, Dev. Biol., 2024, 506, 1–6.
154. T. Lecuit and P.-F. Lenne, Nat. Rev. Mol. Cell Biol., 2007, 8, 633–644.
155. M. C. Gibson, A. B. Patel, R. Nagpal and N. Perrimon, Nature, 2006, 442, 1038–1041.
156. D. Ortolan, R. Sharma, A. Volkov, A. Maminishkis, N. A. Hotaling, L. A. Huryn, C. Cukras, S. Di Marco, S. Bisti and K. Bharti, Proc. Natl. Acad. Sci. U. S. A., 2022, 119, e2117553119.
157. D. Napoli, M. Biagioni, F. Billeri, B. Di Marco, N. Orsini, E. Novelli and E. Strettoi, Int. J. Mol. Sci., 2021, 22, 5381.
158. S. T. Islam, C. Cheng, J. Parreno and V. M. Fowler, Invest. Ophthalmol. Visual Sci., 2023, 64, 20.
159. A. Jain, V. Ulman, A. Mukherjee, M. Prakash, M. B. Cuenca, L. G. Pimpale, S. Münster, R. Haase, K. A. Panfilio, F. Jug, S. W. Grill, P. Tomancak and A. Pavlopoulos, Nat. Commun., 2020, 11, 5604.
160. S. B. Lemke and C. M. Nelson, Curr. Biol., 2021, 31, R1098–R1110.
161. J. A. Zallen and R. Zallen, J. Phys.: Condens. Matter, 2004, 16, S5073–S5080.
162. A. Hočevar and P. Ziherl, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2009, 80, 011904.
163. R. I. Johnson, Dev. Biol., 2021, 478, 173–182.
164. O. Campàs, I. Noordstra and A. S. Yap, Nat. Rev. Mol. Cell Biol., 2024, 25, 252–269.
165. M. Luciano, M. Versaevel, E. Vercruysse, A. Procès, Y. Kalukula, A. Remson, A. Deridoux and S. Gabriele, Biophys. Rev., 2022, 3, 011305.
166. P. Gottheil, J. Lippoldt, S. Grosser, F. Renner, M. Saibah, D. Tschodu, A.-K. Poßögel, A.-S. Wegscheider, B. Ulm, K. Friedrichs, C. Lindner, C. Engel, M. Löffler, B. Wolf, M. Höckel, B. Aktas, H. Kubitschke, A. Niendorf and J. A. Käs, Phys. Rev. X, 2023, 13, 031003.
167. R. Alert and X. Trepat, Phys. Today, 2021, 74, 30–36.
168. S. Sadhukhan and S. K. Nandi, eLife, 2022, 11, e76406.
169. D. S. Roshal, K. Azzag, K. K. Fedorenko, S. B. Rochal and S. Baghdiguian, Phys. Rev. E, 2023, 108, 024404.
170. D. Sánchez-Gutiérrez, M. Tozluoglu, J. D. Barry, A. Pascual, Y. Mao and L. M. Escudero, EMBO J., 2015, 35, 77–88.
171. G. Zhang, R. Mueller, A. Doostmohammadi and J. M. Yeomans, J. R.Soc., Interface, 2020, 17, 20200312.
172. D. L. Barton, S. Henkes, C. J. Weijer and R. Sknepnek, PLoS Comput. Biol., 2017, 13, e1005569.
173. L. Yan and D. Bi, Phys. Rev. X, 2019, 9, 011029.
174. J. Huang, J. O. Cochran, S. M. Fielding, M. C. Marchetti and D. Bi, Phys. Rev. Lett., 2022, 128, 178001.
175. X. Li, A. Das and D. Bi, Proc. Natl. Acad. Sci. U. S. A., 2018, 115, 6650–6655.
176. O. Hallatschek, P. Hersen, S. Ramanathan and D. R. Nelson, Proc. Natl. Acad. Sci. U. S. A., 2007, 104, 19926–19930.

---

Footnotes

† Electronic supplementary information (ESI) available. See DOI: https://doi.org/10.1039/d4sm00445k

‡ It should be mentioned that, in the context of condensed matter, the glass and jamming transition are different and have distinct underlying physics. However, here, we chose to include all the nomenclature that is used in the literature of tissue mechanics without judging the accuracy of the specific use.

§ This neighborhood [scr N, script letter N]_c(i) is not necessarily the same as the one used for the computation of the Hamiltonian [scr N, script letter N](i) and may be defined separately.

¶ For this study, we did not need to introduce any specific length unit. Lengths may be subsumed in the definitions of the prefactors. Only for the fast Fourier transforms, we introduce the length scale u, which makes the wave space vectors scale as u⁻¹.

|| This definition is subject to the caveat that MCS can be understood to be proportional to the real time in a system, but they are not the actual time.

** Other definitions of the local hexatic bond order are possible,^131,132 but we considered the definition provided in eqn (8) the most appropriate for our purposes.

---