Signatures of Jamming in the Cellular Potts Model

We explore the jamming transition in the Cellular Potts Model (CPM) as a function of confinement, cell adhesion, and cell shape. To accurately characterize jamming, we compare Potts simulations of unconfined single cells, cellular aggregates, and confluent monolayers as a function of cell adhesion energies and target cell shape. We consider metrics that may identify signatures of the jamming transition, including diffusion coefficients, anomalous diffusion exponents, cell shape, cell-cell rearrangements, and velocity correlations. We find that the onset of jamming coincides with an abrupt drop in cell mobility, rapid transition to sub-diffusive behavior, and cessation of rearrangements between neighboring cells that is unique to confluent monolayers. Velocity correlations reveal collective migration as a natural consequence of high energy barriers to neighbor rearrangements for certain cell types. Cell shapes across the jamming transition in the Potts model are found to be generally consistent with predictions of vertex-type simulations and trends from experiment. Finally, we demonstrate that changes in cell shape can fluidize cellular monolayers at cellular interaction energies where jamming otherwise occurs.


Introduction
Coordinated cell migration is critical to a variety of biological functions including cell migration, morphogenesis, and cancer invasion 1,2 . For example, during morphogenesis tissues experience dramatic changes in structure as cells rearrange, moving in collective swaths as complex tissues form and take shape 3,4 . For these processes to occur, cells must be able to shift between relatively immobile and highly mobile states. The ability of cells to move in a coordinated fashion is necessary for the shaping and reshaping of tissue, while the ability to remain in place is essential to the long-term stability of a tissue. Complex signaling pathways underlie these cell motility programs 5 .
Broadly, a transition between high and low mobility states may be described as a jamming transition, a switch between a liquid-like and disordered solid-like state in which the dynamics of the system slow dramatically. Jamming was initially described for a system of frictionless, deformable spheres where a transition from free diffusion to collective motion preceded the abrupt arrest of the system components as a function of decreasing temperature, increasing density, or decreasing shear stress [6][7][8] . Active systems, including epithelial cells, have been shown to undergo analogous transitions from diffusive behavior to collective motility that precedes arrest to a solidlike state [9][10][11] . In recent years, the jamming transition has been extensively investigated in such active systems, with particular focus on understanding its implications in biological processes 12,13 .
In biological systems, the jamming transition has been understood to describe a precipitous dynamical slowdown and/or arrest of motion of the system components and has been investigated in a number of experimental and theoretical contexts. The transition has been shown to have clinical relevance in asthma, where primary cells from asthmatic donors were characterized by their resistance to jamming relative to non-asthmatic donor cells 14,15 . Unjamming has been implicated in wound healing and cancer progression 16,17 . Collective motility programs observed in wound healing are inherently connected with glassy behavior in cells 18 . In cancer, epithelial cells undergo the epithelial to mesenchymal transition (EMT) in which cell contacts are disrupted and cells become migratory 19,20 . In 2D and 3D cell culture, cells that have progressed through EMT display fluid-like behavior while healthy epithelial cells remain stationary, as in a solid [21][22][23][24] .
However, it remains unclear to what extent unjamming and EMT are distinct, and to what extent one or both are required in cancer invasion 20,25,26 .
For cellular systems, a hypothetical jamming phase diagram has been proposed, with jamming occurring as a function of increasing cell density, increasing cell-cell adhesion, and decreasing myosin dependent-motile forces 27 . Effective cell adhesion is proposed to be governed by an interplay of adhesion and cortical tension that is critical to determining the energy barriers to cell movement that dictate fluid or solid-like behavior 28,29 . Here, cadherins primarily dictate the degree of cell-cell adhesion while actomyosin contractility dictates the strength of cortical tension opposing the extension of cell-cell contacts 30,31 .
The vertex model, and extensions thereof, has proven to be a popular method and powerful tool for probing the jamming transition in biological systems [32][33][34] . In vertex model simulations, cells are composed of a series of vertices and edges, which evolve according to an energy functional. These simulations shed light on the forces at play in epithelial sheets that lead to the cellular dynamics observed during migratory processes as well as jamming. Vertex-type simulations have specifically highlighted the importance of cell shape in the jamming transition 35,36 . In particular, it has been shown that the jamming transition occurs at a characteristic value of cell shape index, the dimensionless ratio of perimeter to area, P/A 1/2 , 34 and the distribution of cell shape indices correlates with liquid-or solid-like behavior of cells 37,38 . Experimental work has further confirmed the value of cell shape index as a structural descriptor of jamming within a collection of cells -one that is both easily accessible and universal across biological systems 37 .
The Cellular Potts Model (CPM), a type of lattice-based Monte Carlo simulation, has been a staple of cell migration investigations since its introduction but has only rarely been used to investigate jamming [39][40][41][42][43] . The advantages of this model are its relative simplicity and access to sub-cellular detail. Despite its simplicity, the CPM has seen use in diverse contexts such as avascular tumor growth and cell migration, including as a tool to study both cell-matrix interactions and forces exerted during migration [44][45][46][47][48] . The CPM is an extension of the Q-Potts model, in which cells are modeled as an assemblage of contiguous like-spins on a lattice. Adjacent unlike-spins are assigned energies of adhesion, and the CPM stochastically evolves to minimize the energy of the system according to a Hamiltonian typically composed of area, perimeter, and adhesion energy terms. Cell motility arises from spin flips during each Monte Carlo step (MCS) that are accepted with a conditional probability associated with the change in total system energy.
Most CPM implementations use similar forms of the Potts Hamiltonian, though variations are sometimes seen in which the area or perimeter term is not included 48,49 or in which additional terms are added to capture more complex cell behaviors 40,50,51 . CPM studies of cell migration have recapitulated aspects of experimental epithelial collective cell migration and, recently, the CPM was used to demonstrate glassy dynamics of cellular monolayers reminiscent of the jamming transition 42,43 . Taken together, these works suggest that the CPM could also be a useful tool to investigate cellular jamming. Towards this end, it is necessary to perform a more rigorous analysis of the model in the context of indicators of jamming identified by and used in vertex-type models and experimental studies, as well as make connections between collective migration and jamming in Potts model studies.
In this work, we employ the CPM to investigate the jamming transition by analyzing the motility of cells under varying degrees of confinement, from unconfined cells, to aggregates, to confluent monolayers. We find cellular monolayers display a distinctive abrupt drop in diffusivity as a function of cell adhesion, with jamming marked by a transition to sub-diffusivity. Cell shape indices in the CPM are found to be similar to and generally consistent with predictions of jamming in other model and experimental systems. Within the CPM, we propose an approach for characterizing cellular rearrangements as a metric of jamming analogous to explicitly defined cell rearrangements in the vertex model. Velocity correlations in cell monolayers show increasing spatial correlations indicative of collective motion as cell adhesion energies increase, with a maximum in this collective motion directly preceding cell arrest. Finally, we demonstrate that increasing cell shape index in dense monolayers is sufficient to fluidize an otherwise jammed monolayer, and that a combination of mobile anisotropic cells and otherwise arrested rounder cells can fluidize a monolayer independently of cell-cell adhesion energy, demonstrating the potential impact of EMT-induced shape changes on cell motility.

Model Cellular Potts Model
To investigate jamming, we utilized a two-dimensional CPM. It consists of cells, existing on a square lattice of dimensions x (L = 100), with periodic boundaries, each of which has a unique spin integer, ( , ) ∈ {1, 2, … , }, where ( , ) identifies the position of a lattice site. The environment is represented by the spin of ( , ) = 0. Sub-domains with identical spin ≠ 0 represent a distinct cell, and they can be divided into cells of a particular type ( ). The simulation runs by Monte Carlo methods, using a modified Metropolis algorithm 52 , with iterative sequences of spin flips that attempt to minimize the total energy of the system. Each time step, or Monte Carlo sweep (MCS), is composed of flips. For each flip, the spin of a randomly selected lattice site ( , ) is changed to the spin of one of its neighboring sites ( ′, ′) of different spin integer ( ( , ) ≠ ( , )). A flip is accepted with probability given by: where is the Hamiltonian of the system. The Boltzmann weighted probability function guides the simulation so that there is a stochastic reduction in the total energy of the system. The Hamiltonian at any time, , can be divided into three components: The adhesion component is based on the differential adhesion hypothesis 53 Table S1.
For simulation of isolated cells, individual cells were deposited randomly on an empty lattice. Simulations were run for 5 MCS to allow cells to equilibrate to reasonable shapes after deposition. For simulation of cell aggregates, cell centers were randomly deposited in a circular region in the center of the lattice. A Voronoi tessellation was generated from the centroids to divide the region into cells. The resulting tessellation was discretized to a 100 x 100 lattice for simulations. As in isolated cells, simulations were run for a short time (100 MCS, in this case) to allow cells to equilibrate to reasonable shapes after deposition. Equilibration times are short but sufficient to attain reasonable cell shapes for each simulation in the cases described above. After this initial equilibration, simulations were run for 1000 MCS or 100,000 MCS, for isolated cells and aggregates, respectively. For monolayer simulations, cell centers were randomly deposited across the full 100 x 100 lattice, and a Voronoi tessellation was used to construct cell shapes from these centroids, as in aggregates. Monolayers were similarly initially run for normalization of cell shape, here for a longer period of 1000 MCS due to the more crowded system, then run for a total of 100,000 MCS. Isolated cell simulations were run such that there were 20 separate individual cell trajectories for analysis. Aggregate and monolayer simulations were run in triplicate (150 and 300 total cells, respectively) for each parameter set.

Cell Adhesion and Sub-Diffusion in the Potts Model
Cell trajectories were analyzed in greater detail by examining mean squared displacements (MSDs). MSD was calculated via 〈MSD〉 = 〈r(t) − r(0)〉, where brackets denote averaging over cells but not time lags, as such temporal averaging may mask sub-diffusion 56 . Without temporal averaging, ensemble averaged MSDs can suffer from noise. To improve the confidence of subsequent fits, we utilized the mean maximal excursion (MME) method, described in more detail in Ref. 57 . MME curves report similar dynamics and give better quality fits than their MSD counterparts (Fig. S1). For extraction of the anomalous diffusion exponent, MME curves were fit   S2). As such, J acts as an effective temperature, though it is important to note that the CPM separately includes a temperature parameter that, when decreased, also results in decreased acceptance of spin flips. The temperature parameter is fixed at 1 in these simulations (Table S1) (Fig. S3).
Distributions of diffusion coefficients broaden with increasing Jcell in monolayers, a hallmark of glassy systems indicating heterogeneous dynamics 59 . Taken together, it is clear that crowding effects from monolayer confinement are critical to provoking sub-diffusive and frustrated, heterogeneous cell dynamics associated with glassy behavior and jamming.

Cell Shapes in the Potts Model
Cell shapes are widely used as an indicator of jamming in a variety of experimental and simulated systems. Cell shape is typically defined as a dimensionless perimeter to area ratio (cell shape index = / / ), which is a readily accessible quantity in both simulations and experimental work. The jamming transition is expected to occur at a shape index of 3.81, with values at or below 3.81 associated with jammed systems 34 . We considered cell shape index to assess whether this commonly used metric is associated with jamming in the CPM. To minimize grid effects in the CPM, cell perimeters were approximated using a previously described method (See Supporting Text 1 for more detail) 46,47 , and area was determined by counting the number of lattice sites comprising a given cell. For all simulations discussed so far, the target shape index was set to that of a circle (3.55), as shown in Table S1.
We first analyzed cell shapes by examining the fluctuations in shape over time (Fig. S4).
In all cases, cell shape does not evolve appreciably after the early stages of the simulation; thus, we consider the distribution of cell shapes at the final timestep to be representative of the simulation as a whole. For aggregates, cell shape decreases with increasing Jcell but does not fully mirror the trends seen in Fig. 2b. Shape indices for simulations with Jcell = 0 -0.75 lie above the expected jamming threshold, while shape indices for Jcell ≥ 1 fall below this cutoff index (Fig. 3a).
In cases where aggregates disaggregate (Jcell ≥ 1.75, Fig. 3a), we do not expect cell shape to remain a useful structural metric of jamming, since the simulated cells are no longer in contact.
Nevertheless, we show the complete set of shape indices for each simulated set. For aggregates, we also considered the possibility that the total number of cells might alter cell shape index, due to differences in cellular environment between interior and exterior cells, specifically in cases where Jcell < Jenv and aggregates remain intact. However, we see no notable shift in cell shape indices as a function of aggregate size (Fig. S5), indicating negligible dependence on this variable.
In monolayers, the shape index drops below the jamming threshold at Jcell = 1 (Fig. 3b), which closely corresponds with the observed transition to sub-diffusion, which occurs at Jcell = 1.25 (Fig.   2c). In addition, a narrowing of shape distributions at Jcell ≥ 1.25 is seen, consistent with the expectation that in deeply jammed states cell shapes become more uniform 37 .  are exchanging neighbors beyond single lattice site contact fluctuations. Indeed, we find this behavior to hold for both aggregate and monolayer simulations (Fig. 4). We first interrogated the number of neighbor exchanges across time intervals spanning 100 -3x10 4 MCS in CPM aggregates (Fig. 4a). Here, we see the appearance of a crossover point, where the number of cells undergoing only single exchanges is exceeded by cells undergoing multiple changes of neighbors.

Cellular Rearrangements in the Potts Model
Since the appearance of this crossover point occurs at the longest timescales, we next investigated Importantly, the crossover metric described here accurately captures sub-diffusivity trends and does not misidentify contact fluctuations as functional cell motility.

Effects of Cell Shape on Jamming in the CPM
We next varied target cell shape to probe the effects of shape on the onset of jamming within the CPM as well as to further explore the ability of the shape index to predict jamming in the CPM. As discussed above, the target area and target perimeter of cells were previously set to a circular target shape index of approximately 3.55. In experimental contexts, more diverse cell morphologies and, in particular, highly elongated cells are regularly observed, especially among aggressive cancer cells. Therefore, we altered target perimeter to probe the effect of shape variation on the fluid-like or solid-like behavior of cells in these simulations. Given the clearer signs of jamming in cellular monolayers relative to unconfined simulations, we restricted this analysis to monolayers where unjamming would be more challenging.
We hypothesized that increasing the target cell perimeter, while keeping all other variables constant, would increase cell shape index as well as cell membrane fluctuations, allowing cells to remain mobile and pushing the jamming transition to higher values of Jcell, the adhesion parameter.
Thus, we performed a set of monolayer simulations in which the target perimeter of cells was increased such that target shape index increased by a factor of ≈ 2 (Table S1), corresponding to more elongated cells. To evaluate the onset of jamming, again the diffusion coefficient, anomalous diffusion exponent, shape index, and neighbor exchanges were considered (Fig. 5). The first notable difference between these simulations and those with cells with a round target shape is a more gradual decrease in diffusion coefficient as Jcell increases (Fig. 5a), whereas previous monolayer simulations showed an abrupt drop in D associated with the onset of sub-diffusion (Fig.  2b,c). The onset of sub-diffusion here is not seen until Jcell = 2.75 (Fig. 5b). Like previous monolayer simulations (Fig. 2c), the drop from normal diffusion (⍺ ≈ 1) to sub-diffusion (⍺ ≈ 0.4) is abrupt. A steady decrease in observed cell shape index is seen as a function of increasing Jcell (Fig. 5c), approaching the jamming threshold of 3.81. As in cell shape quantification in earlier simulations (Fig. 3), the onset of sub-diffusion does not perfectly correspond to the expected jamming threshold as assessed by cell shape index -however, there is a slight drop and then plateau in cell shape index at Jcell ≥ 2.75. The appearance of the crossover point in neighbor exchanges shifts to nearly match the Jcell value associated with the onset of sub-diffusion, with significant neighbor exchanges resulting in crossover at all Jcell < 2.5 (Fig. 5d). Taken together, these metrics demonstrate that the increase in target cell shape index pushes the onset of jamming to high values of Jcell associated with significant penalties for disrupting initial cell-cell contacts.
The shift of this transition can be further rationalized by considering the effect of target perimeter variation on the interaction between cells. As in other work 46 , we can define an interfacial tension, , by taking the derivative of the CPM Hamiltonian with respect to perimeter: (For a detailed discussion of this equation, see Supporting Text 2). For round cells (target shape index 3.55), the second term of this equation is approximately 0, as target perimeter is effectively reached in all cases, and the interfacial tension is given simply by ≈ J. However, for cells with a larger target perimeter, the second term is typically negative, lowering the interfacial tension between cells to values that permit cells to remain motile. Examining the average interfacial tension per cell across both elongated and round cells reveals that there is a common value for interfacial tension associated with the onset of sub-diffusion (Fig. 5e).
To further characterize the dynamics of CPM monolayers, we more closely examined the cell migration seen in such monolayers. In round cells, it is apparent that the transition from mobile where v and v' denote velocity vectors of different cells, r is the distance between cells, and brackets denote averaging over all lag times of 1000 MCS and all pairs of cells. In monolayers consisting of round cells, an increase in spatial correlation between neighboring cells is observed as the jamming transition is approached (Fig. 6a). Spatial correlations vary in accordance with the average distance between neighboring cells and become longer range in the vicinity of the jamming transition, indicating collective cell motion. The correlation is maximized at Jcell = 1, closely matching the adhesion parameter value associated with the sharp drop in diffusion coefficient and abrupt transition to sub-diffusion (Fig. 2b,c). At higher Jcell values, spatial correlations in velocity decrease as the system moves deeper into the jammed state (Fig. 6a). For cells with a larger target cell shape index, no spatial correlations in velocity are seen in the vicinity of the jamming transition nor are any strong correlations found across the range of Jcell values explored (Fig. 6b). Indeed, the mixed monolayer remained fluid, as indicated by larger diffusion coefficients relative to round cell monolayers (Fig. 7a, Fig. 2b) an anomalous diffusion exponent near 1 (Fig. 7b,  Even at the smallest fraction of elongated cells explored here (20%), cells are only mildly subdiffusive, distinct from the deeply sub-diffusive behavior observed in cell monolayers in which the target shape is a circle.
To further probe cell dynamics in these mixed monolayers, spatial correlations in cell velocities were assessed. Correlations in velocities decrease across the monolayer with increasing percentage of elongated cells (Fig. 7c). In a mixed monolayer, strong correlations are only observed for nearest neighbor cells, with weak correlations observed for next nearest neighbors at lower percentages of elongated cells. At all elongated cell percentages, correlations in space more quickly decay to zero relative to the purely round cell case. It is apparent that the ability of cells to elongate inhibits jamming while decreasing the degree of correlated motion, abolishing the apparent transition to collective motion that precedes the jamming transition.  (Table S2), supporting the existence of a jamming transition in this context.
The effects of confinement are reflected in the magnitude and character of mean square displacements in each system (Fig. 2). Isolated cells remain normally diffusive at all adhesion energies (Fig. 2a), while aggregate cell sub-diffusion is dependent on the strength of environmental confinement (Fig. 2b). In contrast, monolayers experience an abrupt transition to deeply sub-diffusive behavior and multiple orders of magnitude decrease in diffusion coefficients across this same range of Jcell values (Fig. 2c). These differences highlight that abrupt cell arrest is unique to monolayers and emphasizes how degree of confinement affects cell dynamics.
Importantly, sub-diffusion and arrest are not solely linked to energy barriers dictated by the adhesion parameter J, as shown by the normally diffusive character of isolated cell migration.
Sub-diffusion has previously been used as the primary readout for the onset of jamming within the CPM 42 , and here we again find it to be a generally reliable indicator in the confluent monolayer. However, in the set of simulations described here, we employed additional metrics to characterize the onset of the jamming transition to provide a fuller framework for analysis of CPMbased jamming, as well as avoid mischaracterization of situations in which sub-diffusion is not associated with jamming, as observed in aggregates. We considered not only the onset of subdiffusion but also cell shapes, cellular rearrangements, and spatial velocity correlations. These metrics support the existence of jamming-induced sub-diffusion in monolayers and were generally in excellent agreement as predictors of the jamming transition (Table S2). Equilibrium cell shape in the CPM is moderately dependent on the value of Jcell, as cells tend to alter their perimeter to either minimize or maximize contacts with the surroundings. Therefore, as adhesion penalty is increased, a decrease in cell shape index occurs as cells reduce contact length, pushing the cells towards jamming. Shape index distributions follow trends seen in other systems, with the shape index of cells with a round target shape reaching and dipping below the previously identified jamming-associated shape index of 3.81 at the onset of sub-diffusive behavior (Fig. 3a) 34 .
Simulations of cells with a higher target shape index, elongated cells, do not show perfect correspondence between the onset of jamming and cell shape index dipping below 3.81 (Fig. 3b), though trends in average shape index and distribution are consistent with previous findings, with a narrowing distribution and stabilization of the mean cell shape index at a value of ~3.9 for cells that exhibit other features associated with jamming.
While cell shape analysis did not perfectly conform to expected thresholds of jamming identified in other systems, trends in cell shape were in line with previous results. As expected, increasing cell shape index via an increase in the target perimeter resulted in very effective fluidization of otherwise jammed cells (Fig. 5). Although jamming in the CPM is highly dependent on cell-cell adhesion as set by Jcell, as shown in this work and others 42 , increase in target perimeter (without any change in J values) provides cells with sufficient degrees of freedom to overcome otherwise prohibitively high energy barriers associated with cell-cell rearrangements in the monolayer. Through neighbor exchange analysis, we provided a clear measure of cell rearrangements within the CPM and confirmed that the absence of such rearrangements accompanies the abrupt transition to sub-diffusion that was observed (Fig. 5e). Notably, there is a common interfacial tension associated with jamming that applies across cells of different morphologies (Fig. 5f). Such a universal descriptor of jamming could be used to predict jamming across a variety of contexts not limited to the CPM and opens the door for further comparison to experimental data and other in silico models of epithelial systems.
Collective migration is inherently intertwined with the jamming transition 10,18,27 , in that it allows cells to continue to remain mobile when local caging sets in and energy costs to neighbor exchanges become high [61][62][63] . In this work, such motility manifests in migration of round cells in epithelial monolayers (Video 3, Fig. 6a). The results presented in this work highlight the ways in which cell shape influences the onset of jamming. When forcing cell elongation through an increase in cell target perimeter, a significant shift in the onset of jamming as a function of adhesion energy was observed. It is important to note that, while the onset of jamming is clearly pushed to higher values of Jcell, encoding a large target perimeter does not lead to a large increase in magnitudes of monolayer MSDs or diffusion coefficients. As can be seen in Figs. 2-4, the diffusion coefficients are very similar across elongated and round cell type monolayers of a given Jcell value. Rather, the main difference between the two cases is seen in the velocity correlations. Velocity correlations are effectively abolished in the monolayer with cells of higher target cell shape index, which gives cells the ability to deform (Fig. 6b). This same effect is seen in the mixed monolayers, where even relatively small percentages of elongated cells in the monolayer significantly reduces long-range velocity correlations (Fig. 7c). Unlike rounder cells, monolayers containing elongated cells display greatly reduced collective dynamics, maintaining more disordered motion until the onset of complete arrest at very high Jcell values, as they retain sufficient freedom of mobility to not have to rely on collective motion.
The results presented here have implications for the study of metastasis, particularly when considering cancer invasion as an unjamming transition. Cancerous tumors contain a diverse mixture of cell types, and the complex interplay between cells in this heterogeneous environment can influence the occurrence of jamming and unjamming transitions. Prior works have considered the interactions between potentially invasive, mesenchymal cells, and non-invasive epithelial cells, which will have tendencies towards unjamming and jamming, respectively 16 . Notable features of invasive cancer cells, which have progressed through the EMT, are elongated morphologies and softer, more deformable cell bodies. It is well-documented that transformed cells are consistently softer than healthy cells, which aids in their invasion by facilitating more efficient cell shape change that allows squeezing around obstacles and/or through tight spaces 47,64-66 . Our results on elongated cells within the CPM suggest one way in which a deformable cell body can facilitate motility and invasion in highly crowded environments such as tissue.