Relating cell shape and mechanical stress in a spatially disordered epithelium using a vertex-based model

Abstract Using a popular vertex-based model to describe a spatially disordered planar epithelial monolayer, we examine the relationship between cell shape and mechanical stress at the cell and tissue level. Deriving expressions for stress tensors starting from an energetic formulation of the model, we show that the principal axes of stress for an individual cell align with the principal axes of shape, and we determine the bulk effective tissue pressure when the monolayer is isotropic at the tissue level. Using simulations for a monolayer that is not under peripheral stress, we fit parameters of the model to experimental data for Xenopus embryonic tissue. The model predicts that mechanical interactions can generate mesoscopic patterns within the monolayer that exhibit long-range correlations in cell shape. The model also suggests that the orientation of mechanical and geometric cues for processes such as cell division are likely to be strongly correlated in real epithelia. Some limitations of the model in capturing geometric features of Xenopus epithelial cells are highlighted.


Introduction
Many essential aspects of cell behaviour are controlled, both directly and indirectly, by mechanical cues (Huang & Ingber, 1999;Wozniak & Chen, 2009). For example, cell density and substrate adhesion have been shown to affect cell proliferation (Huang & Ingber, 2000;Streichan et al., 2014), while cell division orientation appears to be regulated by mechanical feedback (Théry & Bornens, 2006;Fink et al., 2011;Minc et al., 2011;Wyatt et al., 2015). Many morphogenetic processes, such as gastrulation and convergent extension (Martin et al., 2009), are mechanical processes inducing significant changes to the stresses within the tissue (Lecuit & Lenne, 2007). However, despite its significance in development, the mechanical state of tissues remains poorly characterized in comparison to some aspects of genetics and biochemical signalling.
The geometric properties of cells are governed by cell adhesions and cytoskeletal mechanics (Kiehart et al., 2000;Käfer et al., 2007), which in turn feed into global tissue dynamics (Shraiman, 2005;Martin et al., 2009;Guillot & Lecuit, 2013). The mechanical state of an individual cell is largely dependent on its interaction with its neighbours and adhesion to the extracellular matrix. Experimental techniques such as laser ablation (Hutson et al., 2003;Campinho et al., 2013;Mao et al., 2013) and atomic force microscopy (AFM) (Hoh & Schoenenberger, 1994) have been used to characterize cell mechanics; laser ablation reveals cell-level forces by making small slices in the tissue and observing the recoil velocity of cells, while AFM attempts to deduce the local mechanical properties of a tissue by performing small indentations using a mechanical cantilever. While revealing, such experimental techniques are invasive and typically require modelling for the interpretation of measurements. Live fluorescent imaging combined with high resolution microscopy offers alternative insights into developmental processes such as gastrulation (Rauzi et al., 2008;Heller et al., 2016). Measurements of cell shape over time allows inference of mechanical stress (Chiou et al., 2012;Ishihara & Sugimura, 2012;Xu et al., 2015Xu et al., , 2016, based on an underlying mathematical model. This non-invasive approach has led to significant growth in mathematical modelling of epithelial cell mechanics in 2D and 3D (Hilgenfeldt et al., 2008;Brodland et al., 2010;Okuda et al., 2013;Hannezo et al., 2014;Collinet et al., 2015;Bielmeier et al., 2016;Sugimura et al., 2016;Tetley et al., 2016). However without direct measurements of stress, mechanical predictions taken from geometric data alone are only as good as the constitutive models from which the predictions are derived.
Theoretical models of epithelial mechanics fall into a number of classes, including cellular Potts (Graner & Glazier, 1992), cell-centre (Osborne et al., 2010), vertex-based (Nagai & Honda, 2001;Farhadifar et al., 2007;Staple et al., 2010;Fletcher et al., 2014) and continuum models (Edwards & Chapman, 2007;Nelson et al., 2011). Vertex-based models exploit the polygonal shape commonly adopted by tightpacked cells in a monolayer, characterizing the monolayer as a network of cell edges meeting (typically) at trijunctions. Typically, vertices are assumed to move down gradients of a mechanical energy, often subject to a viscous drag; the network topology changes intermittently as cells intercalate, divide or are extruded. It is of interest to relate such cell-level models, describing cells as individual entities that can evolve at discrete time intervals, to continuum models describing the smooth changes of a tissue in space and time. Some progress has been made in upscaling spatially periodic cell distributions in 1D (Fozard et al., 2010) and 2D (Murisic et al., 2015) using homogenization approaches, or by direct coarse-graining (Ishihara et al., 2017). Simulations have revealed striking properties of more realistic disordered networks in 2D (Staple et al., 2010;Bi et al., 2015), such as a rigidity transition characteristic of a glassy material. Abundant imaging data makes parameter estimation feasible, allowing models to be tested quantitatively and used to explore new biological hypotheses.
In this article, working in the framework of a popular vertex-based model describing a planar monolayer of mechanically (but not geometrically) identical cells, we derive expressions for the stress tensor at the cell and tissue level, and use these results to understand the relationship between a cell's shape and its mechanical environment, showing that the principal axes of the cell's stress and shape tensors align. We parameter-fit simulations to images of Xenopus embryonic epithelia, using cell area over polygonal classes as a measure. Of particular interest is the manner in which mechanical effects constrain the spatial disorder that is intrinsic to epithelial monolayers, which we characterize using simulations, highlighting the appearance of spatial patterns reminiscent of force chains in granular materials. We also discuss the role of the stress acting on the monolayer's periphery in determining the size and shape of cells.

Experiments
Experimental data were collected using tissue from the albino Xenopus laevis frog embryo. Animal cap tissue was dissected from the embryo at stage 10 of development (early gastrula stage) and cultured RELATING SHAPE AND STRESS IN EPITHELIUM i3 Fig. 1. Experimental setup and data analysis. (a) Animal cap tissue was dissected from stage-10 Xenopus laevis embryos and cultured on PDMS membrane. (b) Side-view confocal image of the animal cap (top:apical; bottom:basal), stained for microtubules (red), beta-catenin (green) and DNA (blue). A mitotic spindle is visible in the centremost apical cell. The animal cap is a multilayered epithelial tissue; we analyse just the outer, apical, cell layer. (c) The apical cell layer of the animal cap tissue is imaged live using confocal microscopy (green, GFP-α-tubulin; red, cherry-histone2B). (d) The cell edges are manually traced and cell shapes are derived computationally, being polygonized using the positions of cell junctions. (e) Mean normalized area as a function of polygonal class showing mean and one standard deviation, from experiments (solid and shaded) and simulation (dashed) with parameters , as shown with P ext = 0. Cell areas were normalized relative to the mean of each experiment. (f) Circularity as a function of polygonal class showing mean and one standard deviation, from experiments (solid and shaded) and simulation (dashed) using the same parameters as in (e). (g) Proportions of total cells in each polygonal class in experiments (left bar) and simulations (right bar). Error bars represent 95% confidence intervals calculated from bootstrapping the data. (Colour in online.) on a 20 mm × 20 mm × 1 mm, fibronectin-coated, elastomeric PDMS substrate (Fig. 1a). The animal cap tissue is a multi-layered (2-3 cells thick) epithelium (Fig. 1b), which maintains its in vivo structure when cultured externally for the time period of our experiments (up to five hours). This system has the advantage of closely resembling in vivo tissue whilst also giving the ability to control peripheral stress on the tissue. For this work, a 0.5 mm uniaxial stretch was applied to the PDMS substrate, which ensured that it did not buckle under gravity or the weight of the animal cap. This small stretch was found to have no measurable effect on cell geometry (data not shown) and we therefore assume that there is negligible i4 A. NESTOR-BERGMANN ET AL.
peripheral stress on the tissue. The apical cell layer of the animal cap tissue was imaged using a Leica TCS SP5 AOBS upright confocal microscope (Fig. 1c) and cell boundaries were segmented manually (Fig. 1d), representing each cell as a polygon with vertices coincident with those in images. The vast majority of vertices were classifiable as trijunctions.
Letting a cell, α, have Z α vertices defining its boundary, we characterize the shape of the cell using its areaÃ α and shape tensor,S α , defined with respect to cell vertices as whereR i α is the vector running from the cell centroid to vertex i andẑ is a unit vector pointing out of the plane.S α has eigenvalues (λ α1 , λ α2 ) with λ α1 ≥ λ α 2 > 0. The eigenvector associated with the larger (smaller) eigenvector defines the major (minor) principal axis of cell shape, the two axes being orthogonal. The circularity parameter C α = λ α2 /λ α1 ∈ (0, 1] indicates how round a cell is. The variation of cell area and circularity across an individual monolayer is illustrated in Fig. 1(e and f), distributed across the cells' polygonal class Z α (number of neighbours). The distribution of cell number across polygonal class is shown in Fig. 1 (Fig. 1f). As explained below, we used A exp to fit parameters of the vertex-based model (Fig. 1e).

The vertex-based model
In this section, we derive expressions for cell and tissue stress using the vertex-based model and describe our simulation methodology. We explain relationships between cell stress and cell shape and discuss the mechanical properties of the monolayer.

Geometry of the monolayer network
We represent an epithelial monolayer as a planar network of N v vertices, labelled j = 1, . . . , N v , connected by straight edges and bounding N c polygonal cells, labelled α = 1, . . . , N c . The vector from the coordinate origin to vertex j is given byR j (t); here tildes denote dimensional variables andt is time. Quantities specific to cell α are defined relative to its centroidR α . Cell α has Z α vertices labelled anticlockwise by i = 0, 1, 2, . . . , Z α − 1 relative toR α . We defineR i α as the vector from the cell centroid to vertex i, such unit vectors along a cell edge byt i α and outward normals to edges byñ i α =t i α ×ẑ. The lengthl i α of an edge belonging to cell α between vertices i and i + 1, and the cell perimeterL α , are given bỹ The cell area (assuming convex polygons),Ã α , and shape tensor,S α , are given by ( coordinate origin, and have unique latin superscripts, j i.e.R j . The matrix capturing the mapping from the vertex labels, j, to the vertex labels, i, within every cell, α, is defined as such that, for an internal vertex, Nv j=1 c ij αR j =R α +R i α . For trijunctions, there will exist α, α , α for respective i, i , i such that c ij α = c i j α = c i j α = 1, for a given j. A visual representation of this geometric arrangement is given in Fig. 2.

Cellular forces and energies
We adopt a well-established and widely used vertex-based constitutive model (Honda & Eguchi, 1980;Nagai & Honda, 2001;Farhadifar et al., 2007;Mao et al., 2013;Fletcher et al., 2014;Bi et al., 2015). We consider a monolayer of cells with identical physical properties but differing in general in size and shape. Every cell is assumed to have a mechanical energy,Ũ α , defined bỹ 3) The first term in (3.3) models the cell's bulk compressibility, in terms of a preferred areaÃ 0 and a stiffness K. The remaining terms represent the contractility of the cell periphery, via cortical actomyosin bundles and cell-to-cell adhesion. The parameter˜ represents the contractile strength while˜ tunes the effective preferred cell perimeterL 0 = −˜ /2˜ , such that the energy associated with the peripheral forces is of the form 1 2˜ (L −L 0 ) 2 . The quadratic contributions to the energy as a function of perimeter and area could in principle be extended with higher-order non-linearities. At the tissue level, the system is assumed to evolve down gradients of the bulk energy Nc α=1Ũ α from an initial disordered state. We model the deterministic i6 A. NESTOR-BERGMANN ET AL. evolution by assigning a drag force (relative to the substrate on which the monolayer sits), to each vertex of cell α, of the form −η(Ã/Z α )dR j /dt. The drag magnitude is chosen to scale with the cell's area rather than its number of vertices (a natural assumption if the drag arises from physical interactions distributed across the base of the cell) and viscous resistance to internal shear or extension is neglected. For the time being we do not consider topological rearrangements of the network of cell edges, but return to this when discussing simulations in Section 4.
We non-dimensionalize by scaling lengths on Ã 0 , using where L 0 = − /2 is the dimensionless preferred perimeter. The total energy, U, of the monolayer may now be written as the sum where U 0 = 2 /4 2 is a constant that may be discarded as the dynamics are driven by energy gradients. For later reference we define an associated pressure and tension for each cell as Cellular forces can be computed directly from the mechanical energy, using the fact that δ i U α = ∇ i U α · δR i α . The first variation of the energy with respect to the position of vertex i is given by can be interpreted as the force required to shift vertex i through δR i α to do work δ i U α ; equivalently, f i α represents the restoring force exerted at vertex i by cell α. This force can be calculated explicitly by differentiating the mechanical energy term by term. Considering first the area contribution we find where P α is given by (3.7a) and gives the direction of the bulk compressive force at node i. The perimeter term gives where T α (see (3.7b)) represents a tension and q i α ≡t i α −t i−1 α represents the direction of the inward force due to stretching of the cell perimeter. Thus the force at vertex i can be written The analogous force for a vertex model lacking theL 2 α term in (3.3) is given in Spencer et al. (2017). f i α represents the force generated when perturbing the vertex of a cell in isolation. For the case of a monolayer, each vertex will have a contribution from the three cells attached to it (or fewer, if the cell is at the periphery of the monolayer). Thus the net force on vertex j, f j , will be given by the sum of the contributions from each cell attached to it as where c ij α ensures that, although the summation is over all cells, we count only the contributions from the cells connected to vertex j. More specifically, if cells α, α and α meet at junction j, with anticlockwise tangents t, t , t emerging from the vertex with normals (pointing clockwise) n, n , n orthogonal to each tangent, the net force at the vertex can be written The tangential forces show how each edge is a composite structure with tension contributions from two adjacent cells. The factor of 1 2 in the pressure terms reflects the fact that the force due to pressure acting on any edge is distributed equally between each vertex bounding the edge. The tensions and pressures depend on the total area and perimeter of each neighbouring cell via (3.7). For vertices at the periphery of the monolayer, bordering cells α and α , we write P α = P ext (an imposed isotropic stress) and set T α = 0, so that (3.14) We use this relationship below when considering the boundary conditions at the edge of a monolayer.
When the system is out of equilibrium, the net force at each internal vertex is where the term proportional toṘ j is a viscous drag having contributions from the three cells at the trijunction; the dot denotes a time derivative. WritingṘ j =Ṙ α +Ṙ i α , the drag can be considered as representing an internal dashpot within each cell connecting the cell centre to the vertex plus a drag on each cell centre. Thus the net force on cell α becomes Since inertia is negligible, the net force on any vertex and on any cell must vanish, F j = 0 and F α = 0. The former condition defines the N v coupled evolution equations of the network vertices. When the system is in equilibrium, this simplifies to f j = 0, f α = 0. Likewise the net torque on cell α,

The stress tensor of a cell
For a tensor σ that is symmetric and divergence-free, defined over an area A with perimeter S , we have where R is an arbitrary position vector. Thus taking an area integral and applying the divergence theorem gives (Norris, 2014) We use this weak formulation to derive the stress tensor of the monolayer, taking the stress to be uniform over each cell. The forces acting on cell α are distributed around the vertices, so that taking A = A α (the domain of cell α), (3.18) motivates the definition of the cell stress σ α as This reveals conservative (elastic) and dissipative (viscous) contributions to the stress. The former is If the cell is in equilibrium and under zero net torque, then Zα −1 i=0 R i α × f i α = 0 (see (3.17)), ensuring that this contribution to σ α is symmetric; the symmetry of (3.20) is confirmed below. Likewise the absence RELATING SHAPE AND STRESS IN EPITHELIUM i9 of torque on a cell due to drag in (3.17) requires the dissipative component of the stress to be symmetric, allowing us to redefine the final term in (3.19b) as where S α is the dimensionless shape tensor based on vertex location. We simplify (3.20) by making use of two geometric identities, established in Appendix A, namely we can then express the stress of cell α as Here the elastic components of the stress have been written in terms of an isotropic and deviatoric component. The former defines the effective cell pressure, which has contributions from the cell's bulk and the perimeter (in Young-Laplace form, with an effective radius of curvature 2A α /L α ) as (3.24) We will see below how the competition between bulk pressure and cortical forces can stiffen the monolayer. The traceless contribution to the cell stress is (3.25)

Relating cell stress and shape
We can now explore the relationship between the principal axes of cell shape and stress by considering the commutativity of σ α and S α . The tensors will share an eigenbasis, implying that their principal axes align, if and only if they commute. Having separated the stress tensor (3.23) into an isotropic and deviatoric component however, we require only that S α J α = J α S α and S αṠα =Ṡ α S α , which is established via direct algebraic manipulation in Appendix B. Figure 3 provides a computational illustration of this mathematical result for a disordered monolayer in equilibrium; details of the simulation scheme are given in Section 4. Thus, for an individual cell, the principal axes of stress and shape align (both quantities being defined directly in terms of cell vertex locations). Equivalently, within the present model, cells that are elongated experience a local stress field that is oriented exactly with the direction of elongation. The consequences of this observation are discussed below.  The initial cell array was generated using a Voronoi tessellation and then relaxed to equilibrium using periodic boundary conditions. The eigenvectors corresponding the the principal eigenvalue of σ α and S are plotted in black and yellow, respectively. Darker cells have P eff α > 0 (net tension); lighter cells have P eff α < 0 (net compression). (Colour in online.)

Stress of the monolayer
We now return to (3.18), taking the domain A in (3.18) to cover multiple cells. The area integral can be evaluated over each cell to give a formulation for the 'tissue' stress over a simply connected region R of the monolayer σ R as summing over cells in R. The components of the first two terms on the right-hand side of (3.23) that are proportional to T α at the cell level, and their area-weighted sum in (3.26), are analogous to an expression derived by Batchelor (1970) for a suspension of particles having interfacial tension. Equivalent expressions for the equilibrium stress of the present model based on Batchelor's formulation have been given by Ishihara & Sugimura (2012) and Guirao et al. (2015). For now let us take R to be the whole monolayer. The line integral in (3.18) can be evaluated by setting since F j = 0 at all internal vertices. Let k = 0, 1, . . . , N p − 1 label the peripheral vertices, let peripheral normals n k−1 and n k border vertex k and let p k = 1 2 (n k−1 + n k ). Since the periphery is a closed curve, its sum of tangents vanish, hence its sum of normals vanish, hence Np−1 k=0 p k = 0. Let R 0 be the centroid of the monolayer, and write R k = R 0 + R k 0 , so that Assuming the pressure is P ext uniformly around the periphery, the force balance at the peripheral vertices (3.14) gives where the final expression results from (3.22a) and A = Nc Taking the trace of this sum gives which describes the relaxation of the area of the monolayer to its equilibrium. Once in equilibrium, the system must satisfy A disordered distribution of cells within an equilibrium monolayer will have a range of values of P eff α , and non-isotropic cells will have deviatoric contributions to their stress, but the whole population must satisfy the weighted sums (3.31). For an isolated monolayer that is in equilibrium under zero external loading (the condition relevant to Section 2), we must therefore impose (3.32)

Elastic moduli
It has been shown that, when the cells are identical hexagons, the stress at the tissue level under the present model (neglecting friction) is equivalent to that of linear elasticity when considering small perturbations about the unstressed state (Murisic et al., 2015). We can therefore use the expressions for stress at cell (3.23) and tissue (3.26) level to recover expressions for the associated elastic moduli.
Taking P ext = 0 in a base state, imposing (3.32), we consider an isotropic expansion of a disordered monolayer of magnitude 1 + where 1, so that L α maps to (1 + )L α , A α maps to (1 + 2 )A α and so on. Linearizing about the base state, the dimensional bulk modulus,KÃ 0 K, of the monolayer is given by using (3.24) and (3.31). This prediction holds for a disordered network of cells, and therefore provides a direct means of determining the variability of bulk modulus over different realizations of the monolayer. When simplified to the special case of a hexagonal monolayer, for which L α / √ A α = μ 6 ≡ 2 2 √ 3 ≈ 3.72 for all α, (3.33) reduces in dimensionless form to in agreement with Murisic et al. (2015) and Staple et al. (2010). K remains positive for < 0, but can become zero at = (8/μ 6 ) 2/3 when A = 1. The dimensional shear modulus,KÃ 0 G, for the special case of a monolayer of identical hexagonal cells, is shown in Appendix C to be given by which is also equivalent to the shear modulus derived by Murisic et al. (2015) (but differs, as they showed, with Staple et al. (2010)). Equation (3.35) illustrates how L must exceed L 0 , i.e. cell walls must be under tension, in order for the monolayer to resist shear. Prediction of the shear modulus for the disordered monolayer is much less straightforward; estimates (for a disordered dry foam) are reviewed in Kruyt (2007).
= −2μ 6 , < 0. (3.38) Analysis of the cubic √ AP eff 6 as a function of √ A reveals that it is monotonic (implying a single root of P eff 6 = 0) for > 2/μ 2 6 ; a positive root exists provided < 0 for > 0 that satisfies A * 6 = 1 along (3.38). For > 0, the cubic has repeated roots along As a consequence the parameter map shown in Fig. 4 can be drawn (Farhadifar et al., 2007), with the boundary between regions I and II a defined by (3.38), that between regions II a and III by = 0 and > 2/μ 2 6 and that between regions II b and III by (3.39). We will focus attention below on region II, in which at least one stress-free equilibrium state exists (for hexagons) with positive shear modulus. Along the region I/II a boundary, hexagons have P = 0 (A = 1) and T = 0 (L 0 = L = μ 6 ) and the monolayer loses any resistance to shear (from (3.35)). (In a disordered monolayer, the rigidity transition to a floppy region-I state has been shown to arise closer to L 0 = μ 5 ≈ 3.81 (Bi et al., 2015).) Approaching the region II a /III boundary, the equilibrium cell area approaches A = 0; two possible equilibria exist in region IIb, coalescing at positive A along the region II b /III boundary.
For later reference, we note that for a periodic array of hexagons under an external load P ext (for which P eff α = P ext in (3.31)), we may define (for P ext > −1) such that if P eff 6 (A; , ) = P ext in (3.37) then P eff 6 (A † ; † , † ) = 0. This simple scaling symmetry of (3.37) allows the axes of Fig. 4(a) to be replaced with † and † in order to encompass externally-loaded monolayers subject to non-zero P ext . Figure 4(b and c) illustrates four distinct classes of equilibrium cell shape and stress that arise in simulations of disordered monolayers, distinguished by the signs of the eigenvalues (σ α,1 , σ α,2 ) of the cell stress tensor; recall that the corresponding eigenvectors align with the principal axes of the shape tensor S α . When P eff α = −Tr(σ α ) ≡ −(σ α,1 + σ α,2 ) > 0 (represented by darker cells, Fig. 4(b), the cell is enlarged and under net tension: both eigenvalues of the stress tensor are negative when the cell is rounder, although one can be positive when the cell is more elongated. Likewise when P eff α < 0 (lighter cells, Fig. 4(c), the cell is smaller and under net compression: both eigenvalues of the stress tensor are positive when the cell is rounder, although one can be negative when the cell is more elongated.

Simulation methodology
The majority of computational modelling was performed in Python, with some processes sent through C where Python struggled with performance. The cells were described as an oriented graph using the graph-tool module for Python (Peixoto, 2015). The algorithms and core data structures of graph-tool are written in C++, thus its performance in memory and computation is comparable to that of pure C++. The energy minimization was performed using a conjugate gradient method from the scipy library.
Simulations were performed in a square box of side L , imposing periodic boundary conditions. A Matérn type II random sampling process was used to identify N c initial cell centres within the box, giving mean cell areaĀ = L 2 /N c , chosen to match A * 6 (given that hexagons are the most frequently observed polygonal class in monolayers (Gibson et al., 2006)). A Voronoi tessellation was constructed between the points (and their periodic extensions) to define an initial network of edges and vertices. The system was then relaxed towards the nearest energy minimum. If the length of any edge fell beneath 0.1 A * 6 (taking the larger value ofÃ * 6 in Region II b ), a T1 transition (or intercalation) was implemented and relaxation proceeded further (see Spencer et al., 2017 for a more refined treatment of this process). If the area of a 3-sided cell fell beneath 0.3A * 6 (again taking the larger value ofÃ * 6 in region II b ), the cell was removed via a T2 transition (extrusion). A small isotropic expansion or contraction of the network and the bounding box was used to satisfy the zero-load condition (3.32) within an prescribed tolerance. The initial disorder produced a distribution of values of P eff α across the cell population.

Results
Simulations for > 0 and < 0 are illustrated in Fig. 5(a-d) respectively. In both examples, the P eff α for individual cells in the disordered monolayer lie close to P eff N , the values for perfect polygons, suggesting that P eff α can be well predicted by a cell's area and its polygonal class. P eff N is monotonic in cell area when < 0 (L 0 > 0), whereas it has a turning point for > 0. Despite the potential for bistability in the latter case, cells in a disordered array lie on both branches of the P eff N curves. In both examples, the mean cell area over the monolayer lies below unity, implying that cells lie below their equilibrium area: each cell is held at this level by cortical tension, as the cell perimeters exceed the target value L 0 . Simulations show that pentagons are smaller on average than heptagons; when < 0 pentagons have P eff α < 0 and heptagons have P eff α > 0 (Fig. 5c); in contrast, for < 0 both sets of cells cluster around P eff α = 0 (Fig. 5a).
The inherent disorder in equilibrium monolayers is illustrated in Fig. 6. The variance of P eff α (about mean zero) within a monolayer of 800 cells is mapped at discrete locations across ( , )-parameter space in Fig. 6(a). For each simulation, L was incrementally adjusted to enforce (3.32). The variability weakens near the region I/II a boundary and increases with . Two individual realizations (Fig. 6b and c) reveal mesoscopic patterns that emerge across the monolayer: shading identifies cells with positive or negative P eff α and line segments characterize the orientation of cell shape and stress. The example closer to = 0 (Fig. 6b) reveals slender patterns that are correlated over many cell lengths. Cells that are larger (smaller) than their equilibrium area, with P eff > 0 (< 0), tend to align with their principal axis of shape (and stress) parallel (perpendicular) to the line of cells, in structures that are reminiscent of force chains in jammed systems (Majmudar & Behringer, 2005). In particular, chains of darker cells are elongated parallel to the chain and exert a net tensile force along each chain, whereas lighter cells are compressed along their chain axis and exert a net compressive force along each chain. Further visualization of these structures is provided in Appendix D (Fig. D1a). In contrast, nearer the Region I boundary (Fig. 6c), the correlation length of patterns increases and there appears to be less alignment of neighbouring cells.  Figure 7 illustrates the impact of varying parameters ( , ) (with P ext = 0) on the shape and size of cells when partitioned into polygonal classes. The mean circularity of cells increases with as one moves across region II a ( Fig. 7a and b): near the region-I boundary, cells with more sides become highly distorted (see inset), whereas near the region IIa/III boundary (where L 0 → 0) cells become more uniformly round. Increasing for fixed near this boundary increases the cortical tension and promotes rounding, while reducing the mean cell area (Fig. 7c and d). Moving back across region II a towards the region-I boundary, L 0 increases, reducing cortical tension and allowing cells to enlarge. In comparison to the size of hexagons, the area distribution across polygonal classes (Fig. 7e) is much more uniform near the region I/II a border than near the II a /III border. The non-linearity in P eff N implies that changes in parameters influence circularity and areas among different polygonal classes non-uniformly. In contrast, the total area occupied by different polygonal classes shows surprisingly little parameter variation (Fig. 7f).
In addition to the model parameters ( , ), the density of cells (controlled by P ext in (3.31)) also induces changes in the equilibrium cell packing configurations. As Fig. 8 illustrates, monolayers under uniform net compression (for which P ext < 0 on average) will tend to produce more round cells, closer to perfect polygons. In contrast, monolayers under uniform net tension (for which P ext > 0 on average) exhibit more disordered arrays, with cells tending to be more elongated. In parameter fitting below, we initially impose the constraint P ext = 0.   6. (a) A map of the variance of P eff α at discrete locations within region II of ( , )-parameter space. Lines show the boundaries for a hexagonal network, as in Fig. 4(a). The dark squares along the region II b /III boundary are artefacts, reflecting the co-existence of cells with small and large areas near this boundary. Each datapoint is taken from 5 realizations of a monolayer with 800 cells. (b) An individual monolayer realization for = −0.1, = 0.1, P ext = 0 with 800 cells. Darker (lighter) shading denotes cells with P eff > 0 (< 0). Line segments indicate the principal axis of the shape and stress tensor for each cell, coincident with the heavy arrows in Fig. 4(b), i.e. aligned with the stress eigenvector associated with the eigenvalue of larger magnitude. (c) A similar example for = −1.11, = 0.15. (Colour in online.)

Parameter fitting
Of the features described in Fig. 7, the total area per polygonal class (panel f) is a poor candidate for parameter identification, while the mean area (panel c) requires a dimensional measure of area and the mean normalized area (panel d) shows limited variation. In contrast, the mean circularity (panel a) shows strong parameter variation without the additional requirement of a lengthscale measurement. However, searching across parameter space we found it difficult to capture simultaneously both the distribution of mean area and the distribution of mean cell circularity. Given the key contribution of cell area to the stress tensor, we therefore chose to use cell area (following Farhadifar et al., 2007) to parameterize the model to the Xenopus laevis animal cap explants introduced in Section 2; we return to circularity below.
Using simulations of monolayers under P ext = 0, we generated datasets A sim ( , ), the mean areas of cells in each polygonal class, to compare with experimental data A exp = {Ā Evaluating (4.1) across a grid of parameter samples in region II (Fig. 9a), the posterior was maximized with ( , ) ≈ (−0.26, 0.17), for which L 0 ≈ 0.76. While there are other credible parameter regions near the region III boundary, we can be confident that the monolayer in this experiment is far from the rigidity transition at region I, and reasonably certain that it falls outside region II b (where L 0 < 0). The distribution of area across polygonal classes is captured well by the model (Fig. 1e). For best-fit parameters, cells i18 A. NESTOR-BERGMANN ET AL. which are larger than average (shaded dark in Fig. 9b) tend to align in slender structures or, in some instances, to be isolated at the centre of a rosette of smaller (pale) cells.
Despite matching area distributions well, the circularity distribution is over-estimated across all polygonal classes (Fig. 1f). Figure 8(b) suggests that the circularity can be reduced by putting the monolayer under net tension. To investigate the possibility that the thin basal tissue layer of the animal cap ( Fig. 1a and  b) might induce such a tension in the apical epithelium, we ran additional simulations for which P ext > 0 (see (3.31)), maintaining fixed values of † and † (see (3.40)) in order to remain in an equivalent region of parameter space. A demonstration of the changes in cell area and circularity across polygonal classes as P ext for ( † , † ) = (−0.259, 0.172) is given in Fig. 9(c and d). While the area distribution maintains close agreement with experiment as P ext increases, the circularity moves towards the experimental range but does not fall comfortably within it, even for very large P ext . We conclude that additional refinements to the model (such as higher order non-linearities in the energyŨ α , see (3.3)) may be necessary to ensure quantitative agreement of both area and circularity distributions.

Discussion
We have investigated a popular vertex-based model of planar epithelia, addressing features associated with cell packing rather than division or motility. We focused on a simple version of the model, neglecting refinements such as representations of internal viscous forces (Okuda et al., 2015), non-planarity (Hannezo et al., 2014;Murisic et al., 2015;Bielmeier et al., 2016), descriptions of curved cell edges (Brodland et al., 2014;Ishimoto & Morishita, 2014), internal anisotropy, multiple cell types and so on. We first derived an expression (3.23) for the stress σ α of an individual cell, expressed in terms of its shape. The isotropic component of stress reveals the cell's effective pressure P eff α (3.24), which is set by a balance between the internal pressure associated with bulk (cytoplasmic) forces that regulate cell area and cortical tension that regulates the cell perimeter. With the area below and the perimeter above their respective targets (A α < 1 and L α > L 0 ), the bulk forces push outward against the stretched perimeter, giving the cell some rigidity. The traceless tensor J α in (3.23) characterizes asymmetries in the cell shape that might arise from an imposed shear stress or, in the absence of an external load, internal asymmetries associated with intrinsic disorder. A simple representation of viscous forces associated with drag from the underlying substrate leads to a further contribution to the stress associated with dynamic shape changes. Crucially, the principal axes of the shape tensor S α (defined in terms of the vertex locations) align exactly with the principal axes of the cell stress, as illustrated in Fig. 3. This result may have implications in cell division, where it is postulated that there may be shape-and stress-sensing mechanisms guiding the positioning of the mitotic spindle (Théry & Bornens, 2006;Minc et al., 2011). If the vertex-based model is accepted as a leading-order description of cell mechanics, it follows that it will not be possible to separate these mechanisms by looking solely at cell geometry, since the orientation of any inferred stress will necessarily align with the cell shape. Instead, the system must be perturbed, either mechanically or biochemically (using biological knockdowns, for example), such that the mechanisms can be disrupted and separated. In this context, it is worth highlighting the distinction between the orientation of external stress that may be imposed on a monolayer, and the heterogeneous stress field at the individual cell level (e.g. Fig. 3). Observations show cell division in a stretched monolayer to be aligned with cell shape rather than the i20 A. NESTOR-BERGMANN ET AL. external stress orientation (Wyatt et al., 2015); the present model suggests that the cell-scale stress would be aligned with cell shape, even if the average stress at monolayer level has a different orientation.
The distinction between individual cell stress σ α and tissue-level stress σ R is evident in the expression (3.26) for the stress over a patch of cells, derived as an area-weighted average of the individual cell stresses. For a monolayer under an isotropic external load of magnitude P ext , we derived a constraint (3.31) on the area-weighted P eff α ; furthermore, the averaged deviatoric stress must vanish in this case. When simulating a monolayer that is not subject to lateral forcing, the constraint of zero mean effective pressure (3.32) is important in determining the appropriate cell density within the simulation domain. One can then examine the properties of the monolayer when this configuration is perturbed by small compressive or shear deformations. We derived an exact expression (3.33) for the monolayer's bulk elastic modulus (generalizing results obtained previously for hexagonal cell arrays) and recovered directly an expression (3.35) for the shear modulus in the hexagonal packing limit. The mechanical properties of the tissue can therefore be tuned by varying the relative strengths of the bulk and cortical forces. As shown previously (Bi et al., 2015), a phase transition arises when L 0 ≈ 3.81, which bounds a region of parameter space in which the monolayer loses resistance to shear deformations. Fitting our model to data from embryonic Xenopus laevis tissue, by maximizing a likelihood function derived from the mean area per polygonal class, suggests L 0 ≈ 0.76 in the embryonic tissue, substantially distant from the rigidity transition. The model fit is imperfect however, as we were not able to capture circularity distributions even when varying the peripheral load on the monolayer (Figs 1f and 9). This suggests further constitutive refinements of the model are needed, such as including higher-order non-linearities in (3.3). We also examined how cell shape (and of course size) can be influenced by an external load P ext , with cells becoming rounder when tightly packed (Fig. 8). The bulk isotropic stress (or equivalently the mean cell density) is likely to be a significant parameter when simulating confined tissues, and is an example of a mechanical signal that can be communicated over long distances. Future studies should address anisotropic external loading, which has the capacity to promote more ordered cell packing (Sugimura & Ishihara, 2013).
The present descriptions of the stress tensor are appropriate for small-amplitude deformations close to equilibria, and in future should be extended to account for irreversible cell rearrangements (such as T1/T2 transitions) that endow the material with an elastic-viscoplastic character, as well as accounting for cell division. Kinematic and geometric quantities (such as the texture tensor) characterizing large deformations of cellular materials have been developed that are based on connections between centres of adjacent cells (Graner et al., 2008;Blanchard et al., 2009;Etournay et al., 2015;Guirao et al., 2015;Tlili et al., 2015;Blanchard, 2017), the dual network to that considered here. While it is straightforward to repartition the stress (3.26) over the network of triangles connecting cell centres, it is less clear how to relate it to strain measures defined with respect to cell centres rather than cell vertices, without for example assuming that vertices are barycentric with respect to cell centres (Barton et al., 2016). In particular, the relationship between the tissue-level stress postulated by Etournay et al. (2015) to that emerging from the vertex-based model remains to be established.
While the monolayer can be stress-free at the bulk scale, individual cells can have non-zero P eff α : those for which P eff α > 0 (< 0) are larger (smaller) than the equilibrium area at which bulk and cortical forces balance. Each simulation of a spatially disordered monolayer describes an equilibrium configuration of this very high-dimensional dynamical system, subject to the constraint that all edge lengths exceed a defined threshold (smaller edges being removed by T1/T2 transitions). We have characterized some features of the variability of these states, both in terms of the variance in P eff α over the cell population and the spatial pattern of compressed and dilated cells. While soft monolayers near the region I/II boundary show very long-range patterning (Fig. 6c), stiffer monolayers nearer the II/III boundary appear to exhibit chains of force (and cell shape, Figs 6b and D1a), where lines of tension and compression are transmitted along i21 entangled strings. Evidence of force chains has recently been provided in the Drosophila melanogaster embryo (Gao et al., 2016) and the patterns suggested by our model (Fig. 9) motivate ongoing investigations in the Xenopus system. Robust evidence of force-shape chains in real epithelia would raise interesting questions about the role of mechanical feedback on patterning of cell division.

Appendix A. Geometric identities
The contribution to the stress due to cell pressure first involves Taking components, giving (3.22a). Referring now to the contractility term, we find giving (3.22b). This is symmetric becauset i α ⊗ t i α = t i α ⊗t i α .

Appendix B. Proof that S α ,Ṡ α and J α align
To establish that S α J α = J α S α for cell α, we can ignore the pre-factors in the tensors and need only show Let us henceforth assume that the sums over i, j and q are implicit. We also drop the α subscripts, under the assumption that all vectors are relative to the same cell centroid. Considering the left hand side (LHS) first: matching (B.2). Therefore the tensors commute and we have alignment of the principal axes of stress and shape, when the system is in equilibrium.
Let us now establishṠ α S α = S αṠα . Ignoring pre-factors again, we havė which is symmetric. Given thatṠ α , S α are both symmetric, and their product is symmetric, we have a necessary and sufficient condition that they commute. We therefore also have alignment of the principal axes of stress and shape when the system is out of equilibrium.