Couple stresses and discrete potentials in the vertex model of cellular monolayers

The vertex model is widely used to simulate the mechanical properties of confluent epithelia and other multicellular tissues. This inherently discrete framework allows a Cauchy stress to be attributed to each cell, and its symmetric component has been widely reported, at least for planar monolayers. Here, we consider the stress attributed to the neighbourhood of each tricellular junction, evaluating in particular its leading-order antisymmetric component and the associated couple stresses, which characterise the degree to which individual cells experience (and resist) in-plane bending deformations. We develop discrete potential theory for localised monolayers having disordered internal structure and use this to derive the analogues of Airy and Mindlin stress functions. These scalar potentials typically have broad-banded spectra, highlighting the contributions of small-scale defects and boundary layers to global stress patterns. An affine approximation attributes couple stresses to pressure differences between cells sharing a trijunction, but simulations indicate an additional role for non-affine deformations.


Introduction
The vertex model is a powerful tool for describing the mechanics of spatially heterogeneous multicellular tissues (Weliky and Oster 1990;Nagai and Honda 2001;Farhadifar et al. 2007;Staple et al. 2010;Fletcher et al. 2014;Alt et al. 2017).A confluent planar epithelium, for example, is represented as polygons tiling a plane.A mechanical strain energy is attributed to each cell that is a function of geometric invariants (such as the cell's area and perimeter) and the total energy of the monolayer is minimised, at a rate defined via a model of viscous dissipation, by varying vertex locations, potentially allowing for cell neighbour exchanges (so-called T1 transitions).A force balance at each vertex is used to evolve the system to equilibrium; elastic forces are defined by taking the first variation of each cell's mechanical energy with respect to vertex displacements.The changes of a cell's area and perimeter arising from small displacements of its vertices can thereby be used to define the mechanical (Cauchy) stress attributed to each cell.The model predicts a symmetric Cauchy stress tensor associated with each cell (Ishihara and Sugimura 2012;Yang et al. 2017) that aligns with cell shape (Nestor-Bergmann et al. 2018b) and allows viscoelastic moduli for bulk and shear deformations to be evaluated (Nestor-Bergmann et al. 2018a;Tong et al. 2022).Less attention has been paid to the Cauchy stress defined over the network that is topologically dual to cellular polygons, namely the triangulation connecting adjacent cell centres.The stress attributed to each triangle describes the mechanical environment in the neighbourhood of the tricellular junction lying within the triangle.This stress field is of interest given the role of tricellular junctions as potential sensors of cell shape and mechanical stress (Higashi and Miller 2017;Bosveld et al. 2018;Nestor-Bergmann et al. 2019;Angulo-Urarte et al. 2020;Yu and Zallen 2020).
From a multiscale modelling perspective, the vertex model is also of interest as a bridge between descriptions of discrete cells in a tissue and a continuum description of the tissue's mechanical properties (Murisic et al. 2015;Tlili et al. 2015;Ishihara et al. 2017).In two-dimensional (2D) continuum mechanics, it is often convenient to express the Cauchy stress in terms of a scalar potential, the Airy stress function (Howell et al. 2009).However in seeking to construct the discrete analogue of the Airy stress function, we found (Jensen et al. 2020) that the requirement for both forms of the Cauchy stress (that defined over cells, and that defined over tricellular junctions) to be symmetric places severe geometric constraints on cell shape, specifically that cell edges should be orthogonal to links between cell centres and that each vertex should lie at the orthocentre of the triangle formed by its immediate neighbours.These constraints are not met in typical simulations (nor, indeed, in real monolayers).This discrepancy can be explained in part by noting that while forces balance at vertices in the normal implementation of the vertex model, torque balance is not enforced.Here, we consider how the discrepancy can be accommodated by relaxing the requirement for all Cauchy stresses to be symmetric, by incorporating couple stresses within the constitutive framework.This approach is natural given the use of second-gradient or micropolar models to describe materials with microstructure (Trovalusci et al. 2015;Rizzi et al. 2019).
The Cauchy stress attributed to a cell (which hereafter we call the force stress, evaluated as the first spatial moment of the forces acting over a cell) can be partitioned into an isotropic component (defining an effective cell pressure) and a deviatoric component (describing the shear stress experienced by each cell) (Nestor-Bergmann et al. 2018b).Analogous quantities can be attributed to the triangles bounding tricellular junctions, having vertices at cell centres.Couple stress provides an additional measure of the stress arising from in-plane bending deformations that generate curvature in material elements.Using a standard version of the vertex model, we demonstrate here that while individual cells experience zero torque, a couple can be exerted around tricellular junctions.By considering second-order spatial gradients of a virtual tissue deformation, we show how the couple can be explained, in part, by considering the degree to which a cell is 'bent out of shape' via pressure differences creating moments acting across adjacent cell edges.However, our analysis reveals limits to potential analogies between the vertex model and continuum theories of couple-stress materials, likely associated with non-affine deformations occurring at small scales.
Our calculations are facilitated through the use of tools of discrete calculus (Grady and Polimeni 2010).In particular, incidence matrices capture topological relationships between cell vertices, edges and faces and enable the primal network of polygonal cells to be related directly to the dual triangulation connecting adjacent cell centres.Incidence matrices also provide the building blocks of the discrete differential operators needed to represent stresses using vector and scalar potentials.Unlike the three operators needed for normal continuum mechanics (grad, div and curl), we find that up to 16 different operators (4 grads, 4 divs and 8 curls) are required in two spatial dimensions, in the general instance when links between cell centres are not orthogonal to cell edges.These operators permit representations of spatially 2D vectors in terms of scalar potentials, via Helmholtz-Hodge decomposition.For 2D continuum elasticity, two potentials suffice for simply connected domains (the Airy stress function, plus an additional stress function defined by Mindlin for couplestress materials, Mindlin 1962;Hadjesfandiari and Dargush 2011).For discrete networks of cells, we find that up to eight potentials typically emerge, four defined over the network of cells, and four over the dual triangulation, although these reduce in number when edges and links are orthogonal.The potentials facilitate visualisation of stress patterns across a monolayer and their construction in terms of eigenmodes of scalar Laplacians, built using the geometry and topology of the cell network, reveals how stress fields are influenced both by the macroscopic shape of a localised monolayer and small-scale features such as topological defects in the organisation of individual cells.
We briefly review continuous couple-stress materials in 2D in Appendix A, following Hadjesfandiari and Dargush (2011).Key points to highlight are: (i) the Cauchy stress and couple stress vector can be written in terms of continuous Airy and Mindlin stress functions and Ψ in the form thus satisfying force balance div = , a torque balance relating the antisymmetric component of to curl and a compatibility condition (derived from a constitutive assumption) div = 0 ; (ii) the vector potential curl − grad Ψ is here expressed using a Helmholtz decomposition in terms of the two scalar potentials and Ψ ; and (iii) in the prin- ciple of virtual work, the strain 1 2 (∇ + ∇ ⊤ ) of a small- amplitude deformation ( ) is energy-conjugate to while 1 2 (∇ 2 − ∇(∇ ⋅ )) = −2 (a strain gradient), where is the so-called curvature, is energy-conjugate to .We seek discrete analogues of these relationships below (disregarding the compatibility condition div = 0 , as we make an alter- native constitutive assumption), starting by describing the nature of Helmholtz decomposition over a discrete cellular network and its dual triangulation.Key aspects of discrete calculus that we exploit are summarised in Sect.2, with details provided in Appendix B. In particular, we identify four Laplacians associated with the 16 operators, through which scalar potentials can (in principle) be derived to describe any vector field defined over the cell network.The eigenmodes of the Laplacians, which are shaped by the boundary of the monolayer and the organisation of individual cells within it, provide the building blocks for stress fields.In Sect.3, we show how stresses in a vertex model can be expressed in terms of a force potential (following Jensen et al. 2020) and determine the underlying scalar potentials.Using a standard constitutive model, we evaluate in Sect. 4 the couple stress.Here, the analogy between discrete and continuous descriptions is revealing but imperfect, as the couple stress determined as a rotational component of a vector force potential (as in (1)) turns out to differ from the vector that is energy conjugate to (under an affine approximation), for the particular constitutive model that we investigate.Results are illustrated by computations in Sect. 5 and discussed in Sect.6.

Discrete cellular calculus
Before addressing mechanical questions in Sect.3, it is necessary to develop relevant tools of calculus for quantities defined over a disordered cellular monolayer.We use topological (Sect.2.1) and geometric (Sect.2.2) objects to derive operators (Sect.2.3), in particular discrete Laplacians, enabling vectors to be represented in terms of scalar potentials (Sect.2.4) built from eigenmodes of Laplacians that are specific to the monolayer (Sect.2.5).

Cell topology
We consider an isolated cellular monolayer occupying a simply connected domain on the Euclidean plane, as illustrated in Fig. 1a.Adopting notation used in Jensen et al. (2020), vertices, edges and faces of the (primal) cell network are labelled by k, j and i, respectively (Fig. 1b), where i = 1, … , N c , j = 1, … , N e and k = 1, … , N v .Orientations are assigned to each object, and the topological relationships between edges and vertices, and faces and edges, are defined by the signed incidence matrices A jk and B ij , respectively.(Thus, A jk = 1 if edge j points into vertex k, A jk = −1 if edge j points out of vertex k, and A jk = 0 otherwise; B ij = 1 if edge j neighbours cell i and has congruent orientation, B ij = −1 if edge j neighbours cell i but has opposite orientation, and B ij = 0 otherwise.)and also specify topological relationships between cell centres (assumed here to be cell vertex centroids), links connecting cell centres and triangular faces of the dual network (Fig. 1c).Vertices within the interior of a monolayer are assumed to neighbour three cells: vertex/face neighbours are identified by the unsigned adjacency matrix = 1 2 , where and are unsigned incidence matrices (where Neighbour exchanges (leading to plastic deformations, with consequences for cell packing, Hashimoto et al. 2018) are incorporated in simulations but otherwise not considered in the present study, so that incidence matrices remain fixed.The topological identity = (which can be interpreted by say- ing that the boundary of any localised clump of cells is closed and therefore has no boundary) underpins the construction of discrete differential operators.We also identify centroids of each edge: each cell can then be partitioned into kites (labelled by ik, see Fig. 1b).In general, links between cell centroids in this network do not pass through edge centroids, except for Fig. 1 (a) A monolayer of N c = 112 cells, with N e = 345 edges and N v = 234 vertices, 'grown' using the vertex model via a sequence of random cell divisions.The primal network is defined by cell vertices (blue dots) and edges (black lines); the dual network is defined by cell centres (red dots) and links (white lines); white links also connect centres of border cells to peripheral edge centroids (green dots).A cluster of 7 cells, used in Fig. 4a, is highlighted with bolder colours.(b) A sketch defining geometric objects, their orientations and labels.Black lines denote cell edges, passing through vertices (blue dots) including k , k ′ and k ′′ ; blue lines denote links between cell centres (red dots), including i .Yellow: the kite of cell i at vertex k ′ with area K ik ′ , with two of its vertices at edge centroids (green dots) j ′ and j ′′ .Green: the trapezium with area 1 2 F j spanned by edge j and link j (orientations of other edges and links are not shown).Blue: the triangle surrounding vertex k with area E k .Also shown are the outward normals ij ′ to cell i at edge j ′ , and j ′ k ′′ , to triangle k ′′ at edge j ′ .Cell orientations i are opposite to triangle orientations k .(c) A diagram indicating how incidence matrices and map between vertices, edges, faces on the primal network, and cell centres, links and triangles on the dual network.Metric matrices −1 , e , l and −1 map between networks; e and −1 e map from edges to links and l and −1 l from links to edges.Loops indicate how the Laplacians V , T , C , F in (6) are constructed.
those in cells at the periphery of the monolayer.Thus, the faces of the dual network are internal triangles, or kites within cells at the periphery of the monolayer (either a single kite, or a pair of kites in adjacent cells).

Cell geometry
We introduce geometric information as follows.Points in the underlying Euclidean plane have position vector .Where necessary, p, q, r denote subscripts of vectors and tensors, identifying components with respect to a fixed basis in this plane, and a bold font is used to denote vectors in ℝ 2 .On the primal network of cells, we define vertices by k , edges by As indicated in Fig. 1c, acts as a difference operator in mapping vertex locations k to edges j .In contrast, ⊤ and ⊤ act as boundary operators.For example, the nonzero elements of are the edge centroids around the periphery of the monolayer, where c ≡ (1, 1, … , 1) is the N c -vector (a 2-chain, in the language of discrete calculus) denoting all cells in the monolayer.The number of edges of cell i is given by Z i = ∑ j B ij .We define the centre of cell i as Links on the dual network, triangulating cell centres, are defined by so that links either connect cell centres or connect centres of border cells (at the periphery of the monolayer) to peripheral edge centroids.Orientations of cell faces on the primal network i , and triangles (or peripheral kites) on the dual network k , are prescribed as ± , where the matrix (the 2D Levi-Civita tensor) represents a clockwise ∕2 rotation.i and k are taken to be independent of i and k, respectively, and of opposite sense.To ensure that and apply also to the dual network, orientations of edges j and links j are constrained such that where 1 2 F j is the area of the trapezium spanned by j and j (for interior edges, Fig. 1b) or the area of the triangle spanned by j and j (for peripheral edges).Consistent with typical simulations of the vertex model (Jensen et al. 2020), we allow edges and links to be non-orthogonal.
For cell i, the outward normal to edge j is ij = − i B ij j .Likewise the outward normal to the triangle connecting adjacent cell centres is jk = − k A jk j (Fig. 1b).The areas of cells and of interior triangles, A i and E k , satisfy where centroids of internal links are defined by ∑ j jk ⋅ j gives the area of each triangle at interior vertex k; E k is defined as the area of the adjacent kite (or kites) if k identifies a peripheral vertex (Fig. 1a).The total monolayer area A can then be written indicating how the monolayer can be partitioned into cells (labelled by i), trapezia spanned by edges and links (labelled by j) and triangles (plus peripheral kites, labelled by k).
To summarise, all topological information is encoded in and , while metric information is encoded in edge and link lengths t j , T j = | j | and in the areas A i , E k and F j .Using these we d e f i n e t h e m a t r i c e s ∕F N e with which we can define the square matrices eliminates 'orphan' links at the monolayer periphery that connect to only one cell centre.)The construction of these operators, which turn out to be the Laplacians for scalar fields defined over vertices, triangles, cell centres or cell faces, respectively (Appendix B), is illustrated in Fig. 1c and explained in more detail below.

Discrete operators
Variables (so-called co-chains) defined over cells, edges or vertices are written without the i, j, k subscript, so that {A} i = A i , { } j = j , {E} k = E k , etc.The scalar fields A, F and E are used below to define inner products with which discrete differential operators are constructed.In Appendix B, we describe how discrete analogues of grad, div and curl operators for scalar-valued variables defined on vertices or cell centres, and vector-valued variables defined on edges or links, can be defined.Figure 2 illustrates how the 16 operators act.Explicit expressions for the 8 so-called primary operators are given in Table 1.To summarise briefly, vectorvalued functions defined on edges or links sit in the isomorphic vector spaces E and L , respectively, which can be partitioned into subspaces E = E ∥ ⊕ E ⟂ (or L = L ∥ ⊕ L ⟂ ) of vectors parallel and perpendicular to edges (or links).Gradient operators grad v and grad c act on scalars defined at verti- ces and cell centres (in vector spaces V and C , respectively), creating vectors in E ∥ and L ∥ , respectively.Curl operators acting on scalars, curl v and CURL c , are rotated gradients that create vectors that are normal to edges and links, respectively (in E ⟂ and L ⟂ ).Divergence operators div v and div c measure fluxes of vectors normal to edges and links, mapping vectors from E ⟂ and L ⟂ to scalars defined over faces and triangles (in spaces F and T , respectively).Curl operators acting on vectors, CURL v and curl c , are similar, but meas- ure fluxes parallel to edges and links, mapping from E ∥ and L ∥ to F and T , respectively.Via the fundamental relation- ship = , these operators respect the exact relationships curl c • grad v = 0 , div c • curl v = 0 and so on, as summarised in Fig. 2. Superscripts c and v are used to denote objects associated with cells and vertices, respectively, and therefore primarily involve and , respectively.
When edges and links are non-orthogonal ( j ⋅ j ≠ 0 ), as we assume to be the case, a further eight so-called derived operators (adjoints under a suitable inner product, denoted with a tilde) must be considered (illustrated in Fig. 2); definitions are given in Table 1.Thereby we derive scalar Laplacians • curl v on vertices, given in matrix form by (6a).On the dual network, these have ana- triangles, given by (6b).The four scalar Laplacians reduce to two ( T = V , F = C ) in the special case of edge-link orthogonality, when F j = T j t j .

Helmholtz decomposition
A vector defined over edges or links can be represented in terms of potentials defined over each network, via a form of Helmholtz decomposition.(We do not provide a formal  that apply also to peripheral kites are given in (B12).
proof but call on analogous results developed for mimetic finite differences, da Veiga et al. 2014.)Assuming the cell monolayer is simply connected, any vector ∈ E with ele- ments j ( j = 1, … , N e ) has representation over the primal network of the form = ∥ + ⟂ , where for some v and and for some c and Ψ ).In ( 7), has been decomposed into its components parallel and perpendicular to each edge.
(and so on), we see that the potentials are determined by solving the Poisson problems The same vector can be represented over the dual network.Setting = ̌ ∈ L (where ̌ denotes representation on the dual network), the Helmholtz decomposition is given in terms of components parallel and perpendicular to links: thus ̌ = ̌ ∥ + ̌ ⟂ where for some ψ c , Ψc ∈ C and ψ v , Ψv ∈ T .The potentials are again determined from Poisson problems, namely In summary, (7-11) show how, given a vector field , we can determine the 8 corresponding scalar potentials that provide representations relative to the primal and dual networks (7, 10), by solving a sequence of Poisson problems (8, 11) using the four Laplacians given in (6).

Potential theory for monolayers
As demonstrated in Appendix B, the four Laplacians are self-adjoint under suitable inner products.The Poisson problems (8, 11) can then be evaluated directly using eigenmode decomposition, as demonstrated in (B18).We illustrate the basis in which solutions can be expressed by plotting in Fig. 3 the eigenmodes e C ) for the monolayer shown in Fig. 1a.V has identical structure (but slightly different entries) to T (likewise C and F ) and since links and edges are almost (but not exactly) orthogonal, the spectra and modes are qualitatively very similar.Each Laplacian has a zero eigenvalue with a uniform eigenmode c or v , where v ≡ (1, 1, … , 1) has N v elements.The eigenmodes of V and T are orthogonal under the inner product [⋅, ⋅] V and the eigenmodes of F and C are orthogonal under the inner product [⋅, ⋅] F (see (B6)). Figure 3 (top) shows how the lower-order eigenmodes are influenced strongly by the shape of the monolayer but show consistent patterns over cells and vertices; the eigenvalues of the first 20 modes of V and F differ by no more than 8.2%.
The highest-order modes exhibit strong localisation around defects in the monolayer (Fig. 3, bottom).

Discrete forces, stresses and potentials
We now deploy this methodology to formulate discrete potentials of stress fields in a generic vertex model, specifically a vector force potential (Sect.3.1) and its underlying scalar potentials (Sect.3.2).This reveals a couple stress, even when cells experience zero torque (Sect.3.3).

Vector force potential
A standard computational implementation of the vertex model yields forces ik (of cell i, acting on vertex k) that balance at each vertex and around each cell, so that (respectively) These balances can be interpreted geometrically by rotating each force by ∕2 (a form of Maxwell-Cremona force tiling, Bi et al. 2015), so that forces form closed triangles around each vertex (12a) and closed loops around each cell (12b), as illustrated in Fig. 4(b).The network of rotated forces is topologically equivalent to the network of links connecting adjacent edge centroids (Jensen et al. 2020), illustrated in Fig. 4(a).Just as edge centroids j provide a vector potential for these links, so the vertices of the network of rotated forces define a vector potential j for rotated forces, via such that j − j � = ∑ i ik , summing over a path connecting vertex j ′ to j (Jensen et al. 2020).The force stress c i of cell i (and the force stress v k defined over triangle k) can then be written as the first spatial moment of the forces acting on the cell (or triangle) (Nestor-Bergmann et al. 2018b), or equivalently in terms of the force potential (Jensen et al. 2020), as Here, differencing operators and in ( 13) have been transferred from to or , respectively, to express stresses in terms of edges and links.
The force-moment tensors in ( 14) are extensive and satisfy an important conservation principle (Bi et al. 2015): summing A i c i over adjacent cells yields a quantity defined entirely by forces or force potentials at the periphery of the cells (because of cancellation of internal forces for a system in equilibrium).In particular, for a monolayer of total area A that is subject to a uniform external pressure P ext (Nestor-Bergmann et al. 2018b).The periphery of the force network is defined by rotated peripheral forces: as these are assumed to act normal to peripheral edges, the rotated peripheral forces form a closed loop having exactly the same shape as the periphery of the monolayer, but scaled by P ext (Jensen et al. 2020).
The conditions that force stresses have zero divergence can be demonstrated by splitting stresses over cells and triangles into elements associated with edges and links, as explained in Appendix C. Dividing stresses into isotropic and antisymmetric components (Appendix C) also reveals that Thus, the projections of the vector force potential onto edges and links (in curls, (16c, d)) are associated with couples on cells and triangles, whereas the projections onto normals to cells and triangles (in divergences, (16a, b)) contribute to the isotropic stresses, which we express as the effective pressure P ef f .

Scalar stress potentials
We now pursue the discrete analogue of (1), expressing the force stress in terms of scalar potentials and identifying an associated couple-stress vector.The force potential can be expressed in terms of scalar potentials using ( 7) and ( 10), so that ( 16) becomes We construct vectors orthogonal to and ̌ (its dual rep- resentation) such that Then, by analogy with (A3), we can rewrite ( 14) as This representation is verified in Appendix C. Eq. ( 19) shows how the Cauchy stresses are defined in terms of the force potential , which is given in turn in terms of eight potentials (four per network) in ( 18).The relationship with (1a) becomes clear: where subscripts ⟂ and ∥ serve as reminders of the ori- entations of the vectors relative to edges (for c ) and links (for v ), respectively.We note also from (18) that suggesting how ( Ψ v , Ψ c ) and ( Ψv , Ψc ) can be interpreted as scalar potentials for the candidate couple-stress vector , in its representation over cells and over triangles, respectively (the analogue of (1b)).We note that curl c = − F Ψ c and CURL v ̌ = − T Ψv so that

Potentials with zero couple on cells
We now specialise to the case when individual cells experience no couple (the case relevant to the standard vertex model).To do so, we require F Ψ c = 0 in (17a), which implies Ψ c = a c for some constant a, which we can take to be zero without loss of generality.Thus Ψ c = 0 .Then curl c vanishes (so that (a) i = in ( 16)) but CURL v = T Ψv is likely to be nonzero (giving nonzero torque on triangles).We identify (from (21a)) the couple stress vector with −curl v Ψ v , which is normal to edges, satisfying div c = .
In general we expect ≡ ̌ to be described by nonzero Ψv and Ψc in (21b).
To ensure that the Poisson problems such as (8a, 11a) have solutions, we integrate the forcing terms over the monolayer.Making use of ( 15), ( 16) and (B13), we find Writing = −P ext + ̆ , where j is the intersection of edge j with link j , the identities Thus ̆ satisfies the solvability conditions necessary to invert L F for c and L T for ψ v .Thus, on the pri- mal network, we solve F  c = −div c ̆ = −div c ( + P ext ) to determine c .Similarly, on the dual network, we solve When edges and links are almost orthogonal, we find that C ≈ F and V ≈ T (in that the spectra and eigenmodes are almost identical) so that ψ c ≈  c , ψ v ≈  v and Ψv ≈ −Ψ v .We therefore illustrate just three of the seven nonzero potentials below.

Energetics
We now introduce a constitutive model by specifying the general form of the mechanical energy, introducing the principle of virtual work for pure strain (Sect.4.1) and for inplane bending deformations (Sect.4.2).

Virtual work in the vertex model
We introduce an energy per cell U i = U(A i , L i ) (assuming cells have homogeneous mechanical properties) and define a pressure and tension as P i ≡ U∕ A i and T i ≡ U∕ L i , respectively.The total energy of the monolayer is U = ∑ i U i + P ext A , where A = ∑ i A i and P ext is an exter- nal pressure applied to the periphery of the monolayer, as in (15).We therefore assume there is no moment traction at the monolayer periphery and only a normal force traction.U is a function of vertex locations, via the dependence of areas and perimeters on k .Suppose the monolayer is in a stationary equilibrium configuration (denoted with a prime) and consider virtual displacements k of its vertices.The expansion reveals the force U i ∕ k exerted at vertex k by cell i.Given that A i ∕ k = 1 2 ij A jk , it follows that ∑ i A i ∕ k vanishes at all internal tricellular junctions, so that P ext contributes to forces only along the monolayer's periphery, via virtual displacement of edge centroids: The principle of virtual work requires that, in equilibrium, U is unchanged by small independent displacements of each of the vertices.Equivalently, the sum of all forces at each vertex vanishes when the monolayer is at an equilibrium, i.e.
for all k, as in ( 12).The second variation in (24) captures weakly nonlinear effects and establishes the stability or otherwise of the equilibrium (Yan and Bi 2019), including any jamming/unjamming transition (Bi et al. 2015).We work below with the first variation, but consider how the forces organise into stresses acting over cells.
Consider variations that can be expressed as a smooth function of position under a deformation ( ) , i.e. we map vertices from . Suppose first that is linear in , so that the virtual (24) displacements are of the form k = 0 + ( + ) ⋅ ( � k − 0 ) , where = ⊤ and = − ⊤ are a small uniform strain and rotation, respectively, and 0 and 0 are constants.The principle of virtual work can be formulated by noting that variations in cell area and perimeter under this deformation are (Nestor-Bergmann et al. 2018a) where with no contribution from uniform translation ( 0 , 0 ) and rotation ( ).Then, the first variation in energy can be written This reveals the leading-order cell force-stress tensor c i that is energy-conjugate to and symmetric.The virtual work principle (i.e., at equilibrium, the energy does not vary for small but arbitrary strains ) recovers the bulk constraint (15) with c = c .The isotropic component of c gives the cell effective pressure as (Nestor-Bergmann et al. 2018b) Comparison to direct evaluation of c as a first moment of forces in (14) (Nestor-Bergmann et al. 2018b) shows the success of the affine approximation in this instance.

Strain gradients
We now take this argument a step further and consider displacements involving gradients of strain, allowing to be quadratic in .We continue to neglect effects that are quadratic in strain but account for first-and second-order deformation gradients ∇ and ≡ (∇ ⊗ ∇) , and reformulate the first variation in (24) in terms of We interpolate defor- mation gradients evaluated on vertices onto edge centroids and cell centres, using Taylor expansion to capture the leadingorder effect of spatial variations across any single cell.Accordingly, we use subscripts i, j and k to describe fields evaluated at cell centres, edge centroids and vertices, writing i ≡ ( � i ) , j ≡ ( � j ) and k = ( � k ) and so on.We retain second deriva- tives of but discard third and higher derivatives, assuming deformations vary over scales long compared to the size of individual cells.As shown in Appendix D, the changes in cell perimeter and area to this order become (25) In comparison with (25), we note additional terms.That involving � i − � i gives a correction pushing the evaluation of i towards the cell area centroid i .The third-order tensors i and i (see (D28, D29)) characterise the impact of strain gradients on cell perimeter and area, respectively.They are size-dependent, as is appropriate for objects that measure a gradient.∇ does not change perimeter to this order, but it alters cell area through the curvature i .Returning to (24), the energy maps from U 0 ≡ U � + P ext A � to using (28) and neglecting quantities that are quadratic in strains, where this order (by (D25)), so that P ext does not exert any moment on the periphery. (For sec- ond-gradient materials, terms that are energy-conjugate to gradients of and are identified as hyperstresses, Toupin 1964).However, direct comparison of (29) with the principle of virtual work for a continuum (A5) is not straightforward: deformations for which ∇ are eliminated but which retain are not possible via compatibility constraints, and the present affine (or Cauchy-Born) approximation does not account for possible local adjustments of vertex locations that ensure equilibration.Nevertheless, the comparison suggests c i as a candidate couple-stress vector, defined over cells.Given that ∑ j � ij = , c i vanishes for symmetric cells, for which t ′ j is uniform.We can repartition the contribution to the energy associated with c i to define its analogue on edges and links.As gradients in curvature across the monolayer will not play a role in what follows, we take to be uniform, and drop primes, to define the vector j attributed to edges via (29) where area is partitioned into trapezia of area 1 2 F j , associ- ated with edge/link j.P ext makes zero contribution to at all internal edges, but contributes along peripheral edges.has zero curl around cells (because it acts along normals to edges and sits in E ⟂ ) but has nonzero curl around triangles of the dual network: for internal vertices, Here, the pressure difference across edge j, ∑ i B ij P i , is mul- tiplied by t 2 j to give a moment, and the three contributions to the moment at the tricellular junction are summed at the vertex.C k vanishes if pressures are uniform ( ∑ i P i B ij = 0 ) or if the edges are of uniform size (because = ).We do not seek to impose any conditions on div c or div v .
To summarise, we now have two representations of couple stress.The vector in ( 21) is associated with the contribution of the vector force potential associated with curls.It is normal to cell edges (sitting in E ⟂ , so exerting zero couple traction on any cell because ij ⋅ j = 0 ), contributes to tor- ques around cell vertices via ( 22) and sums to zero over the monolayer via (23d).However, it is not energy-conjugate to the curvature , as shown in (D30).In contrast, the vector also sits in E ⟂ , is energy conjugate to via (29) but is not expressible as a curl of a scalar field defined on vertices.
has a direct physical interpretation as pressure differences acting over cell edges meeting at a vertex to generate a torque, but was derived under an affine approximation that needs evaluation.

Computations
We implemented the vertex model using the commonly used cell energy for which cell pressure and tension are linear in area and perimeter: P i = A i − 1 and T i = Γ(L i − L 0 ) .A vertex drag was implemented so that the system could relax to equilibrium under for some  > 0 .We chose Γ = 0.2 and L 0 = 0.75 , values for which the monolayer is in a jammed state (Farhadifar et al. 2007;Bi et al. 2015).An isolated monolayer under uniform external pressure P ext was established by starting with a small number of cells and allowing cell divisions to occur randomly for an interval; we examine configurations in which the monolayer has settled to an equilibrium state (Fig. 1a, Appendix E).
The forces ik acting at each vertex in the equilibrium state were rotated and assembled to form a force network, as illustrated in Fig. 4b for the cluster of cells shown in Fig. 4a.The three rotated forces around each internal vertex form a closed triangle, and the Z i forces around cell i form closed loops, confirming (12).For sufficiently large |P ext | , the force network may form a planar graph.However in general this is not the case, although the force network maintains the same topology as that of connections between adjacent edge centroids (Fig. 4a) (Jensen et al. 2020).The distorted force loops provide a striking illustration of spatially and temporally heterogeneous loading experienced by individual cells as the monolayer grows (Movie 1).
The vertices j of the rotated force network were then used to evaluate predictions of the model.We evaluated − 1 2 {div c } i and confirmed that it recovered P ef f ,i in (27) (Fig. 5a), while − 1 2 div v gives the corresponding effec- tive pressure partitioned over triangles (Fig. 5b).The fields show similar patterns over large scales.We validated the prediction that curl c = 0 (Fig. 5c) but found that CURL v typically is nonzero, being largest at the periphery (Fig. 5d) but heterogeneous over internal triangles (Fig. 5e).For comparison, we show (in Fig. 5f) the predicted couple C k (32).The patterns are distinct and differ in magnitude.The couple acting over trijunctions and at the monolayer periphery (16d), evaluated using simulations that incorporate nonaffine deformations, therefore differs from the couple (32) predicted via an affine approximation.
The potential c (Fig. 6a) underpins variations in the pressure field P ef f ,i , and is built from eigenmodes of F .Its spectrum shows contributions from a high proportion of modes (Fig. 6b).Because high-order eigenmodes are localised around defects (Fig. 3), the spectrum demonstrates the influence of these small-scale structures on the global stress field.The potential v (Fig. 6c) shows a similar distribution over the monolayer (except for a peripheral layer) and also has a broad spectrum (Fig. 6d).These potentials are smoother functions than P ef f , the latter being a second derivative of the former, but show that the pressure fields are not harmonic (unlike classical linear elasticity).Ψv is largest at the periphery (Fig. 6e), reflecting the structure of CURL v , and also has a broad spectrum (Fig. 6f).
Finally, a set of equilibrium monolayers were generated using a random cell division algorithm.c and v were evaluated to illustrate a range of possible patterns (Fig. 7).In each case, the two potentials resemble each other at the macroscopic scale, but also show heterogeneities at the smallest scales.Eigenvalue spectra are typically broad, although low-order modes can have prominent contributions.These examples show diverse patterns of residual stress within the equilibrium monolayers.

Discussion
Continuum mechanical models are widely used to describe biological tissues, and do so successfully over lengthscales that are large in comparison to a tissue's internal heterogeneities.However at scales comparable to individual cells, the inherent granularity of the tissue becomes evident.The vertex model (Weliky and Oster 1990;Nagai and Honda 2001;Farhadifar et al. 2007;Fletcher et al. 2014;Alt et al. 2017) is one of a class of discrete models of tissue mechanics that resolves stresses at the level of individual cells, exploiting the natural partitioning of space that they provide.This offers immediate advantages in modelling growth processes, by allowing cell division, expansion and rearrangement to be represented explicitly, and capturing growth-induced residual stresses.Likewise, explicit representation of individual cells facilitates the description of subcellular processes (such as the cell cycle, or cell signalling) and enables direct comparison with images.Continuum models rely on assumed strain energy functions, expressed in terms of strain invariants; in contrast, the vertex model  relies on a mechanical energy defined in terms of easily measured geometric invariants (such as the area or perimeter of cells in a planar monolayer).Despite these differences in the approach to constitutive modelling, a Cauchy stress can be defined in both instances.
In continua, it is commonly assumed that the Cauchy stress is symmetric, reflecting the absence of net torque on the smallest material elements; accordingly, stresses on material elements depend on local strains.In a discrete model, however, the smallest elements (e.g.individual cells) have finite size: stresses are specified primarily by local geometric measures (deviations in cell area and perimeter from target values serve as strains) but also by spatial gradients of bulk strain, 'measured' across the length of an individual element.Deformations that generate appropriate in-plane bending may thereby generate torques on tissue elements that are accommodated by so-called couple stresses.The present study is the first (to our knowledge) to address this feature in models of multicellular tissues, by evaluating the couple stress associated with the traditional vertex model.Our study of passive torques in epithelia is distinct from that of Yamamoto et al. (2020), who consider active cortical torque generation in a vertex model using a 'disk-shaft' mechanism.
The monolayers addressed here are deliberately simple: they are mechanically passive and confluent (Kim et al. 2021), and do not demonstrate fluctuations, motility or slippage of adjacent junctions (Nestor-Bergmann et al. 2022).The strain energy that we chose to investigate (33) passes a number of basic tests.Imposing forces at cell vertices is sufficient to ensure zero net force on each cell.This is demonstrated by closed loops in the plane of rotated forces (Fig. 4b); similar networks are used in granular flows (Ramola and Chakraborty 2017) and suspensions (Edens et al. 2021), and help visualise heterogeneous stress patterns.The model also ensures zero net torque on individual cells, and a cell force-stress tensor (26) that is symmetric (which we validated numerically in Fig. 5c).The stress of a cell can be constructed by summing contributions from individual vertices (or equivalently, from individual edges).These contributions can be repartitioned to evaluate the stress over the triangulation connecting cell centres.Here, in contrast, we find that the force-stress tensor is asymmetric (Fig. 5d, e), implying that a torque is exerted in the neighbourhood of each tricellular junction.A couple stress must be incorporated in order to accommodate the torque.
Bearing in mind the simplicity of the constitutive model (33), it is perhaps unsurprising that analogies between the vertex model and continuum models are imperfect.This can be anticipated, given the significance of non-affine deformations in fibre networks (Chandran and Barocas 2006), which can limit the accuracy of continuum approximations that assume affine deformations (Stracuzzi et al. 2022).One route to couple stress is to consider the rotational contribution to the vector force potential (21), that generates asymmetries in stress tensors defined over triangles spanning cell vertices.An alternative route considers the couple stress as a quantity that is energy-conjugate to in-plane bending deformations, represented by the curvature vector .Both routes indicate the existence of couple stress, as a vector  c (a, d, g, j, m, p) and v (b, e, h, k,  n, q), in 6 realisations of localised monolayers with corresponding eigenmode spectra (c, f, i, l, o, r) with amplitudes on a log 10 scale ( c spectra blue; v spectra orange).
1 3 defined over edges and links that has zero curl over cells, but the predictions differ in detail (Fig. 5e, f).The former route evaluates torques directly in terms of computed forces, accounting for non-affine deformations.The latter route rests on an affine assumption.The vector ( ) (31) is nevertheless of interest as it suggests a direct interpretation of torques at vertices arising from pressure differences between the three cells neighbouring a junction, acting over edges of different lengths, creating a net moment (32).It also suggests that couple stresses are intrinsically connected to spatial disorder, given that perfectly symmetric cells do not show the arearesponse to in-plane bending of asymmetric cells in (28b).However, it underpredicts the overall torque (Fig. 5e, f), likely because local equilibration at vertices leads to deformations not captured in an affine approximation (Chandran and Barocas 2006;Stracuzzi et al. 2022).
An array of confluent polygonal cells provides a natural unstructured mesh on which to perform computations.The machinery for pursuing such calculations is provided by discrete calculus, combining tools of algebraic topology (Desbrun et al. 2005;Grady and Polimeni 2010) with mimetic finite differences (da Veiga et al. 2014;Lipnikov et al. 2014).Incidence matrices and encode topological relationships between cell vertices, edges and faces, and the equivalent relationships over the dual triangulation.When combined with appropriate metric information, they can be used to construct discrete differential operators.By respecting the need to preserve exact sequences, a full set of operators can be identified, including positivesemi-definite discrete Laplacians defined over the primal (cell) and dual (triangular) networks (Table 1, Fig. 2), having eigenmodes (Fig. 3) with which potentials can be constructed.Helmholtz-Hodge decomposition (for a monolayer with no holes) enables a vector field defined over cell edges (namely, a force potential built from the forces acting on cell vertices) to be represented in terms of scalar potentials on each network.Thereby, we recover the discrete analogues of the Airy stress function of traditional 2D continuum elasticity, and the additional function introduced by Mindlin to describe couple stresses (Fig. 6).In general, the functions derived over networks of cells are distinct from those derived over the dual triangulation, although they share large-scale features (Figs.6a, c & 7).Broad eigenvalue spectra (Figs. 6, 7) implicate small-scale features near topological defects (Fig. 3) in overall stress patterns.
With this framework in place, we return to a question raised previously (Jensen et al. 2020), namely the consequence of neglecting torque balance in computational implementation of the vertex model.If couple stresses are assumed not to exist, so that all stresses are symmetric (over cells and over the dual triangulation), then cell edges and links between cell should, in principle, be orthogonal.Indeed, a stronger condition was identified (that vertices sit at the orthocentre of the triangle formed by their neighbours Jensen et al. 2020), which suppresses shearing deformations and is typically violated in real monolayers.Invoking couple stress relaxes the orthogonality (and orthocentricity) constraint, but reveals distributions of torques across the monolayer.These are largest at the monolayer periphery (boundary-layer features being characteristic of couple-stress materials, Toupin 1964) but distributed also across the interior of the monolayer (Fig. 5d, e).Couples are relatively weak in comparison with other stresses but they highlight tricellular junctions as sites where asymmetries in cell packing may be detected by mechanosensitive processes.
In the present analysis of force and couple stress in cellular monolayers, we have considered only systems at equilibrium, and not have not accounted for transient viscous effects or neighbour exchanges.However, this assessment of the vertex model demonstrates its utility in crossing scales from cell to tissue.Identifying the Laplacians of the cellular network opens the door for spectral methods to investigate global patterns of stress a tissue in a systematic way.Force chains within cell monolayers or associated with cells embedded in matrix (Mann et al. 2019) are an interesting target for investigation, as they may provide a mechanism for long-range mechanical signalling.More generally, this study also highlights a requirement to recognise that disordered multicellular tissues may need to be modelled at the macroscopic level as couple-stress materials, with boundary-layer effects (Fig. 5d) and torques at tricellular junctions (Fig. 5e) emerging as essential features.

Appendix A: Couple stress in 2D continua
In two dimensions, a continuous simply connected couplestress material in plane strain can be characterised by a force-stress tensor (having zero divergence, ensuring force balance) and couple-stress vector , with the antisymmetric component of expressed as a curl of (ensuring torque balance).There are three independent components of the symmetric component of force stress (s) ≡ 1 2 ( + ⊤ ) and two of , constrained by two (scalar) force balances and a torque balance.These constraints are satisfied by expressing and in terms of two potentials, the Airy stress function ( ) plus a second stress function Ψ( ) described by Mindlin (1962), such that (Hadjesfandiari and Dargush 2011) ensuring that p pq = 0 ( div = ) and p p = 0 ( div = 0 ).Here, is the 2D Levi-Civita tensor representing a clockwise ∕2 rotation; note that curl = ∇ for a scalar and curl = ∇ ⋅ ( ) for a vector , so that curl curl ≡ −∇ 2 , allowing (A1) to be written as (1).In Cartesians, (A1) becomes with y = x Ψ and − x = y Ψ .This formulation makes minimal constitutive assumptions beyond material continuity, except that the condition div = 0 is a com- patibility condition for an isotropic linearly elastic material rather than an equilibrium condition (Hadjesfandiari and Dargush 2011).The Airy stress function determines the isotropic component of the force stress via Tr( ) = ∇ 2 , while the Mindlin stress function Ψ deter- mines the antisymmetric force stress via the torque balance The force stress can be decomposed into isotropic, antisymmetric and symmetric-deviatoric parts as = 1 2 ∇ 2 − 1 2 ∇ 2 Ψ + (s) , where Tr( (s) ) = 0 . (s) has real eigenvalues ± , with ≥ 0 measuring shear, which depends on both and Ψ via s) , and In this sense, can be regarded as a vector potential for the force stress, and and Ψ can be regarded as scalar potentials of in a Helmholtz decomposition ( being the sum of a gradient of and a curl of Ψ).
The gradient of a smooth deformation ( ) can be decom- posed into + , where = ⊤ ≡ 1 2 (∇ + ∇ ⊤ ) represents strain and ≡ 1 2 (∇ − ∇ ⊤ ) =  is a rotation, where ≡ 1 2 ∇ ⋅ ( ) .Likewise, ≡ (∇ ⊗ ∇) can be decomposed as = ∇ + ∇ . is symmetric in its first two arguments, while contracting over them gives which defines a curvature vector = 1 2 ∇ (Hadjesfandiari and Dargush 2011).The corresponding principle of virtual work for a continuous couple-stress material occupying a volume V can then be written (Hadjesfandiari and Dargush  2011)   showing that the curvature vector is energy-conjugate to the couple-stress vector.Here, = ⋅ is a force traction at a surface with unit normal , and m = ⋅ is a couple traction.

Appendix B: Discrete operators
The discrete analogues of 2D differential operators are defined over primal and dual networks and act on variables defined on vertices, edges, and faces of each (Grady and Polimeni 2010;da Veiga et al. 2014).These are summarised in Fig. 2, which shows the primary operators on the cell network ( grad v , curl v , curl c and div c ) and on the dual network ( grad c , CURL c , CURL v and div v ).

B.1 Operators on the primary network
We define vector spaces V , E , F of fields defined on vertices, edges and faces, with associated inner products , where is the 2 × 2 identity.E is influenced by both edges and links via (3).Thus, for any , ∈ V , , ∈ E , f , g ∈ F .
The 'primary operators' (in the terminology of da Veiga et al. 2014) over cells are grad v ∶ V → E , curl v ∶ V → E , curl c ∶ E → F and div c ∶ E → F , and are defined as in Table 1.curl c and div c mimic the integrals arising in Stokes and divergence theorems.In matrix form, w h e r e V = , ).The topological relationship = ensures that curl c • grad v = 0 and div c • curl v = 0 .These exact sequences (de Rahm complexes) and definitions (B7) can be represented using the commutative diagrams and (B6a) where E is the direct sum of spaces of vectors parallel and perpendicular to edges Derived operators (denoted with tildes, following da Veiga et al. 2014) are defined as adjoints of the primary operators under the inner products (B6), satisfying for any ∈ V , ∈ E , f ∈ F .It follows from (B6, B8) that the derived operators have the matrix representations These relationships are summarised in the upper half of which van- ishes because ( ) ⊤ = .The derived operators are given in Table 1.
By specifying fields appropriately in (B8), we can use (B8) to write, for any ∈ V , ∈ E and f ∈ V , (B8a) (see Fig. 2), showing that these scalar Laplacians (given in matrix form in ( 6)) are positive-semi-definite.(Regularisation of C and F at the monolayer periphery, shown in (6b), is explained in Sect.B.4.) We can use (B8) to obtain the orthogonality relations which hold for any functions ∈ V and f ∈ F .These results rely on = , rather than geometric orthogonality, and they underpin the Helmholtz decomposition (7).
We also identify two Helmholtzians (closely related to Hodge Laplacians ⊤ + ⊤ ) acting on E, which sit in the spaces spanned by ̂ j ⊗ ̂ j and ( k ̂ j ) ⊗ ( k ̂ j ) , respectively.We assume that dim(ker(H 1 )) = 0 and dim(ker(H 2 )) = 0 , ensuring that there are no additional harmonic contributions to the decomposition (7).A necessary condition is that the domain has no holes (Lim 2020).

B.2 Operators on the dual network
Primary operators on the dual network are defined as follows: L is the direct sum of spaces of vectors parallel and per pendicular to links, L ∥ ⊕ L ⟂ .Here C = , L = diag( ⊤ j ) , ̃ L = diag(−( k j ) ⊤ ) , T = diag(E k ) .T , L and C are vector spaces of fields defined over triangles, links and cell centres.We note the isomorphisms C ≃ F .Derived operators are defined using the inner prod- ucts with metrics T = V , L = E , C = F , via sequences shown in the lower half of Fig. 2.This leads to operators given in Table 1, with actions illustrated in the lower half of Fig. 2, and the Helmholtz decomposition given in (10).
Additional steps are necessary to accommodate faces of the dual network at the monolayer periphery, which are not complete triangles.Writing ik = C ik ( k − i ) as the spoke connecting cell centre i to adjacent vertex k , we can alternatively write A geometric interpretation for an internal triangle is provided in Fig. 8a: the path around the boundary of a triangle can be reformulated into paths around the component kites (doubling back on all internal edges); since j is uniform over each link, the paths are equivalent to running up and down each spoke.Evaluation over spokes allows definition of adjoint operators as (B12a) Eq. ( B12) is advantageous as it holds for peripheral faces of the dual network, made from just one or two kites.This enables the discrete analogues of the divergence and Stokes theorems to be stated.Transforming sums over cells [triangles] to sums over edges [links], making use of derived operators defined in Table 1 and (B12a), it follows for a vector field that where p j are edges at the periphery of the monolayer.

B.3 Validation
To validate (B12), we define j to be the intersection of edge j and link j .For all internal intersections, ) i lie at the base of j and j , respectively.In general, j lies close to, but is distinct from, the edge centroid j .However we define j = p j on peripheral edges.(We must therefore distinguish m-kites with vertices j , shown in Fig. 8a, from c-kites with vertices j , shown in Fig. 1b).Writing ik = ∑ j B ij j A jk , then ik and ik span each m-kite within the monolayer.
∑ k C ik ik ⊗ ik gathers m-kites into cells, while ∑ i C ik ik ⊗ ik gathers m-kites into faces of the dual network (including the single and double kites at the periphery).Making use of ik ⊗ ik = K m ik + k , where K m ik is the area of the m-kite spanned by ik and ik and ik is a (symmetric) fabric tensor that measures asymmetries in kite shape (Jensen et al. 2020), it follows that ik ⋅ ik = 2K m ik and, because over the entire monolayer, providing a useful test for computations.In contrast, curl c and CURL v provide a map of the asymmetry of each face of the primal and dual networks (Fig. 8b,c).This is largest in the peripheral cells, but is distributed across the monolayer.The opposite orientations of cells and triangles account for the opposite signs of curls, for example, in the neighbourhood of the two triangular peripheral cells.

B.4 Inverting Laplacians
We can partition vertices into peripheral and interior vertices, edges into peripheral, border and interior edges, and cells into border and interior cells, so that (Jensen et al. 2020) Then, the unregularised Laplacian F = −1 e ⊤ becomes where The first contribution to F in (B15a) involves 'orphan' links between border cell centres and peripheral edge centroids via bp , whereas the remaining 'regularised' component of the Laplacian, r F , involves links between adjacent border cells via bb , links between border and interior cells via bi and links between interior cells via B ii .Since all such links lie between adjacent cells, r F satisfies r F c = , and we use this in building potentials.
V , and so ( V q − V p )[e V p , e V q ] V = 0 , demonstrating orthogonality (under the inner product) of eigenvectors having distinct eigenvalues.Writing = ∑ k c k e V k and projecting onto e V p , it follows that while (B17b) has forcing satisfying the solv- ability condition [ v , g] V = 0 .This enables its solution to be expressed in terms of the remaining eigenmodes, to obtain  =  + φ where Similar arguments follow for T under [⋅, ⋅] V , and r C and r F under [⋅, ⋅] F .The field, with uniform divergence but nonzero curl (Fig. 8b, c), allows evaluation of c and ψ v via (8a, 11a).

C.1 Microscopic force stress
The conservation principle (15) makes it possible to consider how A i c i or E k v k are built from their component parts.Accordingly, we define the microscopic force stress as Here, [ ] attributes each edge [link] component to a neighbouring cell [triangle] face but maintains it as a distinct entity from the other edge [link] contributions.It follows immediately that ij ⋅ ̃ c i = 0 for each edge j of cell i and jk ⋅ ̃ v k = 0 for each link j of triangle k, ensuring that div c ̃ c i = and div v ̃ v k = .We decompose j into components along j and i j , so that w h e r e h a t s d e n o t e u n i t v e c t o r s .T h e n , A i Tr( so that the deviatoric micro- scopic force stress becomes The final term involves ∑ j rather than ∪ j , ensuring that div c ̃ cD i = .̃ cD i has symmetric component and antisymmetric component where we have used ̂ j ⊗ ̂ ij − ̂ ij ⊗ ̂ j ≡ i B ij (consider its action on a vector  ̂ j +  ̂ ij ).We can therefore interpret j ⋅ j in (C21) as a torque exerted on each edge of the cell.Analogous expressions to (C20-C21) follow immediately for ̃ k , after projecting j onto j and k j .
The cell and triangle force stresses can be recovered from microstresses by replacing ∪ j with ∑ j in (C19), as in Jensen et al. (2020), to give (14).
We can also draw a distinction between c i in (30), the force stress integrated over cell i, and the corresponding microscopic cell stress which retains edge-to-edge variation rather than averaging over the perimeter.This stress has zero divergence, because evaluating div c ̃ c i includes ∑ j ′ ij summed around a closed loop, which vanishes, and , which also vanishes as � ij ⋅ � j = 0 along each edge.This ensures zero net force on each cell as in (12), and as illustrated by closed polygons in the -plane in Fig. 4b.

Appendix D: Non-uniform deformations
Here, we consider how area and perimeter change under affine deformations ( ) that vary with position, over lengthscales long compared to an individual cell.In the following, subscripts i, j, k attached to or its derivatives denote evaluation at i , j or k , respectively.Dropping third (and higher) spatial derivatives of in Taylor expansions, edges, edge lengths and normals map under the deformation to where (D23c) shows how the mapping is referred to an adjacent cell centre.Recall that ∇ = + captures stretching and rotation and = (∇ ⊗ ∇) = ∇ + ∇ captures their gradients.Here ij is the vector connecting cell centre i to an adjacent edge centroid j , i.e. ij = B ij ( j − i ) .The cell centre (the vertex centroid) is also the centroid relative to edge centroids (because Combining (D25) and (D26), links from cell centres to edge centroids map to Using (D25) and (D23b), cell perimeters change according to where Writing the final term as L � i i ⋅ ∶ ( ) i reveals the 3-tensor i characterising the impact of strain gradients on cell perimeter.Rotation gradients do not affect perimeter changes to this order.
Using (D23c) and (D27), and again dropping terms beyond third derivatives, the cell area becomes The terms involving i and i vanish because ∑ j � ij = .Consistent with (4), we note that where the integral over cell i is simplified by recognising that the unit normal is uniform along each edge.By extension, integrating along edge j of cell i, −t j ∕2 ( j + s ̂ j ) ⊗ ( j + s ̂ j ) ds = ij ⊗ j ⊗ j + 1 12 j ⊗ j , so that integrating around cell i where i ≡ A −1 i ∫ i dA is the area-centroid of the cell, which in general will be distinct from the vertex centroid i .It follows that (D27) (D28) This gives Then, we note that taking i = ± and noting that ij ⋅ i = ∓ 1 2 B ij t j,p ( p ) i (with as defined in Appendix A).Similarly, Hence, we can write w h e r e i h a s d i m e n s i o n s o f l e n g t h , w i t h Y i,pqr = ( � i,r − R � i,r ) pq + (24A � i ) −1 ∑ j B ij t � j,p (t � j,q n � ij,r − n � ij,q t � j,r ).
= A i �  pq ( i,r − R i,r ) +  pr ( i,q − R i,q )) � − 1 12 ∑ i B ij n ij,p t j,q t j,r .
t j,p t j,q n ij,r M i,pqr = −B ij t j,p t j,q i,rs t j,s ( p E qr + p W qr ) i = −B ij t j,p t j,q i,rs t j,s ( p E qr + p qr ) i = ∓ B ij t j,p t j,q rs t j,s ( p E qr + p qr ) i = ∓ B ij t j,p t j,q t j,s ( rs p E qr − qs p ) i = ∓ B ij t j,p (t j,q t j,s rs ( p E qr ) i − t 2 j ( p ) i ) = ∓ B ij t j,p t j,q t j,s rs ( p E qr ) i − 2t 2 t j,r n ij,q t j,p M i,pqr = ∓B ij t j,r t j,p t j,s qs ( p E qr (D29)

D.2 Energy and force potential
Expanding the energy using (13), without the boundary component, gives using (D23a).Using (14a), and taking i = ± this becomes For the chosen constitutive model, curl c vanishes and rota- tion has no impact on energy.The final term can be reformulated to give Gradients of strain and rotations will therefore generate hyperstresses conjugate to ∇ and , but that conjugate to is not equivalent to CURL v .

Fig. 2
Fig. 2 Four diagrams showing the action of operators defined on the primal network of cells (top) and dual network of triangles (bottom), involving vectors parallel (left) and perpendicular (right) to edges and links, respectively.In each diagram, primary operators run along lefthand vertical arrows.Derived operators, running along right-hand vertical arrows, are adjoint to primary operators under inner products (horizontal arrows) acting on elements of vector spaces V ≃ T

Fig. 4
Fig. 4 (a) The 7-cell cluster highlighted with bolder colours in Fig. 1a, with lines linking edge centroids j (green dots); (b) the rotated force network of the highlighted cells, with arrow colours matching the cell applying the corresponding force.The vertices of the force network are force potentials j .The topology and colours of the force network in (b) match those of the edge centroid network in (a).The data are from simulations (see Sect. 5) with parameter values Γ = 0.2 , L 0 = 0.75 , P ext = 0.2.

Fig. 5
Fig. 5 For the equilibrium state shown in Fig. 1a, we show (a) − 1 2 div c , giving the effective pressure in cells, (b) − 1 2 div v giving the effective pressure over internal triangles, (c) curl c , which vanishes over cells and (d) CURL v , giving a measure of the couple in the neighbourhood of each internal vertex, accounting for non-affine deformations.The field in (d) is largest at the monolayer periphery; (e) shows CURL v over internal triangles on a finer colour scale.(f) Couples C k computed using (32) under an affine approximation.The three colour bars (right) apply to (a, b), d) and (e, f) respectively.

Fig. 6
Fig. 6 Potentials (a) c , (c) v (analogues of the Airy stress function), and (e) Ψv of the Mindlin stress function), which provide representations of the force potential over the primal network of cells.(b, d, f) show the corresponding eigenmode spectra of c, e), plotting the amplitudes of coefficients in the spectral representation (B18).

Fig. 7
Fig. 7 Airy stress functions,c (a, d, g, j, m, p) and v (b, e, h, k, n, q), in 6 realisations of localised monolayers with corresponding eigenmode spectra(c, f, i, l, o, r) with amplitudes on a log 10 scale ( c spectra blue; v spectra orange).

Fig
both satisfied exactly: for example, − div

Fig. 8
Fig.8(a) A diagram illustrating how a curl around the edges of a triangle centred on vertex k (blue dot) can be modified to the path taken by thin purple lines, which reduce to integrals along spokes ik (orange), as expressed in (B12).Green lines ik connect intersec-

rV
is most easily evaluated by setting all elements of p e to zero when computing −1 e ⊤ , as in (6b).r C can be evaluated similarly, setting to zero all peripheral elements of −1 l j , cell centres map towhere i ≡ Z −1 i ∑ j B ij � ij ⊗ � ij (arising from averaging displacements around the edges of the cell, where

Table 1
A summary of the 16 discrete operators.The top eight are primary, the lower eight are derived.Those on the left-hand (right-hand) half of the table map to, or act on, vectors parallel (perpendicular) to edges or links, respectively.Alternative expressions for CURL v , div v , CURL v and grad v