Geometry of epithelial cells provides a robust method for image based inference of stress within tissues

Cellular mechanics plays an important role in epithelial morphogenesis, a process wherein cells reshape and rearrange to produce tissue-scale deformations. However, the study of tissue-scale mechanics is impaired by the difficulty of direct measurement of stress in-vivo. Alternative, image-based inference schemes aim to estimate stress from snapshots of cellular geometry but are challenged by sensitivity to fluctuations and measurement noise as well as the dependence on boundary conditions. Here we overcome these difficulties by introducing a new variational approach the Geometrical Variation Method (GVM) which exploits the fundamental duality between stress and cellular geometry that exists in the state of mechanical equilibrium of discrete mechanical networks that approximate cellular tissues. In the Geometrical Variation Method, the two dimensional apical geometry of an epithelial tissue is approximated by a 2D tiling with Circular Arc Polygons (CAP) in which the arcs represent intercellular interfaces defined by the balance of local line tension and pressure differentials between adjacent cells. We take advantage of local constraints that mechanical equilibrium imposes on CAP geometry to define a variational procedure that extracts the best fitting equilibrium configuration from images of epithelial monolayers. The GVM-based stress inference algorithm has been validated by the comparison of the predicted cellular and mesoscopic scale stress and measured myosin II patterns in the epithelial tissue during early Drosophila embryogenesis. GVM prediction of mesoscopic stress tensor correlates at the 80% level with the measured myosin distribution and reveals that most of the myosin II activity is involved in a static internal force balance within the epithelial layer. In addition to insight into cell mechanics, this study provides a practical method for non-destructive estimation of stress in live epithelial tissues.

Cellular mechanics plays an important role in epithelial morphogenesis, a process wherein cells reshape and rearrange to produce tissue-scale deformations. However, the study of tissue-scale mechanics is impaired by the difficulty of direct measurement of stress in-vivo. Alternative, image-based inference schemes aim to estimate stress from snapshots of cellular geometry but are challenged by sensitivity to fluctuations and measurement noise as well as the dependence on boundary conditions. Here we overcome these difficulties by introducing a new variational approach -the Geometrical Variation Method (GVM) -which exploits the fundamental duality between stress and cellular geometry that exists in the state of mechanical equilibrium of discrete mechanical networks that approximate cellular tissues. In the Geometrical Variation Method, the two dimensional apical geometry of an epithelial tissue is approximated by a 2D tiling with Circular Arc Polygons (CAP) in which the arcs represent intercellular interfaces defined by the balance of local line tension and pressure differentials between adjacent cells. We take advantage of local constraints that mechanical equilibrium imposes on CAP geometry to define a variational procedure that extracts the best fitting equilibrium configuration from images of epithelial monolayers. The GVM-based stress inference algorithm has been validated by the comparison of the predicted cellular and mesoscopic scale stress and measured myosin II patterns in the epithelial tissue during early Drosophila embryogenesis. GVM prediction of mesoscopic stress tensor correlates at the 80% level with the measured myosin distribution and reveals that most of the myosin II activity is involved in a static internal force balance within the epithelial layer. In addition to insight into cell mechanics, this study provides a practical method for non-destructive estimation of stress in live epithelial tissues.
Cell and tissue mechanics is an important factor that both affects and regulates animal and plant development and thus is a subject of active study in developmental biology and biophysics, reviewed extensively in [1][2][3][4][5][6][7]. Here we focus on the mechanics of animal epithelial cells that compose tissues in the form of two-dimensional monolayers with tight junctions between adjacent cells and adhesion (of the basal cellular surface) to the substrate extracellular matrix (ECM) [8].
In the absence of a rigid substrate, mechanical properties of such monolayers are dominated by the tissue-wide mechanical network formed by cytoskeletal cortices coupled by intercellular adherens junctions [9][10][11][12]. Cytoskeletal cortices are localized to the lateral sides of cells, just below the apical surface, and are made of actin fibers cross-linked by myosin II motors that actively generate tension within the cortex. The shape of cells within the tissue is determined by the balance of local actomyosin cytoskeletal contractility and the intracellular osmotic pressure [2,9], which acts to oppose the decrease in total cellular volume [13,14]. For the purpose of tissue-scale mechanics, the full three-dimensional force balance that shapes individual cells can be approximated by an effective two-dimensional model of the apical cytoskeletal network. In this simplified 2D view, the contractility of the junctional actomyosin "belts" balances against an effective two-dimensional pressure that prevents the collapse of the apical area under cortical tension: this type of an approximation underlies the widely used "vertex model" approach to epithelial cell mechanics [1,[15][16][17].
Measuring mechanical properties of cells and tissues invivo presents a considerable experimental challenge. AFM [18] and optical tweezer contact microscopy [19] have been used to probe the local rheology of individual cellular inter-faces at great resolution but do not provide a direct readout of internal stress. The most common method for detecting stress in vivo is UV laser ablation, in which focused light 'cuts' the cytoskeletal bundle abutting a cell-cell interface and the resultant retraction velocity is used as a proxy for the local cortical tension [20]. This method is convenient, as it does not require any special preparation of the sample, but is destructive and hence does not allow measurement of the global stress distribution across the tissue. Other methods use genetically encoded FRET tension sensors engineered into load carrying proteins [21] or employ measurements of deformation with implanted beads [22,23]. These methods require specially prepared samples and are technically challenging, both in implementation and in quantitative interpretation.
The difficulty of direct experimental measurement of mechanical stress in developing tissues has stimulated alternative approaches that seek to capitalize on the availability of live imaging data [24][25][26]. For example [25,26] introduced a method for inferring cellular stress from observed 2D cell geometry based on the assumption that the tissue is instantaneously in a mechanical equilibrium, described by a model parameterized directly by intercellular pressure and tension. In the simplest versions of the method [25,26] 2D geometry was parametrized by a polygonal tiling generally used in Vertex Models [1,15]. The validity of this approach, of course, rests upon the accuracy of the assumptions and approximations, which varies between tissues and conditions and has to be evaluated in each case. The major challenge to the approach stems from the sensitivity to both noise in the measurement of cellular geometry and the prior distribution imposed on the boundary conditions at the edge of the observed tissue domain. These difficulties have necessitated additional assumptions, introduced to over-constrain the inference problem [25,26]. A particularly direct generalization aims to extract additional information from the image data, most obviously provided by the observable curvature of cell interfaces [25,27]: an example of such an approach is the CellFIT toolkit [27]. Sensitivity to noise and image quality however remains a major issue. Below, we develop a new approach, improving the image-based "mechanical inference" to the extent that makes it broadly applicable in the practice of experimental data analysis.
In this paper, we explore general constraints imposed by mechanical equilibrium on the geometry of 2D cellular arrays that balance arbitrary interfacial tension against differential cellular pressure. Our analysis will identify a certain mathematical duality between the geometry of cells and a triangulation formed by the equilibrium values of interfacial tensions. This duality provides a set of highly nontrivial constraints that can be used to both effectively de-noise image analysis and stabilize the problem of mechanical inference, which is achieved through a variational formulation of the inference problem. Using synthetic data as a comparative benchmark, we show that our algorithm correctly infers mechanics under arbitrary pressure differentials and moderate measurement noise, performing significantly better than existing methods. To illustrate its practical utility, we apply the algorithm to live imaging data from the early stages of Drosophila embryonic development and demonstrate its ability to accurately predict -based on cell geometry alone -the spatial distribution and anisotropy of myosin II, the molecular motor known as a generator of mechanical stress in the developing embryo [2,[28][29][30]. Synthetic and real data tests suggest practical utility for the new mechanical inference algorithm as a tool for quantifying stress distribution in live tissue in the absence of direct measurement of local forces.
Equilibrium properties of 2D cell arrays with non-uniform pressure In order to model the mechanical state of the 2D epithelial tissue, we generalize the standard Vertex Model [15] which represents the epithelium by a planar polygonal tiling parameterized solely by the positions of vertices, the location where three or more cells meet, hereafter denoted r i . Here, we shall approximate the geometry of cells in the epithelial layer by a Circular Arc Polygonal (CAP) tiling, replacing the straight polygonal edges with circular arcs that correspond to tensed interfaces balancing pressure differentials between the adjacent cells as described by the Young-Laplace (Y-L) law [31]. The equilibrium geometry of a CAP tiling is fully specified by the set of effective interfacial tensions (T αβ , where α, β labels the cells partitioned by the given edge) and the set of effective hydrostatic pressures p α representing the contribution of bulk stress from the cell to the 2D apical force balance. In mechanical equilibrium the radius of the circular arc forming αβ interface, R αβ , is given by the Y-L Law, R αβ = T αβ /[p α −p β ]. A graphical example of the resultant CAP network is shown in Fig. 1.
The dimensionality of the space of all CAP tilings is given by the number of internal degrees of freedom: two positional degrees of freedom for each vertex plus the radius of curva-B Height above the plane 0 1 Figure 1. Circular Arc Polygonal (CAP) tiling and its dual tension triangulation. A) Circular arc polygons provide an approximate representation of equilibrium geometry of a cell array. Curvature of interfaces, is controlled by pressure differences between adjacent cells and is related to interfacial tension by the Young-Laplace law [31]. B) Force balance at vertex i requires the three tension vectors tangent to each edge to sum up to zero, which defines a local triangle. As adjacent vertices in A) share edges, the tension triangles of each vertex form a triangulated surface -the dual representation of force balancewith each triangular face corresponding to each vertex. In the absence of pressure differentials the dual surface would be flat: height variation arises from differential pressure.
ture for each edge. The total count is 2v + e = 7c where c, e, v respectively denote the number of cells, edges and vertices within the tiling, which satisfy v = 2c and e = 3c.
(This count is exact on a torus and is approximate up to the boundary corrections, for a large planar array). Conversely, the mechanical state of the network is parameterized by e + c = 4c parameters corresponding to independent interfacial tensions and cell pressures. Hence, the dimension of the space of all equilibrium CAP tilings is 4c and thus, mechanical equilibrium implies, on average, three constraints per cell. Our analysis below will aim to i) define a relation between the tiling geometry and tensions and pressures and ii) uncover the geometric constraints associated with mechanical balance. Mechanical equilibrium of a CAP network is reached when tensions balance at every vertex and interfacial curvature obeys the Young-Laplace Law. Tension acts tangentially along network edges, Fig. 1 depicts an example of force balance for vertex i due to tensions of interfaces that connect adjacent vertices j, k, l. We note that an edge can be labeled in two equivalent ways: either with the labels of cells the edge partitions or the vertex labels that the edge connects, e.g. edge ij in Fig. 1a separates cells α, β, so we shall define T ij = T αβ and use both schemes interchangeably. Hence force balance requires where ρ ij and R ij is defined as the centroid and signed curvature of the circular arc associated to edge connecting vertices i, j. The *-superscript here is shorthand notation for the counter-clockwise rotation by π/2, e.g. r * ≡ẑ ∧ r.
Each of the three terms is recognized as a tension vector, e.g. T ij = T ij (r i − ρ ij ) * /R ij , acting on vertex i. Hence, force balance can be re-interpreted geometrically as a triangle formed by tension vectors acting on vertex i, with adjacent cells α, β, γ, see Fig. 1 being associated with the vertices of the tension triangle. Crucially, as local force balance triangles associated with adjacent vertices share an edge, stitching together force balance conditions on all network vertices defines a tension triangulation 'dual' to the CAP tiling. Dual triangulation vertices correspond to cells and triangular faces correspond to the vertices of the original CAP tiling, as shown in Fig. 1. The tension triangulation dual to a CAP tiling is not planar. From the definition of the dual triangulation, the triangle angles are simply related to the angles at CAP vertices: φ iα = π − φ iα . If the edges of the tiling polygons were straight, one would find that iφ iα = 2π, corresponding to a planar triangulation [17]. As edges are curved, tension vectors acting at either end of an interface are no longer parallel: the force vector acting on vertex i due to interface αβ is rotated by an angle ∆ϕ αβ = αβ /R αβ from the equivalent force acting on vertex j, where αβ is the physical arc length of the circular edge. This results in the non-planarity of the tension-triangulation manifested by the "deficit" angle for each cell, ∆ϕ α , defined as the sum of curvature contributions from all edges that compose the cell: The deficit angle ∆ϕ α associated with cell α is the discrete Gauss curvature for the corresponding vertex in the dual triangulation.
We now proceed to define equilibrium constraints on CAP geometry and explicit relations between tensions and pressures and geometric observables. By applying the Sine Law to a triangle one can relate the ratio of tensions in adjacent edges to the corresponding angles, e.g. for the triangle i: T αβ /T αγ = sin(ϕ i,γ )/ sin(ϕ i,β ). Multiplying such ratios for a set of triangles that share a vertex of the dual triangulation uncovers a non-trivial constraint on CAP tiling angles (where the product is taken over the set V α of vertices i that belong to cell α, while β and γ label other cells adjacent to i in clockwise order). Eq. (3) defines c non-trivial constraints on the angles of an equilibrium CAP network which we recognize as the generalized form of the geometric compatibility condition introduced in [17]. Provided that the geometric constraints on the CAP tiling given by Eq.
(3) are satisfied, the dual tension triangulation specifies all T αβ up to an overall scale -the key property making inference possible. Next, given the set of tensions T αβ from the tensiontriangulation, pressures can be computed on the basis of the Young-Laplace Law by solving a discrete Poisson equation on the dual triangulation: Eq. (4) represents c equations that define c pressure unknowns up to the homogeneous solution which has to be fixed by the boundary conditions on p α . Importantly, the Y-L law implies that T αβ /R αβ + T βγ /R βγ + T γα /R γα = 0 must be satisfied at each vertex (e.g. vertex i shared by cells α, β, γ). Using the Sine Law on the dual triangulation, this constraint is recast in the purely geometric form which defines v = 2c (one per each vertex) geometric constraints that account for the difference in the number of pressure variables (c) and the number of R αβ variables (e = 3c). Together, Eqs. (3,5) impose 3c constraints on a equilibrium CAP network so that the dimensionality of the latter is indeed given by 4c, which is equal to the total number of tension and pressure variables. In principle, one can infer the underlying balanced mechanical state -the set of equilibrium interfacial tensions {T αβ } and pressures {p α } -by building the tensiontriangulation: utilizing local network angles at each vertex to build each dual triangular face, and subsequently solving Eq. (4) to obtain cell pressures. This approach is contingent upon the network satisfying compatibility conditions Eqs. (3,5). In practice it suffers from two major problems: (i) real cellular networks undergoing morphogenesis are expected to be close to, but not exactly in, mechanical equilibrium, (ii) the measurement of network geometry from imaging data will always be noisy and imperfect. As a result, an algorithm attempting to stitch together a global tension triangulation would rapidly accumulate errors that would dramatically impair the resultant inference. Furthermore, the solution of Eq. (4) requires the knowledge of boundary conditions which are generally not known.  To achieve a robust implementation of the mechanical inference we shall recast it as a variational problem on an unconstrained set of variables. The idea is to fit the observed network geometry to the closest equilibrium CAP tiling in the least-squares sense by utilizing a set of variables that both natively satisfies the constraints of equilibrium and fully characterizes the mechanical state of the tissue. We begin by parameterizing the CAP tiling by circular arc centroids ρ αβ and the associated radii of curvature R αβ and minimize the variational pseudo-energy function r αβ (n) denotes the position of n th pixel on edge α, β obtained directly from the segmentation and N αβ denotes the number of pixels segmented per edge. Eq. (6) has a simple geometric interpretation -it penalizes the Euclidean distance between the estimated and measured circular arc for each pixel of a segmented edge. The 9c fitting variables in Eq. (6) correspond to an arbitrary CAP tiling -we want to restrict to the subset of equilibrium CAP tilings. To define the reduced set of degrees of freedom we revisit Eq. (1) and note that the prefactor for each tension vector can be reinterpreted using the Y-L Law, leading to a relation between the arc centroids of the three edges connected at each vertex Remarkably, at equilibrium, the centroids of all three interfacial circular arcs meeting at a vertex are constrained to be collinear, with distance along the line controlled by relative pressure differences, depicted graphically in Fig. 2a. This is a practically useful observation that provides a direct method to assess the compatibility of any cellular network with mechanical equilibrium solely on the basis of image analysis of cell morphology. Additionally, in order for all three edges, corresponding to collinear centroids ρ αβ , to intersect at a single point defining a CAP vertex, their respective radii must obey which imposes v constraints on e curvature variables, R αβ , reducing the set to a e − v = c independent degrees of freedom (as expected from the fact that edge curvatures are generated by intracellular pressures which make c variables). Given any solution, we can generate another geometrically compatible cell array by transformingR 2 , where z α provide the explicit degrees of freedom for the c dimensional manifold of solutions of Eq. (8) that share the same set of edge centroids.
We next observe that Eq. (7) admits the general solution where we have introduced a vector variable q α for each cell. Substituting this form into Eq. (8) we obtain an explicit general expression for R αβ : Together Eqs. (9,10) provide an explicit parameterization of equilibrium CAP tiling by an arbitrary choice of 4c independent variables {q α , z α , p α }. These variables explicitly satisfy all geometric constraints and provide an explicit local expression for the tensions via the Y-L law. We note also that {q α , z α } variables have a nice geometric interpretation as the set of centroid points in 3D which generate equilibrium CAP tilings via generalized Voronoi construction defined in Appendix A, where we also derive an explicit relation between {q α , z α } and CAP vertices given by , which can be solved to define r i from the reduced variables {q α , z α , p α }. Note that pressures (and hence tensions) are defined only up to an overall scale factor, which leaves geometry invariant.
The Geometrical Variation Method for inferring internal stress distribution.
Eqns. (9, 10) parameterize ρ αβ and R αβ in terms of a reduced set of degrees of freedom, {q α , p α , z 2 α }, that guarantee the geometric constraints of an equilibrium CAP network. Given a grayscale image of cellular boundaries in an epithelial monolayer, comprising a large number of cells  with interfaces resolved at a pixel level, we can find the best approximation of the 2D network formed by cell interfaces. This is achieved by minimization of Eq. (6) with respect to this reduced set of degrees of freedom produces the best fit of image data by an equilibrium CAP tiling. Minimization of Eq. (6) is a non-linear optimization problem subject to linear inequalities on z 2 α 's that ensure positivity of the argument of the square root in Eq. (10). We solve the problem computationally using MATLAB's implementation of the interior-point algorithm. The simple choice of starting the iteration with z α = 0 and the estimate of {p α , q α } from the observed set of network angles (as explained in detail in the Appendix B) was found to produce reliable convergence.  . Validation of GVM for stress inference using synthetic data. A) Comparison of actual and inferred tensions using the "matrix inverse" inference method introduced in [25] made over-determined using measured edge curvature data, as in [27]. 10% noise was added to the synthetic input data (see Appendix C for details). ρ denotes the Pearson's correlation coefficient. B) Same as (A) but using the GVM on the same synthetic data. Note the significant increase of the correlation coefficient.
To evaluate the performance of GVM on mechanical inference, we tested it against synthetically generated data. As shown in Fig. 3, we find that inference of the dual triangulation is quite robust to both measurement noise and large inter-cellular pressure differentials (see the Appendix C for details). Key to the improved performance is the over-constrained nature of the present formulation of the mechanical inference problem which combines the estimation of geometric parameters with stress inference into a single variational analysis of image data. The redundancy does not only stabilize inference in the presence of noise, it also allows us to infer boundary stresses using just information of 2D bulk cellular morphology. The method is local and thus immediately applicable to inference of inter-tissue forces, as will be discussed below.
The inferred set of interfacial tensions and cellular pressures allows us to construct a stress tensor for an epithelial layer in mechanical equilibrium. Specifically, inside the bulk of cell α, the stress tensor is isotropic and constant σ ab = p α δ ab (where a, b index spatial dimensions) and σ ab = T αβr a αβr b αβ δ(|x − ρ αβ | − R αβ ) on edge αβ. One can average this quantity over cell area A α : where {β} α denotes all cells connected to cell α and integration is taken along the cell boundary arcs. In practice it is useful to coarse-grain the tensor by averaging over neighboring cells cells to obtain an approximation to the continuum stress tensor. In the following section, we verify the utility of the proposed mechanical inference by comparing the inferred stress tensor to known biological correlates of stress.
Local and global stress inference and its in vivo correlates.
We now apply GVM to imaging data from living epithelial monolayers. Ideally we would test the GVM inference against direct measurements of stress. Presently, the most reliable readout of local stress in live tissue is provided by observed levels of (fluorescent labelled) junctional myosin, which has been previously demonstrated to correlate with local interfacial tension measured by laser ablation [2]. We shall carry out the comparison between inferred stress and observed myosin level first on the scale of cells, then on the scale of the whole tissue, using the data on Drosophila embryonic development.
During the initial stages of Drosophila embryonic development, the ellipsoidal monolayer of epithelial cells forming the embryo undergoes a series of non-trivial mechanical transformations. Immediately following the formation of the ventral furrow (VF) -the first step of gastrulation -Drosophila embryo undergoes germ-band extension (GBE): a major morphogenetic movement involving a convergent extension of the lateral ectoderm, which approximately doubles its length along the embryo's anterior-posterior (AP) axis. This process has been demonstrated to be driven by the activity of the junctional pool of myosin II, which exhibits a non-uniform and anisotropic distribution on the surface of the embryo, in particular, forming contractile supracellular cables that run along the dorso-ventral (DV) axis of the embryo [28,30,32]. Laser ablation assays have demonstrated that these myosin cables, associated with DV oriented edges, exhibit significantly higher cortical tension than AP oriented cell junctions [32,33]. The quantitative relation between myosin and mechanical stress was further elaborated in our earlier study of morphogentic flow [28] which demonstrated that a symmetric 2D tensor m ab describing coarse grained distribution of myosin is a useful proxy for the stress tensor.
Applying GVM to the embryonic epithelium images [28] we found that cell-array geometry observed over the first 60 minutes of convergent extension, is quite well approximated by equilibrium cell network: E ≈ 1, which means our best fit equilibrium CAP geometry differs from the image segmentation by on average one pixel per edge! Fig. 4 shows results of the analysis of the lateral ectoderm during GBE. Qualitatively, inferred stress exhibits anisotropic stress cables that run along the DV axis in agreement with previous studies [32]. For a quantitative comparison, we compute the correlation coefficient of tension inferred on individual cellular interfaces with the myosin line density measured on same interfaces: a histogram of the calculated results for each time point is shown in Fig.4D. The mean correlation coefficient, ρ ∼ .4, is a two-fold improvement over the earlier matrix inverse method [25].
The observed correlation indicates that the inference method is picking up underlying mechanical effects. There are however numerous sources of noise that weaken correlation. To name two: i) while the analysis was carried out on a single snapshot, cell geometry is fluctuating on the timescale of seconds; ii) linear density of (fluorescentlabelled) myosin is not an exact measurement of line tension. Most importantly, the assumption that cells are in a mechanical equilibrium is at best only an approximation: in the case of GBE, there is a mean morphogenetic flow of cells indicating the presence of unbalanced local forces within the tissue [28]. Below we shall demonstrate that accounting for this systematic deviation from mechanical equilibrium effect improves correlation between inferred stress and measured myosin distribution at the mesoscopic scale.
Edge by edge comparison is the most exacting test as it is sensitive to local fluctuations. Myosin distribution however also exhibits non-trivial variation over the surface of the embryo (see Fig. 4E) and it is informative to compare it with the coarse-grained inferred stress tensor Eq. 11. In constructing the latter, it is helpful that GVM inference can be carried out locally on partially overlapping image frames and "stitched" (see Appendix D for details) into a continuous coarse-grained stress field for the whole surface of the embryo (as illustrated in Fig. 4F). In Fig. 4D we compare the inferred coarse-grained stress tensor with the measured coarse-grained myosin tensor [28]. Both myosin and inferred stress are enriched in the lateral ectoderm and are anisotropic along the DV axis, with the quantitative agreement at the level of about 60% (which corresponds to average ∆ ∼ .4 in Fig. 4F). In Fig. 4D we compare the inferred coarse-grained stress tensor with the measured coarse-grained myosin tensor [28].
To more exactly quantify the difference between the measured myosin tensor m ab (as defined in [28], ) and the inferred stress tensor σ ab we defined a normalized root-meansquare (r.m.s.) deviation: where m, σ abbreviate m ab , σ ab respectively and ... r stand in for averaging over a coarse-graining region at location r and λ is the unknown overall scaling factor relating myosin and stress which we chose so as to minimize the global average of ∆. Hence 1 − ∆(r) is the measure of local agreement (in both magnitude and anisotropy) between of σ and m within a coarse-graining region r.
In general, the inferred stress tensor is consistent with the measured myosin tensor -both exhibit strong anisotropy localized to the lateral ectoderm, with principal axis along the DV axis. The most substantial disagreement between the inferred stress and measured myosin tensors, as shown in Fig. 4F, is localized at the center of the dorsal side of the embryo; myosin is spatially inhomogeneous along the DV axis (high at the lateral sides) in contrast to the inferred stress tensor, observed to be constant along the DV axis. A plausible explanation for this discrepancy is that we are directly comparing the total myosin tensor, which is thought to drive early morphogenetic flow [28] and thus is unbalanced, to our inferred stress tensor which explictly assumes mechanical equilibrium.
Our previous study [28] related observed meso-scale myosin distribution m ab and observed morphogenetic flow, by focusing on the divergence of the myosin tensor ∇ a m ab which corresponds to the "unbalanced" internal stress within the tissue that generates cellular flow. (Note, following [28] we assume σ ab ∼ m ab .) However, only a fraction of myosin contributes the unbalanced stress, the rest generating internal stress obeying force balance: m ab = m ab U +m ab B with the "balanced" fraction of myosin m ab B being divergence-less, ∇ a m ab B = 0. It is the latter component of myosin that is expected to correlate with the predictions of GVM-based inference.
To decompose the measured myosin tensor [28] into "balanced" and "unbalanced" components we note that any 2D symmetric tensor can be written as (where ε ac = −ε ca is the antisymmetric unit tensor). The divergence of the last term is zero and it can be identified as m ab B , the balanced component of the myosin tensor. Taking the divergence of Eq. (13) yields a partial differential equation for the vector field u b : which was solved using the same method as in [28], see Methods. (Strictly speaking Eq. (14) defined u a only up to a harmonic gradients u a → u a + ∇ a ψ + ε ac ∇ c ω with ∇ 2 ψ = ∇ 2 ω = 0. Since the only solution to the latter equations on a closed surface of genus zero is a constant, our solution for u a is unique.) Eqns. (13,14) provide an explicit determination of the balanced component of myosin tensor: Fig. 5 displays the distribution of m ab B on the surface of the embryo and compares it with the total myosin distribution and the inferred stress. As shown in Fig. 5D balanced myosin dominates accounting for more than 80% of the total, but the unbalanced component increases with time, especially upon the onset of the GBE (10min post CF). We also find that removing the unbalanced component from the myosin distribution being compared to the inferred stress, substantially increases the agreement between the two (see Fig. 5D). Specifically, during GBE the alignment of our inferred stress tensor relative to the total myosin tensor decreases to ∼ 60%, consistent with the higher fraction of unbalanced myosin and thus morphogenetic flow.
Another advantage to the GVM inference is it can be applied to cell arrays with arbitrarily large pressure differences. Let us now provide an example of how our variational stress inference can be applied to study interesting questions concerning mechanical control of biological phenomena in the context of cell division orientation.
It is known that the spatio-temporal patterning and orientation of cell divisions plays an important role in morphogenesis. Mitotic domains of synchronously dividing cells partition the Drosophila embryo in a highly regular manner that directly shape eventual larval segments [35]. Additionally, the patterning of mitotic spindle orientation has been suggested to contribute to elongation of the posterior region of the lateral ectoderm [36]. While the upstream signal that instructs orientation of cell cleavage plane is unknown, studies suggest that mechanical tension within the tissue contributes to spindle alignment [37,38] in the dividing cell.
To test the hypothesis that cell cleavage axis tends to align with local tension within cells, we analyzed 70 tracked divisions in mitotic domains 6 and 11 (as defined in [35]) during the late phase of GBE (20-35min post CF). Fig.  6AB, provide an example of using the GVM method to infer tensions in individual interfaces of a cellular network. The cell cleavage axis was compared directly to the orientation of the tension axis determined from the inferred stress tensor, Eq. (11). We found that cell cleavage indeed correlates strongly with inferred tension and that the principal axis of stress is a much better predictor of spindle orientation at the time of division than the commonly used "long axis" defined directly by cell elongation. Furthermore, the GVM-based inference was more accurate (based upon 95% confidence intervals using t-test) than the earlier "matrix inverse" approach of [25]. The improved accuracy is due to the GVM's ability to capture large pressure differentials between cells in a morphologically heterogeneous cell arrays, as exemplified by Fig. 5A.

Discussion
The mechanical inference method described here was based on the model which assumed that (i) the 2D epithelial cell array is instantaneously in an approximate mechanical equilibrium, (ii) cell mechanics can be approximated by the balance of cytoskeletal tension localized at cell interfaces (and varying from one edge to another) and the effective areal pressure (preventing collapse of the apical surface of cells). Together, these assumptions place a nontrivial constraint on cell geometry that is readily testable on the basis of imaging data. The existence of local geometric constraints facilitated the formulation of a local mechanical inference scheme that combined estimation of stress with the simultaneous determination of the best-fit cell geometry from imaging data. The method was observed to be a significant improvement over similar methods, both at the scale of individual cells and the mesoscopic scale.  [34]. Time corresponds to roughly 26 minutes after cephalic furrow formation. C) The correlation between myosin and inferred tensions comparing the matrix inverse method defined in [25] and the new GVM-based inference algorithm during the first 40 minutes of GBE. D) An example of stress inference, combining two overlapping image frames used for the analysis, shown directly on the embryo and on a cylindrical projection with the ventral line cut and mapped onto the top and bottom edges of the image with the dorsal side along the midline. Inset shows color-coded edge tension inferred by the GVM algorithm. The inferred stress displays stress cables as expected for the lateral region of the embryo at the time considered. E) The mesoscopic anisotropy of myosin, 26min post CF formation, shown on the embryo and its cylindrical projection. Anisotropy is largest in the lateral region. The principal axis of myosin in this region points along the DV axis. F) The spatial distribution of the normalized r.m.s. difference (∆(r)) between the GVM inferred stress tensor and the measured total myosin tensor (at 26min post CF). Both tensor fields are represented as ellipses for direct comparison. The large discrepancy (∆ ∼ .8) found in the vicinity of the anterior and posterior poles (mapped to the left and right edges of the cylindrical projection) is due to the poor imaging of the poles. The large discrepancy in the center of the dorsal region is real and can be explained by the difference between the total and "balanced" myosin distributions as explained in the text.
Interestingly, equilibrium CAP tilings, which we used to approximate 2D epithelial cell geometry, are closely related to the familiar Voronoi tessellations of the plane. By virtue of its duality to a Delaunay triangulation, a Voronoi tiling satisfies the geometric compatibility constraint given by Eq.
(3). Equilibrium CAP tilings generalize the Voronoi construction to circular arcs instead of straight edges and correspondingly have a non-planar dual triangulation (see Appendix A. for the explicit construction). In the limit of constant pressure equilibrium CAP tilings reduce to standard Voronoi polygons with one additional degree of freedom per cell associated with the "isogonal" deformation previously decribed in [17].
In addition to testing the validity of the GVM-based stress inference, our analysis of myosin and inferred stress distributions in Drosophila embryo has revealed that despite the dynamical nature of GBE, the epithelial shell of the embryo maintains approximate mechanical equilibrium, in the sense that mechanical stress associated with the observed myosin distribution is mostly (at the 80% level) balanced internally and does not contribute to cellular flow. This conclusion is reached by a direct analysis of the measured myosin tensor, Eq. (15). Quite remarkably, the presence of this balanced internal stress is also correctly inferred from the GVM-based analysis of cell geometry across the surface of the embryo. The conclusion that tissue flow coexists with approximate internal force balance within a rearranging array of cells, provides an interesting insight into mechanics of tissues.
Also of biological interest in the analysis above, is the quantitative mapping of the unbalanced myosin, which according to [28] acts as a driver of global morphogenetic flow. Disregarding the errors at the anterior and posterior poles that arise from image resolution issues, the largest deviation from mechanical balance is found (Fig. 5 ABC) along the dorsal surface, where total myosin falls below the level needed to for internal force balance. This suggests an important role for DV patterning in the convergent extension flow. We believe the ability to estimate global patterns of balanced and unbalanced stress on arbitrary two-dimensional surfaces opens up a novel method in which one can identify the factors that drive morphological change.
We expect the GVM-based stress inference to be immediately useful for experimentalists studying tissue mechanics and the mechanics of morphogenesis of entire organs.

Materials and Methods
Confocal microscopy Raw data shown in Fig. 4 was taken on a Leica SP8 confocal microscope equipped with two HyD detectors, a 40x / NA 1.1 water immersion objective, and 561 nm laser line.
Light sheet imaging In toto images for Figs 5 and 6 where taken on a custom-built multi view light sheet microscope described in [39]. Briefly, the setup consisted of two excitation and two detection arms. On each detection ∆θ = θ stressθ div  [35]) and the principal (extension) axis of inferred local stress (using different methods). Stress for mitotic cells was estimated 2 minutes before the registered time of division. While GVM-based analysis finds significant correlation, predictions based on the extended "matrix inverse" [25,27] [40], and custom written MATLAB code. Fusion of individual views taken at 45 degrees angles was carried out using FIJI multi-view fu-sion plugins [41]. Cartographic projections where generated using ImSAnE [34].
Fly stocks Sqh-GFP; membrane-mCherry Numerical solution of Eq. (14) A 2D triangulated mesh of the embryo surface was constructed using ImSaNe [34] and FELICITY [42] -a finite element software package for MATLAB -was utilized to solve the PDE and compute surface derivatives. We refer the reader to [28] for a detailed description of the method.
Image segmentation In-vivo data was segmented using a custom pipeline implemented in MATLAB, available on Github at nnoll/TissueAnalysisSuite. ilastik, a supervised machine learning classifier, was used to as a pre-processing step for each image, followed by the application of a Laplacian of Gaussian filter. The resultant image was segmented using the watershed algorithm.
Vertices were defined as branch points of the resultant skeletonization -edges are segmented as the set of boundary pixels that run between two such branch points. Our CAP is parameterized by not only vertex position but also edge curvature. Each edge was fit to a circular arc using the Pratt method, which is robust for small angle samplings of the underlying circle. Interfacial myosin concentration was measured by dilation of each segmented edge by 2 pixels and averaging over the resulting set of pixels in the myosin channel. All segmentation information is stored within a custom data structure and can be immediately used for the GVM inference.
Cell divisions during late germ band extension were registered by tracking cells. Cell tracking was achieved by computing pixel overlaps between segmented cells in subsequent time points -cells were paired based upon the cell they most overlap with in the succeeding frame. Mitotic cells were defined as tracking events where two cells overlapped with one in the previous time point. The tracking was manually curated to ensure no false divisions were called.
3D reconstruction ImSaNe [34] was used to measure, parameterize, and store the surface and embedding coordinates of the Drosophila embryonic surface. Segmentation of cells was done using the cylindrical mapping of the embyro. The 3D vertex positions were subsequently estimated using the embedding grids obtained from the ImSaNe algorithm.
We computed the mesoscopic myosin distribution in the same way as detailed in [28]. The output from the automated segmentation of myosin is a summation over microscopic nematic tensors of the form subject to the constraint that α p α z 2 α = 1. R αβ is segmented from the image. To test the robustness of the GVM inference, synthetic data was generated by initializing a triangular lattice of ∼ 120 generating points q α within a rectangle of size [1, √ 3/2]. Using the Delaunay/Voronoi correspondence outlined in Appendix A, one can easily generate arbitrary CAP networks -i.e. for any set {q α , p α , z α } -in mechanical equilibrium. Edges within the CAP network are defined by Eq. (A2) and thus d 2 α must be computed from {q α }. Distance from q α is calculated using MATLAB's bwdist function and then scaled according to Eq. (A1) to obtain d 2 α (r). This procedure is repeated for each point q. The minimum of {d 2 α } for each spatial location is taken -the net result is a scalar field that measures the minimum weighted distance away from any triangulation vertex, d 2 (r) = min α d 2 α (r). An example of the above procedure is shown in Fig. A.1A. Edges are 'ridges' representing local maxima d 2 (r) that are found easily using the watershed algorithm - Fig. A.1B shows an example. The resultant equilibrium network can then be immediately calculated from the original parameters for {q α }, {p α }, and {z α }.
Intercellular pressures p α were sampled from a uniform distribution to generate lattices of varying curvature. Additionally, values for {z α } were pulled from a Gaussian with standard deviation .2. White noise, denoted δ in Fig.A.2, was added to both vertex position r i as well as position of edge centroids ρ αβ in order to assay sensitivity of the algorithm.
We benchmarked the relative efficacies of four mechanical inverse algorithms as a function of measurement noise and the contribution of pressure to the mechanical balance: (A) the matrix inverse with constant pressure [25], (B) the extended matrix inverse of [27], (C) a variant of the GVM inference constrained to constant cell pressures, (D) the GVM algorithm. The resultant correlation between inferred and the known generated tensions and pressures are shown in Fig. A.2(A-D). Both the variational constant pressure (C) and GVM inferences (D) are significantly less sensitive to measurement noise than their matrix inverse counterparts (algorithms (A) and (B) respectively). Furthermore, we see algorithm (B) loses accuracy at moderate curvature values (shown on the vertical axis, normalized to the average chord length between vertices) due to the assumptions of small pressure differentials required to linearize the underlying force-balance equations. Conversely, GVM is robust to both, showing high correlation in the entire tested parameter regime.
Appendix D: GVM algorithm applied to curved surfaces As proposed, the GVM can be easily extended to formulate a tractable inference scheme for the balanced mechanical stress within a curved tissue's tangent plane. Due to the inclusion of edge curvature information, the inverse is extensively over-determined which allows one to simultaneously infer both bulk and boundary stress using information of just bulk geometry. Consequently, the global mechanical state can be 'stitched' together by inferring stress on local patches of cells that can, with good approximation, be treated as planar. The blue cell array depicted in Fig. A  All cells are segmented in 3D using ImSaNe [34] Each local cellular patch is projected from the 3D embedding space onto the 2D plane that minimizes the sum of squares of deviation. The resulting projection is used as an input into the GVM method. (C) Scatter plot between the inferred tensions using the workflow outlined and the known tensions shown in (A). The inset displays the dependence on the number of patches used to cover the sphere. As was expected, correlation is monotonic with sampling resolution.
vided the area of interest is much smaller than the surface's radius of curvature, we can fit a well defined tangent plane to the patch. Let R n i denote the 3D position of the i th vertex within patch n -it is a matrix of size 3 x v n where v n is the number of vertices contained in patch n. The best fit triad of vectors is obtained easily via an SVD decomposition R n i = U n ΣV T n The approximate planar graph of patch n, shown as the black cell array in Fig. A.3B, is obtained by projecting R n i onto the two principal components; GVM is applied on the distortion. Importantly, the set of inferred tensions and pressures within each patch is unique up to an overall scale and thus there exists an undetermined relative scale between each patch -denoted λ n . We fix such scales by definining each patch to overlap by 1/4 their linear extent such that a subset of edges are involved in multiple inferred regions. Hence, the scale λ n for each patch is found by minimizing the squared difference between inferred tensions of edges shared by adjacent patches globally, subject to the constraint that the average scale is 1 to ensure a non-trivial solution λ n = 0.
This was used to define the patch size used in the empirical measurements during Drosophila embryogenesis. The outlined procedure was validated in-silico for synthetic spherical embryos containing roughly 3000 cells, with mechanics patterned by a vertex model minimized on the surface of a sphere. An example of a simulated embryo with azimuthal pattern of tension is shown in Fig.A.3A, both in the embedding space and in the cylindrical unwrapping of the sphere, analogous to the data shown for the Drosophila embryo. As shown in Fig. A.3C, excellent agreement between the inferred and known tensions was found provided the patch size was small compared to surface curvature. The inset shows this occurred when the defined patches contained 100 or less cells. This was used to define the patch size used in the empirical measurements during Drosophila embryogenesis.