Force networks, torque balance and Airy stress in the planar vertex model of a confluent epithelium

The vertex model is a popular framework for modelling tightly packed biological cells, such as confluent epithelia. Cells are described by convex polygons tiling the plane and their equilibrium is found by minimizing a global mechanical energy, with vertex locations treated as degrees of freedom. Drawing on analogies with granular materials, we describe the force network for a localized monolayer and derive the corresponding discrete Airy stress function, expressed for each N-sided cell as N scalars defined over kites covering the cell. We show how a torque balance (commonly overlooked in implementations of the vertex model) requires each internal vertex to lie at the orthocentre of the triangle formed by neighbouring edge centroids. Torque balance also places a geometric constraint on the stress in the neighbourhood of cellular trijunctions, and requires cell edges to be orthogonal to the links of a dual network that connect neighbouring cell centres and thereby triangulate the monolayer. We show how the Airy stress function depends on cell shape when a standard energy functional is adopted, and discuss implications for computational implementations of the model.

OEJ, 0000-0003-0172-6578; SLW, 0000-0001-5680-4030 The vertex model is a popular framework for modelling tightly packed biological cells, such as confluent epithelia. Cells are described by convex polygons tiling the plane and their equilibrium is found by minimizing a global mechanical energy, with vertex locations treated as degrees of freedom. Drawing on analogies with granular materials, we describe the force network for a localized monolayer and derive the corresponding discrete Airy stress function, expressed for each N-sided cell as N scalars defined over kites covering the cell. We show how a torque balance (commonly overlooked in implementations of the vertex model) requires each internal vertex to lie at the orthocentre of the triangle formed by neighbouring edge centroids. Torque balance also places a geometric constraint on the stress in the neighbourhood of cellular trijunctions, and requires cell edges to be orthogonal to the links of a dual network that connect neighbouring cell centres and thereby triangulate the monolayer. We show how the Airy stress function depends on cell shape when a standard energy functional is adopted, and discuss implications for computational implementations of the model.

Introduction
Multicellular biological tissues have an intrinsically granular structure, associated with the mechanical integrity of individual cells. While cells may be sufficiently soft for many tissues to deform like continuous media described by smooth strain fields [1], stress fields can remain heterogeneous [2] and may display features that are not captured in terms of smoothly varying (homogenized) variables. Accordingly, the vertex model of tightly packed cells [3][4][5][6][7][8][9][10][11][12] has become a popular framework with which to model plant and animal development, particularly of tightly packed epithelial monolayers. The vertex model captures cell geometries efficiently, enables straightforward computation that resolves individual cells, and is based on simple mechanical assumptions. Integrating over regions, it can be used to derive tissue-scale properties such as elastic moduli [13][14][15]. In addition to capturing a jamming/unjamming phase transition, with resistance to shear vanishing as cells lose cortical tension-a topic of much current attention [13,[16][17][18][19]]-the vertex model also exhibits inherently discrete mechanical structures (such as force chains and correlated patterns of stress [20,21]), which have the potential to influence biological behaviour. Despite its popularity, however, the mechanical constraints underpinning the vertex model have not yet been fully articulated.
In classical elasticity, materials are defined with respect to a reference state, using a strain energy function defined in terms of strain invariants. The vertex model differs in using cell area and perimeter as intrinsic measures of shape (for systems such as epithelia that are well described by two-dimensional models), and the concept of a reference state is not employed. In many ways, the manner in which cells pack together under an external load instead resembles a granular material, which can accommodate multiple configurations under given boundary conditions [22]. Here, we exploit this analogy to identify the force network associated with a planar cell configuration, and derive the corresponding force potential and Airy stress function. We show that the Airy stress function is defined over kites that tile individual polygonal cells. Whereas stress components can be expressed as second derivatives of the Airy stress function in a planar elastic material, here stress is constructed using discrete derivatives, as deployed for granular media [23][24][25] and in models for self-equilibrated frameworks [26]. Accordingly, we exploit some machinery from graph theory and discrete calculus, making extensive use of incidence matrices, which serve as analogues of finite-difference (coboundary) operators (or, when transposed, as boundary operators) [27][28][29][30], while avoiding the full formalism of exterior calculus.
The Airy stress function serves as a discrete scalar potential for the vector force potential, and its existence guarantees that intra-and intercellular stress tensors are symmetric, i.e. that there is a torque (or moment) balance across a monolayer. We show in the present case that this condition places a geometric constraint on the intercellular stress in the neighbourhood of cellular trijuctions. This stress-geometry condition is provided by a fabric tensor resembling that described by Ball & Blumenfeld [31] for granular materials; to our knowledge it has not been used previously in the context of the vertex model. We show how the fabric tensor can be used to determine the orientation of stress in the neighbourhood of trijuctions. Furthermore, we show that a torque balance in intercellular stress (not normally considered in biological studies that focus on intracellular stress, nor imposed in simulations that only apply a point-wise force balance on vertices) reveals the requirement that links between cell centres (appropriately defined) should, within the framework of the vertex model, be orthogonal to the cell edges that they intersect and, crucially, that each cell vertex should lie at the orthocentre of the triangle connecting adjacent edge centroids. We show how these constraints can be used to identify a consistent triangulation of the monolayer that is dual to the network of cell boundaries.
The vertex model is of course a simple idealization of a complex biological system. The geometry of a typical epithelium (figure 1a) is summarized by the locations of its trijunctions (vertices), combined with topological information identifying the cell edges connecting vertices, and the cells that are bounded by edges ( figure 1b)  a degree of intrinsic disorder, captured for example by a distribution of edge lengths (figure 1c). Figure 1b illustrates one possible dual network, constructed in this instance by links connecting the centroids (defined with respect to cell vertices) of adjacent cells. The links also show variability in length (figure 1b). The angles at which links intersect their corresponding cell edges are quite tightly distributed around π/2 (figure 1d), but show some evidence of non-orthogonality. We discuss this observation in light of theoretical predictions below.
In this study, we ignore neighbour exchanges (T1 transitions), cell extrusion, cell division and intrinsic cell motility, focusing simply on monolayer configurations with fixed topology. For simplicity, we also assume that all internal vertices are trijunctions. In §2, we implement the planar vertex model using incidence matrices and lay out some relevant geometric and topological results before representating intra-and intercellular stress fields in terms of potentials in §3. These results are intrinsic to the vertex-based description and independent of a constitutive model, which we introduce in §4. Adopting a widely used approximation for cell elastic energy, we show how intracellular variations in Airy stress function are proportional to the cell's cortical tension, and can be expressed directly in terms of cell shape. Findings are summarized in §5, where we propose a potential computational strategy that respects torque balance and discuss the relevance of the model to real epithelia.

The planar vertex model
We consider a localized monolayer of N c confluent cells, represented as tightly packed polygons covering a simply connected region of the plane. We assume that an external isotropic stress P ext is applied around the periphery of the monolayer. In computations, starting from some (typically disordered) initial condition, vertex locations either evolve under a local force balance until the system reaches equilibrium, or they are adjusted directly to minimize a global energy. In either case, each vertex in the monolayer can be assumed instantaneously to be under zero net force (inertial effects are neglected). We wish to understand the impact of imposing, additionally, a torque balance across the monolayer.

(a) Cell topology and geometry
Given the nature of the vertex model, and the quality of available imaging data, we take cell boundaries as the primal network, which we assume is embedded in a Euclidean space. The cellular monolayer is, therefore, defined by a set of vertices (position vectors) r k ∈ R 2 , k = 1, . . . , N v , a set of oriented cell edges t j (of length t j ), j = 1, . . . , N e and a set of oriented cell faces (that we simply call cells) a i (of area A i ), 1 = 1, . . . , N c . Here a i = A i i where i = ±ε represents a clockwise rotation by ±π/2. (ε is the 2D Levi-Civita symbol satisfying ε T = −ε, εε = −I; the summation convention is not adopted here.). Orientations of edges and faces are prescribed but arbitrary; here we will assume that all cells have the same orientation. We collect vertices, edges and faces into vectors r = (r 1 , . . . , r N v ) T , t = (t 1 , . . . , t N e ) T and a = (a 1 , . . . , a N c ) T but for clarity use matrix notation sparingly below, writing sums explicitly in many cases.
The topology of the monolayer is defined using two incidence matrices [28]. The N e × N v matrix A has elements A jk that equal 1 (or −1) when edge j is oriented into (or out of) vertex k, and zero otherwise. The N c × N e matrix B has elements B ij that are non-zero only when edge j is on the boundary of cell i, taking values +1 if the edge is coherent with the orientation of the cell face and −1 if not. Replacing −1 with 1 in each matrix produces unsigned incidence matrices A and B, identifying neighbours but not orientations. Further properties of A and B are given in appendix A. The N c × N v matrix C = 1 2 B A has elements C ik that equal 1 if vertex k neighbours cell i and zero otherwise. Thus Z i ≡ k C ik (summing over all vertices) defines the number of edges (and vertices) of cell i. We let R i represent the centre of each cell, without specifying yet how it might be related to the cell's vertex locations ∪ k C ik r k (where ∪ k denotes collection, without summation, over all vertices).
To account for boundaries of the monolayer, vertices (and all other functions defined on vertices, with subscript k) are partitioned as N p peripheral and N v − N p interior vertices so that r = [r p , r i ] T , edges (and relevant functions with subscript j) as N p peripheral, N b border and N e − N p − N b interior edges so that t = [t p , t b , t i ] T , and cells (and functions with subscript i) as N b border and N c − N b interior cells so that a = [a b , a i ] T . A peripheral edge has two peripheral vertices; a border edge has one peripheral and one interior vertex; an interior cell has only interior edges. Internal vertices always represent trijunctions. Figure 2a illustrates this for a small monolayer of seven cells. We may then partition the incidence matrices as where A pp is an N p × N p matrix, etc., so that Edges are defined by t j = k A jk r k , with lengths t j = t j · t j . This defines the unit vectorst j = t j /t j . The perimeter of cell i is L i = j B ij t j (summing over all edges). It follows (for later reference) royalsocietypublishing.org/journal/rspa Proc. R. Soc. A 476: ∂L i /∂r k is therefore the sum of two unit vectors aligned with the two edges of cell i that meet vertex k, pointing into the vertex.
To define cell areas, we construct n ij defines the outward normal of cell i at edge j and c j defines the centroid of edge j. Let φ = 1 2 x · x, where x is a position vector in R 2 and integrate (∇ ⊗ ∇)φ = ∇ ⊗ x = I over cell i, where ⊗ denotes the dyadic outer product. Applying the divergence theorem to an integral over cell i, The oriented area of cell i can therefore be written as The trace of (2.6) gives j B ij t j · c j = 0. This can be understood by recognizing φ as the potential for position x = ∇φ; its discrete form φ k = φ(r k ) jumps by k A jk φ k = c j · t j along edge j, and the net change in potential vanishes around a closed loop because BA = 0 (appendix A), a device we will exploit later on. Also, as shown elsewhere (e.g. [19,21]), ∂A i /∂r k is, therefore, the sum of two inward normal vectors associated with the edges of cell i meeting at vertex k, with length equal to half of each edge.

(b) Dual networks
There are multiple networks that are dual to the (primal) network of cells. The simplest is the triangulation (a simplicial complex) connecting adjacent cell centres. Assigning orientation k to all triangles (opposite to that in all cells), the orientations of links between cell centres are induced by the choice of A and B (appendix A), with link T j = i B ij R i dual to edge t j . For a localized monolayer, peripheral triangles and links are truncated; complete links are given by We will also make use of a second dual network, formed by links between cell centres and edge centroids (figure 2a). This partitions each cell into kites (described in more detail below), with three kites surrounding each vertex. The resulting six-sided tristar at each vertex shares three vertices with the triangle connecting cell centoids, but their edges in general are distinct. We denote the area of the tristar at vertex k as E k .
A more fine-grained edge-centroid network is built by connecting adjacent edge centroids around each cell. Thus defines links between adjacent edge centroids running clockwise as polygons around cells and anticlockwise as triangles around vertices (figure 2a; appendix B). To invert (2.8), we may use where P(j, j ) denotes the set of paths over the edge-centroid network connecting j and j , demonstrating how c j is a discrete vector potential for s ik . As loops around any interior vertex k or any cell i are closed, it follows that to give tr(C bpT S bp ) = 0, because the boundary of the centroid network is closed, while diagonal elements of C biT S bi + C iiT S ii vanish (representing closed loops around interior vertices); all diagonal elements of CS T vanish (representing closed loops around cells). Finally, dual to the edge-centroid network is the network of spokes connecting cell centres to vertices. The outward radial spokes of cell i satisfy C ik (q ik − r k + R i ) = 0 (figure 2c).

(c) Kites
We combine spokes and links between edge-centroids to build kites (figure 2a,c). The links between the cell centre and the edge centroids defining the boundaries of kite ik within cell i are where u ik = j B ik c j A jk is the sum of the two edge centroids bounding kite ik, so that c 3ik and c 4ik run anticlockwise around the kite (figure 2c). The area of the kite is given by q ik ⊗ s ik ) (see appendix C). Following [31], we can write where c 1ik and c 2ik run anticlockwise along cell edges. The area of tristar k is therefore E k ≡ i C ik K ik . Summing kites over the tristar, the internal edge contributions (involving c 1ik and c 2ik ) cancel leaving only boundary contributions, giving (2.13) The fabric tensor F k measures the asymmetry of each tristar [31]; it can be written (appendix B) as Constructing a cell from kites, edge contributions cancel as well (because kites are defined on edge centroids), giving an alternative formulation of the cell area as

Representations of cell and tissue stress
Let f ik be the force on vertex k due to cell i. The requirement that the net force at interior vertex k and the net force on any cell i both vanish is representing two discrete divergences of f ik . Stating (3.1) more generally to account for boundary forcing, we require the diagonal entries of C T (F − F ext ) to vanish (balancing forces at each vertex, including the periphery), and the diagonal entries of CF T to vanish (an internal force balance on each cell), where the matrices F and F ext share the structure of C in (2.2) and {F} ik ≡ f ik . Now summing over cells and vertices, respectively, and the external force (imposed pressure around the monolayer periphery) has matrix blocks . Thus zero diagonal entries of C biT F bi + C iiT F ii give (triangular) force balances at interior vertices (including contributions from peripheral cells where appropriate). Zero diagonal entries of CF T give (polygonal) force balances over interior and peripheral cells, and zero diagonals of C bpT (F bp − F bp ext ) give the force balance on peripheral vertices. For a monolayer satisfying (3.1), the first moment of the force defines the stress σ i associated with cell i via We call the isotropic component of the stress in each cell the effective pressure, P eff,i ≡ 1 2 tr(σ i ). The stress σ of the monolayer as a whole may then be written as  (3.6) showing that the total stress must be isotropic, internal shear stresses must cancel and therefore that [21] i A i P eff,i = AP ext . (3.7) Equation (3.6) also ensures zero net torque on the monolayer due to P ext . We now consider how the force balances (3.1) can be represented geometrically, with a view to identifying the (intercellular) stress σ k defined over tristars.

(a) The force network
The connection between the force network and the edge centroid network becomes clear if we rotate each force anticlockwise by π/2 (via a Maxwell-Cremona construction [22]): then i C ik (−εf ik ) = 0 and k C ik (−εf ik ) = 0, implying that the rotated force vectors form a network that is topologically equivalent to the edge-centroid network (2.10), with closed triangles around vertices and closed polygons around each cell (figure 2b). While the edge-centroid network is planar (by construction), the force network may not be. In particular, the peripheral forces (3.3) map to which collectively form a closed loop, matching the shape of the perimeter of the edge-centroid network (connecting all the peripheral centroids). Fixing the location of one peripheral edge centroid at the origin, the loop is clockwise if P ext > 0, anticlockwise if P ext < 0 and collapses onto the origin if P ext = 0. The centroids c j form a discrete potential for the edges s ij via (2.2), (2.9). Similarly, we can identify the vertices of the force network h j (figure 2b) as a potential for the forces, by writing The stress over cell i can then be written in terms of the force potential as noting that C ik becomes redundant when A jk and B ij both appear in the sum, and using t j = k A jk r k . Taking a transpose gives σ i should be symmetric for cell i to be under zero torque. This requires and allows us to write the cell stress as a discrete curl of h around its periphery via , giving the stress in rotated basis. It follows that (3.14) Comparison with A i I = j n ij ⊗ c j (from (2.6)) suggests that stress can be associated with a mapping from c j to h j . Equation (3.14) also shows that the isotropic component of cell stress is a discrete divergence,  Figure 3. Three cells, labelled i, i and i , sharing vertex k and edges j, j and j (taken anticlockwise). In cell i, spoke q ik connects cell centre R i to vertex r k , intersecting the link s ik between neighbouring edge centroids c j and c j . Kite ik is spanned by q ik and s ik ; the three kites neighbouring vertex k, that together form a tristar, are shaded. The Airy stress function ψ ik (see below) is defined on kites. Jumps in ψ ik between neighbouring kites in the same cell, sharing vertex c j , are defined by the projection of the force potential h j on the cell edge t j . Jumps in ψ ik between neighbouring kites in the same tristar, sharing vertex c j , are defined by the projection of h j on the link T j between cell centres (not shown) that intersects edge t j . (Online version in colour.)

(b) Stress as a map between networks
We can compare (2.15), which constructs cell area from kite areas, to stress written as (3.4), suggesting that stress can also be interpreted as a mapping between s ik and −εf ik . An explicit construction for such a map was provided in [31]. The mapping M k between vertices c j , c j , c j , running anticlockwise around a triangle surrounding vertex k (figure 3), to h j , h j , h j , also ordered anticlockwise, is (3.16) where the triangle area a k satisfies 2a With the map defined, we can write f ik = εM k s ik . Then from (3.4), This shows how the cell's stress is built from the kite shape tensor q ik ⊗ s ik in (2.13), weighted by contributions −M T k ε from the cell's vertices. Accordingly, the stress over a tristar built from the three kites surrounding vertex k is (3.19) where the fabric tensor F k is given by (2.14). The stress over the tristar at vertex k can also be written in terms of spokes q ik and vertex forces f ik . Replacing the latter with the force potential h j and reordering, we find that multiplying each h j is the difference between neighbouring spokes, i.e. the straight link T j between cell centres. In general, this does not pass through c j , so contributing to non-zero F k . Explicitly, using (3.9), (3.20) where T j = i B ij R i and k is the orientation of the triangle surrounding vertex k. (As explained in appendix A, we assume that all cells have the same orientiation i , and that k = − i uniformly; with fixed j and k, B ij + B i j = 0 for all options in figure 5 below, so that the r k terms cancel in (3.20).) The outward normals to triangle k are N jk = − k A jk T j . It follows that, analogously to (3.15), the isotropic component of tristar stress is (c) Expressing stress in terms of the Airy stress function To enforce zero torque on cell i, given by (3.12), we define a discrete potential ψ ik (the discrete Airy stress function, which assigns a scalar value to each kite of the monolayer) satisfying for either of the cells neighbouring edge j, i.e. with B ij = 1, which automatically satisfies j,k B ij A jk ψ ik = 0 (because BA = 0-see appendix A). Likewise, zero torque on tristar k requires, from (3.20), j A jk h j · T j = 0, which we satisfy with potential ψ ik satisfying for both pairs of kites bounding edge j (i.e. with A jk = 1). Pursuing the analogy with planar elasticity, we seek to define h j as a discrete curl of ψ ik , here evaluated over the spokes q ik of the four kites surrounding c j (e.g. the path q ik − q i k + q i k − q ik in figure 3) in order to recover the stress in terms of ψ ik .
To illustrate the definition of ψ ik , consider the three kites surrounding vertex k (figure 3), noting that links between neighbouring cells can be expressed in terms of spokes. With i, i , i and j, j and j ordered anticlockwise around vertex k, we require from (3.22) that Likewise, at neighbouring vertex k , we require showing that the jump in ψ across edge j is symmetric between neighbouring kites. Accordingly, we can define averages of ψ ik over neighbouring elements, φ ij = 1 2 k A jk ψ ik and θ jk = i We now express h j in terms of the ψ ik in the four neighbouring kites (i.e. inverting expressions such as (3.23)), using the network of spokes. Equation (C 3) (appendix C) demonstrates how a vector g can be constructed as a discrete curl of a potential defined across a diamond spanned by non- parallel vectors a and b. Assuming there are two jumps in potential, when crossing a and  respectively, with the jumps proportional to g · a and g · b but not a linear combination of the two (as is the case for h j , t j and T j in (3.22)), it is necessary for a · b = 0. We therefore require t j · T j = 0, (3.26) i.e. each link between adjacent cell centres must intersect the corresponding cell edge orthogonally. Equation (3.26) is therefore necessary for both σ i and σ k to be symmetric (equivalently, for each to be expressed in terms of ψ ik ). For the jumps in h j · t j and h j · T j to align appropriately with t j and T j , a rotation and rescaling are necessary as in (C 4), to give Given (3.26), we can also express the force potential directly in terms of edges and links as Recalling that i t j defines a normal to edge j relative to cell i, we see that noting that, for all four cases in figure 5 below, ( i t j ) · T j = t j T j . Thus from (4.3), and using (3.25), we obtain an alternative to (3.15) As might be expected from classical elasticity, the isotropic component of the stress is given as a discrete Laplacian (over the primary network) of the Airy stress function, involving (for an interior cell) 3Z i kites and 2Z i independent values of ψ ik . Likewise, noting that ( k t j ) · T j = t j T j , the isotropic stress over tristars is given by a Laplacian over the dual network involving (for an interior cell) 9 kites and 6 independent values of ψ ik , namely providing an alternative to (3.21). Finally, we can write the tristar stress in terms of links and edges using (3.20, 3.28) as in each of the four cases illustrated in figure 5 below. Thus, making use of (3.25) and (3.31), showing how the shear stress is captured by differences in the ψ ik field between neighbouring kites intersecting the tristar. Likewise, using the identities t T j i = −T T j (t j /T j ) and T T j i = t T j (T j /t j ), we find products of the unit tangent and the inward unit normal, weighted by v ij and taken anticlockwise, such that j v ij = 0. It is not immediately obvious that the stress tensors in (3.33, 3.34) are still symmetric. However, writing the outer product of unit vectors ast j ⊗T j = (cos α j , sin α j ) T (− sin α j , cos α j ) (when i = −ε), where α j is the orientation of edge j with respect to a fixed axis, then the final sum in (3.34) is (for i = ±ε) which we will make use of shortly.

(d) Tristar stress and the fabric tensor
We now reconcile the two expressions for σ k in (3.19) and (3.20). First, consider the condition Link T j crosses edge t j , bounded by two vertices, each surrounded by a triangle of vectors s ik of the edge centroid network having area a k ; two such triangles are illustrated in figure 3. Depending on the chosen orientation of cells and edges, (3.37) implies that link T j is parallel (or antiparallel) to the furthest edge of each triangle (such as s i k in figure 3), with the magnitude of T j relative to the edge s ik given by the ratio E k /2a k . In other words, (3.37) implies that each vertex bounding t j lies at the orthocentre of the triangle of edge-centroid-links surrounding each vertex. Direct substitution of (3.37) into (3.20), giving σ k as an outer product of links with the force potential, recovers −εM T k ε with M k defined in (3.16), given as an outer product of edge-centroidlinks with the force potential. Thus (3.37) is equivalent to the condition σ k = −εM T k ε. (3.38) Symmetry of σ k is ensured by the existence of the Airy stress function and the orthogonality condition (3.26); (3.37) extends this symmetry to M k . Furthermore, (3.19) then implies that F k M k = 0, while (3.38) gives M k = −εσ k ε, yielding the stress-geometry condition [24,31,32] F k εσ k = 0. The role of stress as a mapping between networks is also evident via f ik = σ k εs ik , showing how a force balance can be turned into a divergence of stress (via (3.1)). The mapping can also be used to show that the area of the triangle k in the force network is det(σ k ) times that in the edge centroid network (appendix D). We can use (3.39) to infer stress orientation in the neighbourhood of a vertex. As long as det(σ k ) = 0, we can write σ k in terms of its principal axes and eigenvalues as σ k = σ k1 e k1 ⊗ e k1 + σ k2 e k2 ⊗ e k2 , where e k1 · e k2 = 0. Likewise, as long as det(F k ) = 0, we may express it in terms of its principal axes as This will be satisfied when the orthogonal axes of each tensor align, so that the cross products vanish. The fabric tensor, therefore, provides a direct mechanism for inferring stress orientation in the neighbourhood of vertices, except when there is sufficient symmetry for the fabric tensor to vanish.
The R i correspond to 2N c scalar quantities. Given a set of vertices, for 2N c > N e − N p the system is underconstrained and one expects to find many possible cell centre locations for which (3.26) is satisfied (i.e., for a small number of cells, it is easy to construct a triangulation of cell centres for which links are orthogonal to edges). However, for larger monolayers, with 2N c < N e − N p (anticipating that N e ∼ 3N c and N p ∝ √ N c for N c 1), then (3.41) becomes overconstrained. In other words, a set of vertex locations emerging from a simulation that does not impose a moment balance cannot be expected, in general, to admit a triangulation satisfying (3.41). Similarly, figure 1d shows how (3.41) is violated when cell centres are chosen to be cell centroids, satisfying The constraint (3.37) can be interpreted from a cell-based perspective (in which cell vertices r define cell centres R), as illustrated in figure 4a. The vertices ∪ k C ik r k of cell i define its edges ∪ j B ij t j , edge centroids ∪ j B ij c j and the internal links between them ∪ k C ik s ik . Equation (3.37) requires that the edge radiating outwards from cell i at vertex k must be orthogonal to s ik ; the triangle of links ∪ i C ik s ik around vertex k is then fully defined since r k is at its orthocentre. This in turn specifies the centroid (and therefore length) of each edge radiating from the cell. A triangulation of the plane is then constructed (in principle) via each triangle of the edge-centroid network at each vertex being rotated by π and expanded by a factor 1/λ k , say, under the constraint (3.37), to cover the plane (figure 4a), with R i the common vertex of all triangles covering cell i. Noting that 2a k is the product of the base and the height of the edge-centroid triangle at vertex k, link j should have length T j satisfying (from (3.37)) where h ik is the altitude of the triangle at vertex k (figure 4a). Given that the area U k of the triangle spanned by ∪ i C ik R i maps to a k via a k = λ 2 k U k , with altitudes related by a result that can be verified by (2.13).
The covering of the plane by rotated and expanded triangles requires the internal angles at the vertices of the triangles to sum to 2π where they meet at R i . (To demonstrate that this is feasible, consider the closed triangles linking edge centroids around the vertices of cell i. The outermost vertex of each triangle is bisected by an edge separating two neighbouring cells. The resulting angle α, marked in figure 4a appears also within another vertex of the triangle, indicating that π/2 − α contributes to the internal angle of the polygon of links s ik within cell i. Since all such internal angles sum to (Z i − 2)π , it follows that all angles α sum to 2π .) The covering also requires consistent scaling between neighbouring triangles. We show in appendix F how this condition is satisfied (up to a translation and uniform scaling of the triangulation), and how such a covering may extend to the whole monolayer. Figure 4b illustrates such a covering for three cells, satisfying the orthocentric property (3.26), (3.37). In summary, the zero-net-torque constraints described above, specifically the requirement that each internal vertex lies at the orthocentre of the triangle formed by its neighbours (or equivalently its neighbouring edge centroids), can be used to define a self-consistent set of cell centre locations.

Constitutive modelling
So far, we have assumed that the mechanical load applied to (or generated by) a cell can be approximated by forces applied at its vertices, without specifying how these might be related to the size and shape of the cell. Commonly, cell i, with area A i and perimeter L i , is assumed to have mechanical energy U i = U(A i , L i ), so that cells have identical mechanical properties but distinct shapes and sizes. U typically includes a quadratic area-dependent term penalizing departures from a reference area, that measures the resistance of the cytoplasm to expansion or contraction, and a quadratic perimeter-dependent term that penalizes departures from a reference length, capturing the resistance to stretching of the cell cortex as may take place under shear. These contributions define a pressure P i and a tension T i for each cell via Equations (2.3) and (2.7), showing how the length and perimeter of cell i change when vertex k moves, can then be used to evaluate f ik = δU i /δr k , the first variation of the energy of cell i with respect to a small displacement of vertex r k . This determines the restoring force at r k due to cell i as We can evaluate the cell stress using (2.6) and (4.1), noting that C ik becomes redundant, as j B ij t jtj ⊗t j , a result consistent with prior studies (e.g. [8,19,21,33]). These terms can be interpreted by noting that under an imposed uniform strain E, A i changes by A i I : E ≡ A i Tr(E) and L i changes by L i Q i : E [15]. The cell structure tensor Q i commutes with the cell shape tensor k q ik ⊗ q ik , implying that the principal axes of stress and shape align at the cell level [21]. The isotropic component of the stress in each cell shows that the effective pressure has contributions from both the cell's interior and its periphery. (3.34) shows that the intracellular differences in ψ ik are proportional to T i L i . More specifically, writingt j = (cos α j , sin α j ) T with respect to some fixed Cartesian axes, recalling that L i = j B ij t j . Comparison with (3.36) then enables us to express the intracellular jump in Airy stress function between kites explicitly, as given in appendix E, revealing that the variation of Airy stress function within a cell is T i L i times a nonlinear but dimensionless measure of cell shape. We note that, since U is defined in terms of areas and perimeters, it clearly satisfies material frame indifference. Under a change in reference frame, in which r → Yr + Z, where Y is a rotation and Z a translation, it is straightforward to show that t → Yt (as do other vectors contributing to f in (4.1)), while the identity ε = YεY T ensures that other tensors, such as stress, transform under σ → Yσ Y T , and hence are also frame indifferent.
The constitutive model can also be extended to accommodate viscous dissipation, either within the cell itself or as a result of substrate drag. In the former case, at time t, P i in (4.1) is replaced with with P i + γ dA i /dt, and T i with T i + μdL i /dt, corresponding to a dissipation rate Φ i ≡ γ (dA i /dt) 2 + μ(dL i /dt) 2 in cell i for positive parameters γ and μ, imposing that total dissipation Φ = i Φ i is minimized subject to Φ = −dU/dt [15], where U = i U i . In the latter case (which is much more widely implemented in the literature), a drag force imposed on vertex k by cell i of magnitude 1 3 ηdr k /dt is added to f ik at each internal vertex for some η > 0, leading effectively to N v coupled ODEs for r k (t) of the form with f ik given by (4.1). Both instances lead to identical conclusions in terms of the structure of the force network and of the Airy stress function in (3.33), (3.34), for example, but differences in detail once the stress is expressed in terms of pressures and tensions. Crucially, however, (4.5) alone is insufficient to ensure moments are balanced across the monolayer. In summary, tracking variations of energy (and possibly dissipation rate) in terms of displacement of individual vertices (rather than in terms of strains, as in conventional elasticity), and imposing force balances alone, are insufficient in general to guarantee a torque balance. Extra constraints must be imposed on the evolution of the total energy U as it moves towards equilibrium. Conditions might be used, introducing Lagrange multipliers λ ik that ensure that each internal vertex lies at the orthocentre of the triangle formed by adjacent edge centroids (figure 4a). Conveniently, (4.6) involves cell edges and links between edge centroids that can be directly expressed in terms of vertex locations. Following the construction in appendix F, we can construct a dual network that identifies cell centres and links, up to a translation and scaling. The degree of freedom in scaling is accommodated by jumps in the Airy stress function across cell edges, but otherwise there is no impact on representations of stress.

Discussion
The planar vertex model describes cells as a network of polygons that tile a region of the plane.
We have shown that a natural dual network is one that connects cell centres (suitably defined) via the mid-points of cell edges, forming tristars around each vertex (figure 2). To represent force balances geometrically, further subdivision of these networks is required, into the links s ik between adjacent edge centroids and spokes q ik within each cell. The building blocks of the primal (cellular) and dual (tristar) networks are kites, defined by q ik ⊗ s ik . The antisymmetric part of this outer product gives the oriented area of the kite in cell i neighbouring vertex k; the symmetric part characterizes asymmetries in tristar shape via the fabric tensor F k , defined in (2.14). Stress in two dimensions (in continuum linear elasticity) can be written as the curl of a vector potential, which itself can be written as the curl of a scalar (the Airy stress function). In the present problem, we have shown that the Airy stress function ψ ik is defined on kites and curls are discrete: the vector force potential h j on edge j can be constructed as a curl of ψ ik taken over adjacent spokes q ik via (3.27), while cell stress σ i is a curl of h j taken around cell edges t j via (3.13); likewise, tristar stress σ k is a curl of h j around links T j between adjacent cell centres via (3.20). Jumps in ψ ik between neighbouring kites capture the projection of h j onto t j or T j : jumps across cell edges contribute to the isotropic stress, and jumps within cells across links contribute to shear stress. We find that t j · T j = 0 is a necessary condition for h j to be defined as a discrete curl of a potential having the appropriate jumps; equivalently, it is a necessary condition for a torque balance on cells and tristars. However, t j and T j need not intersect at edge or link centroids, and so networks differ in general from a classical (or radical) Voronoi construction [34]. We also identified the fundamental constraint (3.37) requiring that each cell vertex should be the orthocentre of the triangle formed by adjacent edge centroids (or, equivalently, of the triangle formed by the three adjacent vertices), from which we were able to develop a selfconsistent dual network (appendix F). Our strategy of using polygonal cell boundaries to define the primal network, and using physical constraints to identify an appropriate dual triangulation (specifically via an orthocentric construction), differs from many other studies in the discrete calculus literature in which a simplicial complex is taken to be primal and a priori barycentric or circumcentric constructions are used to build a dual network of polygons [28,35,36].
The force network and Airy stress function both provide mechanisms for visualizing stress. Stress can be interpreted as a map between the centroid network (with edges s ik and vertices c j ) and the force network (with edges −εf ik and vertices h j ). However, this map can be distorted, with the periphery of the force network, for example, shrinking to zero as the external load P ext tends to zero. The isotropic stress fields, P eff,i over cells or P eff,k over tristars (3.30), (3.31), are defined as discrete Lapalacians materials, one of a class of potentially significant mechanical heterogeneities [38]. However, further work is needed to identify the analogue in the present problem (if it exists) of the Beltrami-Michell equation (which leads to the Airy stress function satisfying a biharmonic equation in continuous planar elasticity), which would make P eff harmonic. A further useful visualization arises from the constraint that the orientation of intracellular stress σ k in the tristar that surrounds vertex k must share its principal axes with the fabric tensor F k , provided there is sufficient asymmetry for F k to be well defined. An analogous construction in granular materials has been connected to the orientation of force chains [39]. Remarkably, the stress-geometry condition does not depend directly on the choice of constitutive model.
In simulating the vertex model, it is common to either minimize a total mechanical energy U(r) directly by displacing the vertices r k , or to apply a drag η to each vertex, so that an equilibrium is reached by timestepping N v coupled ODEs for r k (t) of the form (4.5). In both approaches, the cell stress σ i can be evaluated and, happily, it is symmetric (4.2), ensuring local torque balance. However, this condition is not sufficient to ensure global torque balance, as consideration should also be given to the stress σ k around vertices. In other words, our study shows that cellular materials described by the vertex model should also be subject to a stress-geometry condition (3.39) equivalent to that described for granular materials [24,31]. Our study, therefore, suggests that it is necessary to constrain the optimization of U (for example, via candidate algorithm (4.6)), to ensure that appropriate geometric constraints are satisfied as the system approaches a final equilibrium state. A secondary construction identifying cell centres (appendix F), allows imposition of (3.26). We will address computational approaches with which to implement the torque balance conditions (3.26), (3.37) elsewhere.
This study is based on two fundamental assumptions: first, that the forces acting on each vertex can be partitioned into contributions from each cell, and that these constitute all the forces in the system; second, that there are no intra-or inter-cellular torques. From these assumptions, we deduced orthogonality of links and edges, and orthocentricity of vertices with respect to their neighbours. An alternative strategy was taken by [40,41], who partitioned forces at vertices (such as (4.1)) into contributions from each edge, deriving a triangulated dual network embedded in R 3 . The relationship between these approaches is discussed briefly in appendix G. An interesting further consequence of vertex orthocentricity is that all internal angles of polygonal cells should exceed π/2, implying that a T1 transition (a neighbour exchange) will arise as soon as the internal angle between adjacent cell edges becomes too acute. This is in contrast to the standard vertex model, when a threshold condition is often needed to trigger such a transition [42].
As figure 1d illustrates, orthogonality between links (connecting cell centroids) and edges (connecting vertices) is imperfect in real systems. There are obvious epistemic reasons: there are errors in the measurement and segmentation of cell boundaries; cell walls are not straight; and additional forces acting on some cells (due to division or motility), that are not easily partitioned at vertices, may be missing from the force balance (3.1). Furthermore, while cell centroids can be determined directly from images (using (3.42)), these will typically deviate from the cell centres that enable conditions such as (3.37) to be satisfied. Careful optimization strategies will, therefore, be needed to align self-consistent models that respect torque balance with data such as figure 1. It also remains to be seen to what extent the neglect of global torque balance has influenced predictions of previous computational realizations of the vertex model. The discrepancy may be subdominant to many of the other approximations implicit in modelling complex biological cells with simple physical models. For example, models that impose a Voronoi structure on the monolayer, solving only for the motion of cell centres, gain computational efficiency at the cost of some fidelity [19,43], at a level that has previously been judged acceptable for the purposes of the studies in question. Nevertheless, it is desirable to ensure physical balances are properly and fully respected, particularly as models grow in sophistication, and we argue that an orthocentric construction is more appropriate. At a more fundamental level, the appearance of a stress-geometry condition also raises intriguing questions about the role of microstructure in homogenized models of biological tissues. In summary, by identifying the underlying structure of the stress field implicit in the vertex model in terms of an Airy stress function and by identifying geometric constraints arising from torque balance, this study supports the development of more robust simulations, facilitates deeper understanding of the mesoscopic structures in disordered cellular monolayers, and provides a secure foundation for future upscaling approximations.
Ethics. All Xenopus work was performed using protocols approved by the UK Government Home Office and Acknowledgements. Conversations with Alex Nestor-Bergmann are gratefully acknowledged.

Appendix A. Incidence matrices
Treating the network of cells as primal, then a dual network is the triangulation between adjacent cell centres [28]. An orientation imposed on one network induces an orientation on the other. As figure 5 illustrates, consistency requires cells to share the same orientation, and triangles to share the opposite orientation. Consider cell edge t j , connecting vertices k and k , that is dual to the link T j between the centres of cells i and i . All possible values of the entries in A (indicating the orientation of t j with respect to vertices k and k and of triangles with respect to T j ) and B (indicating the orientation of T j with respect to vertices i and i and of cells with respect to t j ) are enumerated in the figure.
A T and B T have interpetations as boundary operators. Since all cells have edges that form closed cycles, A T B T = 0, and hence BA = 0 [28]. The rank-nullity theorem gives dim(ker(A T )) + dim(im(A T )) = N e and dim(ker(B T )) + dim(im(B T )) = N c . ker(A T ) identifies sets of edges with no boundary, i.e. edges that form closed cycles. The N c cell boundaries form a linearly independent basis for all such cycles, and hence dim(ker(A T )) = N c . For a localized monolayer, there is no combination of cells that has no boundary, and so dim(ker(B T )) = 0. It follows that B T , and hence B, is full rank (dim(im(B T )) = N c , i.e. the set of all cell boundaries is spanned by N c independent cycles), whereas A has rank dim(im(A T )) = dim(im(A)) = N e − N c (giving the size of the set of independent boundaries of edges).

Appendix B. Links between edge centroids
To establish ( Then B i j = 1 and B i j = −1, so that B ij A jk again takes the values −1 and 1 when summed over cells in the anticlockwise direction. Thus εε k i,j B ij A jk produces the signature (1, −1), (1, −1), (1, −1) when taken over the three cells surrounding vertex k, looping around the six edges of the tristar. Assuming instead that k = ε leads to reversals in the sign of εε k and of B ij , which cancel, leading to the same pattern, hence establishing (2.14).

( C 3 )
For the jumps in potential P and Q to depend independently on g · a or g · b then requires a · b = 0, in which case (C 3) gives the familiar orthogonal projection P = (g · b)b/b 2 and Q = −(g · a)a/a 2 . Finally, for a vector h, we seek the potential such that the jump P in the a direction is h · a and that the jump Q in the b direction is h · b. The necessary construction is a discrete curl around the