Interplay of curvature and rigidity in shape-based models of confluent tissue

Rigidity transitions in simple models of confluent cells have been a powerful organizing principle in understanding the dynamics and mechanics of dense biological tissue. In this work we explore the interplay between geometry and rigidity in two-dimensional vertex models confined to the surface of a sphere. By considering shapes of cells defined by perimeters whose magnitude depends on geodesic distances and areas determined by spherical polygons, the critical shape index in such models is affected by the size of the cell relative to the radius of the sphere on which it is embedded. This implies that cells can collectively rigidify by growing the size of the sphere, i.e. by tuning the curvature of their domain. Finite-temperature studies indicate that cell motility is affected well away from the zero-temperature transition point.

Rigidity transitions in simple models of confluent cells have been a powerful organizing principle in understanding the dynamics and mechanics of dense biological tissue. In this work we explore the interplay between geometry and rigidity in two-dimensional vertex models confined to the surface of a sphere. By considering shapes of cells defined by perimeters whose magnitude depends on geodesic distances and areas determined by spherical polygons, the critical shape index in such models is affected by the size of the cell relative to the radius of the sphere on which it is embedded. This implies that cells can collectively rigidify by growing the size of the sphere, i.e. by tuning the curvature of their domain. Finite-temperature studies indicate that cell motility is affected well away from the zero-temperature transition point.
Recent years have seen a growing interest in the way that mechanical interactions between cells play a fundamental role in structural and dynamical processes in biology [1,2]. This connection has been particularly apparent in the context of morphogenesis, where natural connections between mechanical stresses, cellular divisions, and the buckling and bending of epithelial sheets can be seen [2][3][4][5][6]. Simple coarse-grained models have been a natural way to approach such problems, ranging from latticebased models to soft spheres to deformable polygons to phase field models [7][8][9][10].
Here we focus on vertex models, which represent confluent monolayers as polygonal or polyhedral tilings of space; each geometrical unit corresponds to a coarsegrained cell [11] and the degrees of freedom are the vertices of the geometrical units. Vertex models attempt to explicitly represent mechanical interactions between neighboring cells by force laws that depend on the local geometry of the system, and have been used to model biophysical processes covering not only morphogenesis but also wound healing and tumor metastasis [12][13][14][15][16][17][18].
Such models have received attention not only for their appealing geometrical coarse-graining of clearly complex biological systems, but also for the unusual properties such models can support. For instance, two-dimensional vertex models have unusual zero-temperature rigidity transitions [19][20][21][22] with accompanying exotic mechanical states [23,24], their glassy dynamics at finite temperature can be deeply anomalous [25], and they can support unusual interfaces between coexisting populations of cells [26]. Although systematic and testable mappings from real confluent cellular systems to these very geometrical models are challenging, their unusual mechanical and dynamical properties suggest potential ways in which cells could exploit simple physical mechanisms to achieve unusual configurations or motions that may be useful for development.
While experiments on flat cellular monolayers are quite common, epithelial proliferation often takes place in domains where the curvature of the layer is strongly present and may be strongly varying (as in the ellipsoidal shapes of developing embryonic systems or in the regions of both positive and negative curvature in branching morphogenesis). While gradients in curvature surely play an important role, we begin in this work by studying vertex models in domains of constant positive curvature, i.e., on the surface of the sphere. We are particularly interested in the interplay between the curvature of the cellular monolayer and the mechanical or dynamical state of the system. We will see that the curvature of the domain has natural consequences for the zero-temperature rigidity transition in such models, and that the finite temperature "glassy" behavior of cells feels this change in the underlying transition.
In moving to a three-dimensional embedding there are many natural extensions of the vertex model that could be considered [3,[27][28][29][30], for instance expressing the columnar nature of confluent epithelial cells by having separate polygons representing both the apical and basal surface of each cell with corresponding faces on the lateral sides. For simplicity we consider the so-called "3D apical vertex models" [27], which represent each cell only by an apical polygon whose vertices are not restricted to lie in the plane. To demonstrate the tight connection between curvature and rigidity, we further specialize our study to two-dimensional cells constrained to the surface of a sphere, as schematically illustrated in Fig. 1.
We begin by writing down the 2D vertex model energy functional, This energy depends on the area A i and perimeter P i of each of the N cells, indexed by i. The model parameters are the "preferred" geometric values, A 0 and P 0 , along with the area and perimeter stiffnesses k A and k P (here we assume the monodisperse case in which all cells have identical preferences). Biologically, A 0 is commonly taken to represent a combination of cellular incompressibility and the resistance of the monolayer to height fluctuations, and P 0 to represent a competition between tensions and adhesions acting between cells; more broadly this can be viewed as a minimal Taylor series expansion in geometrical properties that describes cellular matter rather than foams [31].
The density dependence of this model can be made transparent by exploiting the fact that in these models the cells completely fill space, i A i = A total , and by choosing the unit of length to be the average area of the cells, A [9,32]. Using a and p to refer to the dimensionless area and perimeter, and letting k r = k A A /k P , the above energy can be written as Thus, if P 0 is a control parameter (and if the stiffnesses k r and P 0 are themselves density-independent), the density dependence of the model in flat space enters via p 0 = P 0 / A = P 0 ρ 1/2 . The parameter A 0 couples to the total size of the system, serving as an offset to the pressure of the system but not affecting the forces between degrees of freedom [32]. Extending the above expressions to a spherical vertex model requires no change of notation (although justifying the geometric coarse-graining likely requires a more biologically informed derivation, as we discuss in the conclusion). We simply interpret the "areas" and "perimeters" to be those measured on the sphere: perimeters are given by sums of geodesic distances as one traverses the vertices composing the cell, and areas are given by the area of the spherical polygons enclosed by those geodesic arcs. The forces acting on the vertices are given by the negative spherical gradients of 4 (explicit expressions are given in the appendix). The statistical mechanics of fluids confined to manifolds of constant curvature is itself a rich topic [33,34], and we benefit from methods built to understand such fluids in cases where the forces depend not on the Euclidean distance between interacting units but rather the geodesic distances.
To implement efficient and highly scalable numerical simulations of the above equations, allowing T1 transitions to facilitate neighbor exchanges between cells and evolving the degrees of freedom under equations of motion ranging from energy minimization schemes to overdamped Brownian dynamics to self-propelled "active" dynamics, we combine the the GPU-accelerated frameworks described in Refs. [35,36]. We constrain the vertices of the apical polygons to lie exactly on the sphere, which we accomplish within the context of a standard projection operator formalism [37,38] as described in the Appendix.
We first directly probe the athermal rigidity transition of the spherical vertex model as a function of N and p 0 ; for simplicity (see the appendix) here we first focus on the k r = 0 limit of Eq. 2. We prepare between 100 and 500 initial configurations for each value of p 0 , seeded by randomly placing cell centers on the surface of the sphere and deriving the initial positions of the vertices from the convex hull of that point pattern. We perform a FIRE energy minimization of these configurations [39] to find the inherent state associate with each initial configuration. Note that, like it's counterpart in flat space, the spherical vertex model described here is extensively underconstrained; as such, we anticipate that the ground states of the model are mechanically stable only in the presence of residual stresses [21,23,32].
Thus, we estimate the rigidity transition for a given value of N by computing the fraction of minimized states, F (p 0 , N ), which minimize to an inherent state of zero energy (within numerical precision). The probability distribution of transition points is given by the derivative of this function; to take this derivative while suppressing noise, we convolve a linear interpolation of the F (p 0 , N ) with the derivative of a Gaussian with standard deviation related to the shape of F (p 0 , N ) (see Ref. [32] for details). We have done this for both for the standard (planar) and spherical vertex models at k r = 0, and the results are shown in Fig. 2.
Our results for the mean value of the transition for the planar vertex model, and the variance of the distribution, are consistent with previous studies [40], although we note that other simulations, based on the Surface Evolver package [41] and minimizing under a different protocol, have reported slightly different results [20]. As might be expected, the primary effect of approaching the thermodynamic limit in the planar case is to develop a more sharply peaked distribution about the N → ∞ limiting value, with very little change in the mean value of the distribution. In contrast, the effects of changing the size of the sphere relative to the typical size of each cell is readily seen in the way the distribution of the transition point not only sharpens but also shifts with N .
We find that the critical value of p 0 separating the mechanically rigid and floppy phases as a function of N closely tracks (but is not precisely equal to) the way in the which the perimeter of a unit-area regular pentagon varies on a sphere of total surface area N . Note that this value of p 0 forms a natural bound for the non-linear rigidity transition: sufficiently large cellular displacements require cells to exchange neighbors, on average cells have six sides, so during the T1 transition a spherical pentagon must be formed. If p 0 (N ) < p penta (N ) this configuration will cost energy, but the precise connection between this bound on the nonlinear behavior of rearrangements and the infinitesimal rigidity calculation shown in Fig. 2 remains unclear (as it does also in the planar case).
To show that this qualitative shift is neither just a result of the k r = 0 limit explored above (i.e., in the absence of a cost associated with deviations in target cell area) nor an artifact of the exotic mechanical states at zero temperature found in vertex models [21,23,24,40], we perform preliminary studies of the finite-temperature dynamics of disordered configurations of the spherical vertex model, with k r = 1. Shown in Fig. 3, we perform overdamped Brownian dynamics at fixed model parameters {p 0 , T, k r , A } while varying N , where the chosen p 0 is above value associated with the peak of the distribution in Fig. 2 for small N but below it for large N . One clearly sees an echo of the underlying change in the rigidity transition point, p * 0 (N ), in the way the dynamics shift from simple, fluid-like, diffusive dynamics at small N to transiently caged, glassy dynamics as the size of the sphere increases and the system more closely resembles a vertex model in flat space whose parameters place it in the glassy regime.
The fact that the mean of the rigidity transition shifts as the relationship between curvature and cell-size varies suggests a novel mechanism by which cells could collectively tune between mechanical phases as a function of their curvature. Models of 3D collections of cells in embryonic zebrafish development have shown the potential for coexistence between fluid-like behavior in regions of high curvature and solid-like behavior in regions of lower curvature [42]. Perhaps more relevant to this explicitly two-dimensional model, developing insect embryos look much like ellipsoidal versions of the right panel in Fig.  1, with regions of high and low curvature. Thus, although we currently neglect gradients in curvature, the curvature-dependent rigidity discussed here might be directly relevant in the modeling of such systems [43].
We also note that the connection between curvature and the ability to support mechanical stresses suggests a relationship between density and jamming that is qualitatively different from the planar case. This comes from the competition between the scaling of the critical perimeter with typical cell size, p 0 = P 0 / A , and the effect of curvature as expressed by ratio of the sphere radius to the typical cell size. The situation is schematically depicted in Fig. 4, which shows an estimate of the shifting of the rigidity transition P c as a function of A . In one limiting case, the number of cells could increase on a sphere of fixed radius. In this scenario, the decrease of A lowers the critical transition point, so cells dividing (at constant other model parameters) could induce the system to collectively unjam via growth. In the other limiting case, the number of cells could increase in a simultaneously enlarging spherical domain (so that the cell number increases at fixed A ). In such a case the system could potentially rigidify via growth.
Recent work has suggested that real monolayers of epithelial cells in curved space may adopt configurations in which the apical and basal surfaces of a cell have very different geometries [5,44,45]. Whether the apical vertex models considered in this and related works are still sufficiently expressive coarse grained models to capture the underlying physics requires further work; it may be that in curved space models written only in terms of a single cross-sectional plane are sufficient only when the individual cells are small enough to not appreciable feel the effects of curvature. Before considering such complications, interesting extensions of the spherical vertex model presented here are anticipated by some existing studies of apical vertex models on curved surfaces [3,29,30], in which the curved space is not a fixed embedding but can itself evolve and deform as the cells collectively exert stresses on their environment. For systems with the topology of a sphere, we expect that gradients in curvature could themselves be a source of the residual stresses that are necessary to rigidify the sorts of extensively underconstrained systems described by Eq. 4.
Additionally, it will be very interesting to investigate the finite temperature glassy dynamics of this model in greater detail. Previous work on planar vertex and Voronoi models identified a deeply anomalous type of "sub-Arrhenius" dynamics, in which the relaxation time of the cells grew more slowly than exponential with decreasing temperature [25]. One speculation relates these unusual glassy dynamics to the unusual, residual-stressdriven rigidity transition those models possess at zero temperature. Embedding the vertex model on a sphere, as we have done here, provides one way to formally probe this hypothesis. We have constrained the vertices to lie exactly on the sphere, but have included no other energetic terms related to the curvature. Previous works on the apical vertex model have taken Eq. 4 and supplemented it with an additional term of the form [30] where i and j run over all neighboring faces andn i is the surface normal corresponding to cell i.
In going to 3D one would have to carefully treat the finite thickness of the monolayer -together with the way curvature may interact with the biological processes that give rise to adhesions, tensions, etc. -to derive the appropriate bending energy for these models, but the expression above is still formally useful to us here. On a unit sphere note that the geodesic distance between two points is |r ij | = cos −1 (n i ·n j ); in the limit where the inter-cellular spacing is small compared to the radius of the sphere the above term can be approximated by E b ≈ B 2 ij r 2 ij , adding an additional quadratic constraint for ever pair of cellular neighbors. Whereas Eq. 4 represents an extensively underconstrained system that can only rigidify through residual stresses, Eq. 2 introduces enough additional constraints to rigidify the system more conventionally. Thus, studying the glassy dynamics of the spherical vertex model as a function of tuning B from zero to unity could test the root cause of the anomalous glassy behavior seen in other simple models of cellular matter.
I would like to thank Matthias Merkel and Michael Moshe for illuminating discussions and for critical readings of this manuscript.

Projection operator formalism
We follow Refs. [37,38] in using the projection operator formalism to enforce the hard constraint that the degrees of freedom lie on the surface of a sphere. To take a specific example of an equation of motion used in this manuscript, to numerically implement overdamped Brownian dynamics at temperature T we write where µ is an inverse friction coefficient, η a normally distributed random force with zero mean and with η iα (t)η jβ (t ) = 2µT ∆tδ ij δ αβ in each of the three Cartesian directions denoted by greek indices (so that the noise has the correct statistics in the tangent plane of the vertex). The operator P T (a, b) = b − (â · bâ) projects the forces and the random noise onto the tangent plane at the location of the degree of freedom. To maintain the spherical constraint small time steps must be used, and degrees of freedom are projected back onto the surface of the sphere of radius R after they are moved via r i (t + ∆t) = P N (r i (t) + ∆r i ) for P N (a) = R a |a| .

Force calculations
For completeness we explicitly write some of the expressions used to compute forces in the vertex model whose energy is given by and where the degrees of freedom are constrained to lie on the surface of a sphere of radius R. We note that on the surface of the sphere there are many equivalent expressions for the geodesic distance between two points (or the included angle between three points, or the area of a spherical polygon given by n points, etc). While analytically equivalent, these expressions typically have different regimes of numerical stability. For instance, given two points on the sphere, n 1 and n 2 , the distance d may be written as The first expression above is perhaps the simplest and least computationally expensive, but it is poorly conditioned for very small distances (as might be relevant when vertices get very close to each other before performing a T1 transition). The second expression is poorly conditioned for large distances, whereas the third is the most computationally expensive but is well-conditioned for all distances. These questions of numerical stability become especially acute when dealing with the forces, and we have found it important to implement self-consistency checks on the force calculations and switch to analytically equivalent but numerically different routes of calculating gradients in the spherical vertex model.
From this one readily appreciates the substantial cost of computing gradients of d b or d c , and whenever possible we opt for the simpler expressions stemming from d a .
Similarly, note that the area of a spherical triangle A( n 1 , n 2 , n 3 ) can be written as a = d( n 2 , n 3 )/R, c = d( n 1 , n 2 )/R.
Clearly, again, care must be taken in choosing distance functions that will lead to well-conditioned expressions for both the area and gradients of the area while also minimizing the complexity of the resulting expressions. Additional considerations include the efficiency and numerical stability of computing cellular areas as either the sum of spherical triangles formed by the cell centroid and consecutive vertices around the cell or via the sum of the included angle at each of the n vertices, A({ n 1 , n 1 , , . . . , n n , }) = R 2 n i α i − (n − 2)π ; (17) this is particularly delicate when the cells are not convex, or when edges cross.

Initial conditions
The initial conditions in this work are chosen to be high-temperature random configurations of cells which we then either quench to zero temperature or perform finite-temperature overdamped brownian dynamics on. We accomplish this by first picking a desired number of cells, N c and distributing N c points uniformly on the surface of the sphere. We use the Computational Geometry Algorithms Library (CGAL) [46,47] to construct the convex hull of these points, and take the initial vertex positions to be the centroids of the resulting facets (projected back onto the sphere) [48].