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.

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-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-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-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-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-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). This primal cellular network generically shows 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.

Figure 1.

(a) An epithelium (animal cap) dissected from a Xenopus laevis embryo and adhered to a fibronectin-coated PDMS membrane, imaged by confocal microscopy; cell edges are identified with GFP-alpha-tubulin (green); cell nuclei with cherry-histone 2B (red). Some cell shapes are mapped out in magenta. (b) The segmented image, with each cell represented as a polygon bounded by vertices at its trijunctions. (c) Distributions of edge lengths (between trijunctions) and link lengths (between cell centroids). (d) The distribution of angles at intersections between links and edges (illustrated by the inset in (b)), peaking at π/2. (Online version in colour.)

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 confluent cells, represented as tightly packed polygons covering a simply connected region of the plane. We assume that an external isotropic stress Pext 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.

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) , , a set of oriented cell edges t (of length t), and a set of oriented cell faces (that we simply call cells) (of area A), . Here where ϵ = ±ε represents a clockwise rotation by ±π/2. (ε is the 2D Levi-Civita symbol satisfying ε = −ε, ; 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 , and 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 × N matrix has elements A that equal 1 (or −1) when edge j is oriented into (or out of) vertex k, and zero otherwise. The N × N matrix has elements B 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 and , identifying neighbours but not orientations. Further properties of and are given in appendix A. The N × N matrix has elements C that equal 1 if vertex k neighbours cell i and zero otherwise. Thus (summing over all vertices) defines the number of edges (and vertices) of cell i. We let R represent the centre of each cell, without specifying yet how it might be related to the cell’s vertex locations (where 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 peripheral and N − N interior vertices so that r = [r, r], edges (and relevant functions with subscript j) as N peripheral, N border and N − N − N interior edges so that t = [t, t, t], and cells (and functions with subscript i) as N border and N − N interior cells so that . 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 is an N × N matrix, etc., so that

Figure 2.

(a) An illustration of a localized monolayer. Blue lines show cell edges, meeting at vertices. This example has N = 7 cells (six border, one interior), N = 30 edges (18 peripheral, six border, six interior), N = 24 vertices (18 peripheral, six interior). Orientations of edges and faces are not indicated. Green dots are centroids c of each edge and red dots illustrate centres R of each cell. The solid orange lines connecting edge centroids form triangles around each internal vertex and polygons around each cell. Each cell is constructed from kites: three kites (shaded) surrounding an internal vertex together define a tristar. A force f due to cell i on vertex k is associated with each kite. (b) Solid purple arrows show rotated forces −εf. The force balances on vertices and cells imply that the rotated force vectors form a network that has the topology of the network containing edge centroids. The centroids c, therefore, map to vertices of the force network h (circular symbols). An imposed uniform pressure is represented by the peripheral forces, represented in part by supplementary links (dashed) that close triangles. (c) Kite ik, spanned by the vector q from the centre of cell i to vertex k and the vector s connecting the centroids of the edges adjacent to vertex k. The vectors bounding the kite are also indicated. (Online version in colour.)

Edges are defined by , with lengths . This defines the unit vectors . The perimeter of cell i is (summing over all edges). It follows (for later reference) that ∂L/∂r 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 defines the outward normal of cell i at edge j and c defines the centroid of edge j. Let , where x is a position vector in and integrate 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 . This can be understood by recognizing ϕ as the potential for position ; its discrete form ϕ = ϕ(r) jumps by along edge j, and the net change in potential vanishes around a closed loop because (appendix A), a device we will exploit later on. Also, as shown elsewhere (e.g. [19,21]), ∂A/∂r 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.

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 ϵ to all triangles (opposite to that in all cells), the orientations of links between cell centres are induced by the choice of and (appendix A), with link dual to edge t. For a localized monolayer, peripheral triangles and links are truncated; complete links are given by and , where and .

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.

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 denotes the set of paths over the edge-centroid network connecting j and j′, demonstrating how c is a discrete vector potential for s. As loops around any interior vertex k or any cell i are closed, it follows that More generally, the N × N matrix with elements can be combined with in (2.2) to give , because the boundary of the centroid network is closed, while diagonal elements of vanish (representing closed loops around interior vertices); all diagonal elements of 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(q − r + R) = 0 (figure 2c).

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 is the sum of the two edge centroids bounding kite ik, so that c3 and c4 run anticlockwise around the kite (figure 2c). The area of the kite is given by (see appendix C). Following [31], we can write where c1 and c2 run anticlockwise along cell edges. The area of tristar k is therefore . Summing kites over the tristar, the internal edge contributions (involving c1 and c2) cancel leaving only boundary contributions, giving The fabric tensor 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 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. Stating (3.1) more generally to account for boundary forcing, we require the diagonal entries of to vanish (balancing forces at each vertex, including the periphery), and the diagonal entries of to vanish (an internal force balance on each cell), where the matrices and share the structure of in (2.2) and . Now summing over cells and vertices, respectively, and the external force (imposed pressure around the monolayer periphery) has matrix blocks where and . Thus zero diagonal entries of give (triangular) force balances at interior vertices (including contributions from peripheral cells where appropriate). Zero diagonal entries of give (polygonal) force balances over interior and peripheral cells, and zero diagonals of give the force balance on peripheral vertices.

For a monolayer satisfying (3.1), the first moment of the force defines the stress σ associated with cell i via We call the isotropic component of the stress in each cell the effective pressure, . The stress σ of the monolayer as a whole may then be written as where , restricting the final sum to peripheral cells and peripheral vertices because interior forces balance via (3.1). Imposing the boundary condition (3.3) gives the conservation law [22] showing that the total stress must be isotropic, internal shear stresses must cancel and therefore that [21] Equation (3.6) also ensures zero net torque on the monolayer due to Pext.

We now consider how the force balances (3.1) can be represented geometrically, with a view to identifying the (intercellular) stress σ defined over tristars.

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 and , 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 Pext > 0, anticlockwise if Pext < 0 and collapses onto the origin if Pext = 0.

The centroids c form a discrete potential for the edges s via (2.2), (2.9). Similarly, we can identify the vertices of the force network h (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 becomes redundant when A and B both appear in the sum, and using . Taking a transpose gives σ 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

Likewise , so that , giving the stress in rotated basis. It follows that Comparison with (from (2.6)) suggests that stress can be associated with a mapping from c to h. Equation (3.14) also shows that the isotropic component of cell stress is a discrete divergence,

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 and −εf. An explicit construction for such a map was provided in [31]. The mapping between vertices c, c, c, running anticlockwise around a triangle surrounding vertex k (figure 3), to h, h, h, also ordered anticlockwise, is where the triangle area a satisfies 2a = (εs)s = (εs)s = (εs)s and s = c − c, −εf = h − h, etc. The action of the map is demonstrated via

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 connects cell centre R to vertex r, intersecting the link s between neighbouring edge centroids c and c. Kite ik is spanned by q and s; the three kites neighbouring vertex k, that together form a tristar, are shaded. The Airy stress function ψ (see below) is defined on kites. Jumps in ψ between neighbouring kites in the same cell, sharing vertex c, are defined by the projection of the force potential h on the cell edge t. Jumps in ψ between neighbouring kites in the same tristar, sharing vertex c, are defined by the projection of h on the link T between cell centres (not shown) that intersects edge t. (Online version in colour.)

With the map defined, we can write . Then from (3.4), This shows how the cell’s stress is built from the kite shape tensor q ⊗ s in (2.13), weighted by contributions from the cell’s vertices. Accordingly, the stress over a tristar built from the three kites surrounding vertex k is where the fabric tensor is given by (2.14).

The stress over the tristar at vertex k can also be written in terms of spokes q and vertex forces f. Replacing the latter with the force potential h and reordering, we find that multiplying each h is the difference between neighbouring spokes, i.e. the straight link T between cell centres. In general, this does not pass through c, so contributing to non-zero . Explicitly, using (3.9), where and ϵ is the orientation of the triangle surrounding vertex k. (As explained in appendix A, we assume that all cells have the same orientiation ϵ, and that ϵ = −ϵ uniformly; with fixed j and k, B + B = 0 for all options in figure 5 below, so that the r terms cancel in (3.20).) The outward normals to triangle k are N = −ϵ AT. It follows that, analogously to (3.15), the isotropic component of tristar stress is

Figure 5.

Four different cases having consistent orientations of the link t (blue) between cell vertices k and k′, and the link T (red) between cell centres i and i′. The corresponding orientations of cells i and i′, and triangles k and k′ are also shown. The tabular inset shows corresponding values of incidence matrices. Cell and triangle orientations are given in terms of the Levi-Civita symbol ε, representing clockwise π/2 rotation. (Online version in colour.)

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 ψ (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 , which automatically satisfies (because —see appendix A). Likewise, zero torque on tristar k requires, from (3.20), , which we satisfy with potential ψ satisfying for both pairs of kites bounding edge j (i.e. with ). Pursuing the analogy with planar elasticity, we seek to define h as a discrete curl of ψ, here evaluated over the spokes q of the four kites surrounding c (e.g. the path q − q + q − q in figure 3) in order to recover the stress in terms of ψ.

To illustrate the definition of ψ, 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 ψ over neighbouring elements, and , so that and

We now express h in terms of the ψ 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 b, 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, t and T in (3.22)), it is necessary for a · b = 0. We therefore require i.e. each link between adjacent cell centres must intersect the corresponding cell edge orthogonally. Equation (3.26) is therefore necessary for both σ and σ to be symmetric (equivalently, for each to be expressed in terms of ψ). For the jumps in h · t and h · T to align appropriately with t and T, 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 ϵt defines a normal to edge j relative to cell i, we see that noting that, for all four cases in figure 5 below, (ϵt) · T = tT. 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 kites and 2Z independent values of ψ. Likewise, noting that (ϵt) · T = tT, 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 ψ, 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 Now and 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 ψ field between neighbouring kites intersecting the tristar. Likewise, using the identities and , we find In cell i, the final sum in (3.34) allocates a scalar (, the pairwise difference in ψ values taken in the same sense as the orientation of cell i) to each edge and then sums the outer products of the unit tangent and the inward unit normal, weighted by v and taken anticlockwise, such that .

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 as (when ϵ = −ε), where α is the orientation of edge j with respect to a fixed axis, then the final sum in (3.34) is (for ϵ = ±ε) This is symmetric because . Thus which we will make use of shortly.

Tristar stress and the fabric tensor

We now reconcile the two expressions for σ in (3.19) and (3.20). First, consider the condition Link T crosses edge t, bounded by two vertices, each surrounded by a triangle of vectors s of the edge centroid network having area a; two such triangles are illustrated in figure 3. Depending on the chosen orientation of cells and edges, (3.37) implies that link T is parallel (or antiparallel) to the furthest edge of each triangle (such as s in figure 3), with the magnitude of T relative to the edge s given by the ratio E/2a. In other words, (3.37) implies that each vertex bounding t lies at the orthocentre of the triangle of edge-centroid-links surrounding each vertex.

Direct substitution of (3.37) into (3.20), giving σ as an outer product of links with the force potential, recovers with defined in (3.16), given as an outer product of edge-centroid-links with the force potential. Thus (3.37) is equivalent to the condition Symmetry of σ is ensured by the existence of the Airy stress function and the orthogonality condition (3.26); (3.37) extends this symmetry to . Furthermore, (3.19) then implies that , while (3.38) gives , yielding the stress-geometry condition [24,31,32] The role of stress as a mapping between networks is also evident via f = σεs, 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(σ) 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(σ) ≠ 0, we can write σ in terms of its principal axes and eigenvalues as σ = σe ⊗ e + σe ⊗ e, where e · e = 0. Likewise, as long as , we may express it in terms of its principal axes as where f · f = 0. Then (3.39) implies 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.

Relating cell centres and cell vertices

We have not yet specified how cell centres R might be related to cell vertices r, so that conditions (3.26) and (3.37) may be satisfied. The orthogonality condition (3.26) applies to all border and internal edges: the links are , and the edges , . Then (3.26) requires, for N − N border and internal edges The R correspond to 2N scalar quantities. Given a set of vertices, for 2N > N − N 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 < N − N (anticipating that N ∼ 3N and for N ≫ 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 , where

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 of cell i define its edges , edge centroids and the internal links between them . Equation (3.37) requires that the edge radiating outwards from cell i at vertex k must be orthogonal to s; the triangle of links around vertex k is then fully defined since r 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/λ, say, under the constraint (3.37), to cover the plane (figure 4a), with R the common vertex of all triangles covering cell i. Noting that 2a is the product of the base and the height of the edge-centroid triangle at vertex k, link j should have length T satisfying (from (3.37)) where h is the altitude of the triangle at vertex k (figure 4a). Given that the area U of the triangle spanned by maps to a via , with altitudes related by h = λ H, (3.43) implies E = Th = 2λ U = 2a/λ, implying a result that can be verified by (2.13).

Figure 4.

(a) Geometric constructions demanded by torque balance. Red triangles connect adjacent edge centroids. Each cell vertex is at its orthocentre, i.e. the (black) cell edge passing through a vertex of a red triangle is orthogonal to the opposite side of the red triangle. The links between adjacent cell centres (blue) are orthogonal to the cell edges they intersect. Thus each red triangle surrounding a vertex is similar to the blue triangle surrounding the same vertex (opposite edges are parallel), differing by a rotation of π and uniform scaling. Dashed lines are orthogonal to cell edges and dotted lines are orthogonal to triangle edges. The angles marked with arcs are equal. h indicates the altitude of the edge-centroid triangle at vertex k. (b) An illustration of three cells (black) satisfying orthocentric constraints and a corresponding network of cell centres (blue), constructed using an algorithm described in appendix F. (Online version in colour.)

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. (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 within cell i. Since all such internal angles sum to (Z − 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 and perimeter L, is assumed to have mechanical energy , so that cells have identical mechanical properties but distinct shapes and sizes. 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 and a tension 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 = δU/δr, the first variation of the energy of cell i with respect to a small displacement of vertex r. This determines the restoring force at r due to cell i as We can evaluate the cell stress using (2.6) and (4.1), noting that C becomes redundant, as where , 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 , A changes by and L changes by [15]. The cell structure tensor commutes with the cell shape tensor , 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.

Writing (4.2) as , comparison with (3.34) shows that the intracellular differences in ψ are proportional to . More specifically, writing with respect to some fixed Cartesian axes, recalling that . 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 times a nonlinear but dimensionless measure of cell shape.

We note that, since is defined in terms of areas and perimeters, it clearly satisfies material frame indifference. Under a change in reference frame, in which , where is a rotation and Z a translation, it is straightforward to show that (as do other vectors contributing to f in (4.1)), while the identity ensures that other tensors, such as stress, transform under , 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, in (4.1) is replaced with with , and with , corresponding to a dissipation rate in cell i for positive parameters γ and μ, imposing that total dissipation is minimized subject to Φ = −dU/dt [15], where . 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 is added to f at each internal vertex for some η > 0, leading effectively to N coupled ODEs for r(t) of the form with f 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 (3.26), (3.37) together suggest that a constrained energy minimization of the form might be used, introducing Lagrange multipliers λ 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 between adjacent edge centroids and spokes q within each cell. The building blocks of the primal (cellular) and dual (tristar) networks are kites, defined by q ⊗ s. 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 , 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 ψ is defined on kites and curls are discrete: the vector force potential h on edge j can be constructed as a curl of ψ taken over adjacent spokes q via (3.27), while cell stress σ is a curl of h taken around cell edges t via (3.13); likewise, tristar stress σ is a curl of h around links T between adjacent cell centres via (3.20). Jumps in ψ between neighbouring kites capture the projection of h onto t or T: jumps across cell edges contribute to the isotropic stress, and jumps within cells across links contribute to shear stress. We find that t · T = 0 is a necessary condition for h 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 and T 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 self-consistent 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 and vertices c) and the force network (with edges −εf and vertices h). However, this map can be distorted, with the periphery of the force network, for example, shrinking to zero as the external load Pext tends to zero. The isotropic stress fields, Peff, over cells or Peff, over tristars (3.30), (3.31), are defined as discrete Lapalacians of the Airy stress function, where , , serve the role of Hodge star operators [28]. Ramola & Chakraborty [37] used the spectral properties of a graph Laplacian as a tool to understand force localization in granular materials. Likewise, the geometrically weighted Laplacians (5.1), defined on the vertices of the primal and dual networks, are promising candidates for determining the structure of mesoscopic patterns of stress in cellular 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 Peff harmonic. A further useful visualization arises from the constraint that the orientation of intracellular stress σ in the tristar that surrounds vertex k must share its principal axes with the fabric tensor , provided there is sufficient asymmetry for 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, or to apply a drag η to each vertex, so that an equilibrium is reached by timestepping N coupled ODEs for r(t) of the form (4.5). In both approaches, the cell stress σ 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 σ 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 . 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.