Siqi Liu

Patrick Lemaire

Edwin Munro

Madhav Mani
email: madhav.mani@gmail.com

A mathematical theory for the mechanics of three-dimensional cellular aggregates reveals the mechanical atlas for Ascidian embryogenesis

The intricate three-dimensional (3D) structures of multicellular organisms emerge through genetically encoded spatio-temporal patterns of mechanical stress. Cell atlases of gene expression during embryogenesis are now available, but connecting these to mechanical design principles that govern the emergence of embryonic shape requires an ability to measure mechanical stresses at single-cell resolution, across embryos, over time -A Mechanical Atlas. Developing a new mathematical theory for the static mechanics of 3D multicellular aggregates involving cell pressures and membrane and line tensions, we present a parameter-free and image-based strategy to constructing spatio-temporal maps of the mechanical stresses driving morphogenesis. We present in-silico and in-vivo evidence in favor of the accuracy and robustness of our approach. The ensuing mechanical atlas, within the context of ascidian gastrulation, reveals the adiabatic nature of the dynamics, its dependencies on the cell-cycle and cell-lineage, and the novel identification of spatio-temporal variations of cellular pressures.

Introduction

A central challenge in biology is to understand how embryonic cells collectively shape tissues, organs and whole embryos to determine organismal morphology. At the heart of morphogenesis are forces that are produced locally by individual cells, transmitted across cell-cell contacts and resolved as global patterns of cell shape change, cell movement, and cell rearrangement. In recent decades, we have learned a great deal about the molecular machinery that cells use to produce and transmit forces within developing embryos. We have also begun to learn how the deployment of this machinery is controlled in space and time by genetically programmed and self-organized patterns of gene expression and intercellular signaling. However, ultimately, the biochemical signals that "control" morphogenesis, and the molecular machinery they control, must speak the common language of mechanics, by controlling the local production or transmission of active forces, or by controlling local material properties that resist deformation or displacement of cells. In addition, there is increasing evidence that developmental and physiological control mechanisms read the local mechanical force and deformation within cells and tissues and transduce these into biochemical signals, to close mechanochemical feedback loops. Thus the ability to measure mechanical forces over space and time in developing embryos is an essential prerequisite to understanding fundamental design principles for robust multicellular tissue morphogenesis.

In response to this challenge, a variety of methods have been developed to measure forces in living embryos. One general class of methods involves local application of force (e.g. through local mechanical indentation (ref) or micropipette aspiration (ref) or optical or magnetic tweezers (ref)) or local disruption of mechanical continuity (e.g by laser ablation (ref)), coupled to measurement of the resulting deformations. A complementary class of methods involves the observations of "force sensors" (e.g. molecular fret sensors, liquid droplets or elastic beads), embedded within a force-producing tissue, whose deformations reflect internally generated forces. Importantly, all of these methods require the use of physical models to infer forces or material properties from observed deformations. These methods have provided valuable insights into the mechanics of cell and tissue morphogenesis at different spatial and temporal scales. However, most are invasive and/or allow only sparse sampling at a only few locations at a time, and this has limits their use in mapping spatiotemporal patterns of force that underlie the collective mechanics of multicellular tissue morphogenesis.

Recent advances in our abilities to image cell and tissue morphologies and deformations in living embryos has spurred the development of a third class of methods, commonly known as image-based force inference. Force inference methods use embryonic cells themselves sensors, using microscopic observations of cell shape and deformation to infer the underlying forces. Like the other approaches, image-based force inference relies on physical models that predict cell shape and deformation, given specific assumptions about the nature, organization and magnitudes of the forces operating within individual cells and across cell boundaries. The central idea is to solve an inverse mechanics problem to build a mapping from the observables in images of living tissues -the geometries of cells and cell-cell contacts -to the hidden variables such as tension along cell contacts and intracellular pressures, that define the mechanical state of a tissue. A key virtue of image-based approaches is that they are non-invasive, and by construction, offer the possibility of inferring forces at many simultaneously observed points within an embryo over time.

Thus far, with a few recent exceptions (ref), methods for force inference have been implemented in two dimensions. Their applications have focused on the analysis of tissues, such as planar epithelia, in which it can reasonably be assumed that the dominant mechanics operate in 2D, or which have sufficiently stereotyped geometries that a third dimension can be absorbed into a 2D description. However, it has become increasingly clear that embryos regulate cell and tissue morphology in 3D -that is the relevant forces are both patterned and resolved in three dimensions, and in these many cases, restriction to 2D fails to capture the relevant mechanical variables. For example, the physical conjugate to hydrostatic pressure -the cell volume -has no physically meaningful counterpart in two dimensions. Likewise, in most embryonic tissues, active contractile forces are organized with respect to surfaces or lines of contact (between 2 and 3 cells respectively) which occupy and explore all possible orientations in 3D, defying a simple 2D representation.

Given these limitations of 2D force inference, and rapid advances in imaging approaches that allow to capture high resolution data on cell boundaries in 3D over time in living embryos, the need to expand the capability for force inference to 3D has become increasingly acute. However, doing so presents a number of serious challenges. First, moving to 3D reduces the information content of the input data relative to the number of quantities that must be determined to infer the mechanical state. Second, existing force inference methods rely on an image segmentation step in which an intermediate description of cell boundaries is extracted from the raw data. Unfortunately, the accuracy of the force inference is highly sensitive to the quality of the segmentation, and the segmentation becomes increasingly difficult in 3D.

Here we describe the formulation, implementation and application of a robust approach to 3D force inference which overcomes these limitations. Our approach is based on a new physical theory for the mechanics of close-packed three-dimensional cellular aggregates in which cell geometries are governed by the balance of hydrostatic pressure forces, and contractile forces operating (differently) on curved surfaces and curved lines of contact between adjacent cells. We derive a mapping from the set of possible equilibrium geometries to values for pressure, surface tensions and line tensions. This mapping is unique up to three scalar values that define absolute scales for pressure, tension, and a "zero-mode" that describes the relative contributions of line tensions and surface tensions to defining a given geometry. We implement a method to project raw image data produced by light sheet imaging of membrane-labelled specimens directly onto the nearest possible equilibrium geometry, without interposing a segmentation step. Applying our approach to synthetic image data representing cellular aggregates at mechanical equilibrium, we show that the global nature of the projection and the underlying mechanics ensures a robust solution, recovering the correct equilibrium geometry, with < 5% error for typical signal:noise ratios, allowing robust inference of the underlying pressures and tensions. Finally, we apply our approach to construct a mechanical atlas for gastrulation in the early ascidian embryo from light sheet data collected from membrane labelled ascidian embryos at multiple timesteps during gastrulation. In addition to confirming some previous predictions about the forces that drive, the mechanical atlas reveals surprising new information, including lineage specific temporal control over intracellular pressure, and the role of mechanical integration of pressure and tension across multiple tissues in driving endoderm invagination (the first major step of gastrulation). Together our results introduce a robust new tool that be applied to infer global mechanics in 3D in many other organismal and embryonic contexts, and highlights the power of this tool to reveal fundamental physical design principles for morphogenesis

2/32
F F F αβγ + F F F βγδ + F F F αβδ + F F F αγδ = 0.

Theory

Overview of Theory

Biologically, there are 3 cell biological factors that we incorporate into our three-dimension static mechanical theory that have broad consensus in the community. These three factors are schematized in Figure 1 A-C. 1) Effective cell pressures that can vary across cells and over time. These pressures are considered to be an effective stress that are isotropic, i.e. lack any directionality, impacting a cell's surrounding in a democratic manner. Phenomenologically speaking, this effective cell pressure produces a curvature of cell-cell contacts. The origins of the pressure itself have contributions arising from hydrostatic effects and due to isotropic effects of active cellular processes, including the activity of the actin and microtubule based systems in a cell. In our model, cells have effective pressures that can vary across cells and over time. Contrasting the two-dimensional theory, the conjugate variable to the effective pressures we infer are cell volumes, not cell areas, thus accurately accounting for any hydrostatic contributions to cell pressure. 2) Effective surface tensions along surfaces of cell-cell contact. This is incorporated in our model based on the overwhelming evidence that the contractile activity of the acto-myosin systems that are relevant to morphogenesis in epithelial-like tissues are localized and regulated at cell-cell contacts. There are an innumerable number of spatio-temporal heterogeneities and anisotropies within a single surface of contact that our model averages into a single number. Rapid timescale dynamics that probe the elastic features of the acto-myosin system and the cell cortex are also ignored. Assuming a viscoelastic model, wherein the long time dynamics are fluid-like, our model incorporates all the manifest biological complexity into a single number representing a homogeneous and isotropic effective surface tension. This number varies across cell-cell contacts in a multicellular aggregate and over time. Contrasting the two-dimensional theory, the surface tensions we account for correspond to two-dimensional surfaces that are curved within a three-dimensional aggregate of cells. 3) Effective line tensions along the curves along which three cells meet. While a theory for foams comprised of soap bubbles need not incorporate such an effect, it must be incorporated into a three-dimensional mechanical theory of an embryo. Following studies in the ascidian embryo itself, and other systems, there is considerable evidence that the contractility of the acto-myosin system along the triple-contact of three cells is independently regulated, and thus we account for it. Again, our theory ignores the potential heterogeneities along curves of triple-contact, and again we restrict our attention to only the long time dynamics where the response of the system is anticipated to be fluid-like. Our theory accounts for a single homogeneous effective line tension, which can vary across the triple-contact curves in a three-dimensional multicellular aggregate, and over time.

3/32

Physically, we incorporate two factors based on significant experimental evidence in embryos and multicellular aggregates, and our fundamental understanding of mechanics. 1) The three biologically motivated mechanical ingredients incorporated into our theory, of effective cell pressures and tensions, relate to each other according to the Young-Laplace (YL) relation. Intuitively, YL relation posits that any pressure difference between the two contacting cells produces a force which will bend, or curve, the surface of contact. The extent of the curvature of the surface of contact is a function of the pressure difference and the effective surface tension along it. A surface with a high surface tension will be less curved due to a given pressure difference than a surface with low surface tension. The YL relation identifies that the ratio of the pressure difference and effective surface tension determines the mean curvature of the surface of contact. Since our theory only accounts for a single number accounting for a homogeneous and isotropic surface tension, the YL relation determines that the only surfaces of contact with a constant mean curvature, i.e. sections of a sphere, are possible within the scope of the theory. In an entirely mathematically and physically equivalent manner, now considering a curve of triple-contact, the three surfaces pull the curve thereby bending it. The extent of its curvature is a function of the three surface tensions bending it, and the line tension along it, that attempts to keep it straight. A second YL relation relates these two forces to each other, determining the curvature of the curve of triple-contact. Since our theory only accounts for a homogeneous line tension, the YL relation determines that only curves of constant curvature, sections of circles, can arise within our theory. Thus, the only geometries that can arise from our theory are spherically-curved-polyhedral (SCP) tessellations -those whose surfaces are sections of spheres and lines are sections of circles. 2) The incorporated mechanical ingredients in our model are close to a static balance with each other. This assumption is founded on significant evidence emerging from laser ablation experiments. Generically, following a laser ablation, the tissue immediately recoils on timescales on the orders of a few seconds. This timescale is orders of magnitude faster than the timescales of morphogenesis itself, which occurs generically is on order of several minutes to hours. This separation of timescale suggests that the overwhelming fraction of mechanical forces are in balance with each other. The observed dynamics of multicellular aggregates in embryos are therefore not due to a mismatch of the mechanical factors but due to slow changes in the mechanical parameters themselves. This physical ingredient is incorporated into our theory through requiring that the line tensions associated with curves of triple-contact balance each other at vertices where they meet.

The simplicity of the theory outlined above is a virtue since it imparts a kind of rigidity. As described above, the theory can only give rise to geometries where surfaces of cell-cell contact and curves of triple-contact are sections of spheres and circles, respectively -a SCP tessellation. This feature is a non-negotiable of the theory. The quality of the applicability, or fit, of the theory to data from a particular system can thus be quantitatively assessed. Significant departures from the theory indicate its over simplification in a given empirical setting, and the structure of the departures indicate the kind of additional mechanical complexity that must be included. We demonstrate the unreasonable quality of the fit of the theory to live-imaging data during ascidian gastrulation, thus given credibility to the inferred spatio-temporal patterns of effective cell pressures and tensions.

The Fundamental Equations of the Theory

This section presents three equations that mathematizes the biological and physical ingredients discussed above at an intuitive level. Based on the above assumptions we are assured a static mechanical balance between effective cellular pressures P α , surface tensions T αβ , and line tension F αβγ . Indices label cells, where membranes and edges are indexed by the two and three cells that define them, respectively. The balance of stresses normal to membranes (Figure 1A) is given by the Young-Laplace equation,

|P α -P β | = 2T αβ H αβ , ( 1a
)
where H αβ is the mean curvature of membrane surface. Homogeneous pressures and tensions induce constant mean curvatures, enforcing that the shape of a membrane is spherical or planar (a sphere with infinite radius of curvature), as shown in Figure 1A. A cellular aggregate comprised of a tessellation of polyhedrons where all membranes are planar corresponds to situations with a uniform pressure in all cells P α = P 0 . When there are pressure differences, the faces of the polyhedral cells must be sections of spheres instead of flat planes, inducing the geometry of cellular aggregate as a spherically-curved-polyhedron (SCP) tessellation.

Membrane tensions act on a junction providing a traction force normal to the direction of the junction (Figure 1B). The junctional balance relation between them is a modified, one-dimensional, Young-Laplace equation:

T T T αβ + T T T βγ + T T T γα = -F αβγ κ κ κ, ( 1b
)
where κ κ κ is the curvature vector of the junction. In a SCP tessellation, the junction must be a circular arc since it is formed by the intersection of spheres. Quantitatively, κ κ κ has a constant magnitude with the direction towards the circle's center. Increasingly complex voronoi diagrams. A) 2D voronoi diagram, the blue dots are voronoi sites q q q, the black lattice are cell regions {r r r|d 2 α (r r r) ≤ d 2 j (r r r), ∀j = α}, under Euclidean distance d 2 α (r r r) = |r r r -q q q α | 2 . B) 2D weighted voronoi diagram, with same sites as A), cell regions under the weighted Euclidean distance d 2 α (r r r) = |r r r -q q q α | 2 -θ 2 α . Size of blue dots represents the weight θ. C) 2D generalized weighted voronoi tessellation, with same sites and weights as B), cell regions under the distance d 2 α (r r r) = p α (|r r r -q q q α | 2 -θ 2 α ). The color scale of the blue dot represents the power p.

A vertex is defined as a point where junctions corresponding to adjacent cells meet (Figure 1C). While higher order vertices are naturally accounted for by our model we flesh it out in the context of 4 intersecting junctions. The line tensions along these four junctions, as defined above, balance at the vertex, satisfying

F F F αβγ + F F F βγδ + F F F αβδ + F F F αγδ = 0. ( 1c
)
Equations ( 1a)-(1c) together determine the correspondence between effective cellular pressures, surface tensions, and line tensions, and the geometry of the SCP tessellation.

An Inference of Mechanics from Geometry

Overview: The above equations, (1a)-(1c), mathematical embody our biological and physical assumptions. Each equation has mechanical aspects pertaining to the numbers corresponding to effective pressures and tensions, and numbers associated with the geometry of various interfaces of contact between cells in three-dimensions. These equations suggest the possibility of inferring the mechanical parameters of a multicellular configuration from imaging of its geometry alone. If possible, we would be able to infer the spatio-temporal patterns of effective cell pressures and tensions! But is it possible -given the geometry of cells in an embryo, can we infer all its effective pressures and tensions? The problem cast in this manner is a classic example of an inverse problem. At first glance the inverse problem seems trivial since equations (1a)-(1c) are linear! However, the compact index notation in these equations hides the interdependence of these relations in a compact cellular aggregate. In particular, the same quantity appears in multiple relations, making the system of equations highly coupled. In fact, a priori, it is not even clear whether the highly coupled systems of equations are solvable -are there a sufficient number of constraints for the number of mechanical parameters that we wish to solve for? Though the fundamental equations, 1a-1c, are local balance relations, solving the problem requires thinking globally. In light of the fact that the only permissible geometries are SCP tessellations we have additional knowledge of nonlocality of the constraints. For example, the curvature of a line in the multicellular aggregate and its two vertices aren't independent since they must define a section of a circle, while the lines that circumscribe the boundary 5/32 of a surface and its curvature can't be independent since they must collectively reside on a section of a sphere. Thus, solving the inverse problem requires a reparameterization of the geometry that naturally admits these global constraints. We therefore need only search for solutions to equations (1a)-(1c) within the space of SCP tessellations, rather than the far more vast space of all multicellular geometries! Thus, the protocol for inferring forces from a real image is to 1) determine the closest fitting SCP tessellation to the observed geometry, and 2) to determine its corresponding mechanical configuration.

Step 1 requires a quantitative parameterization of the entire space of SCP tessellations. Fortunately, a generalization of Voronoi diagrams is sufficient. In particular, the generalized weighted Voronoi (GWV) construction (Figure 2) provides an exhaustive accounting of all SCP tessellations. Step 2 requires an analytical bridge between the parameters of a GWV diagram and its effective pressures and tensions. This requires an analytical solution of the original equations of motion, which is one of the central results reported in this study. We further demonstrate that while the system of equations (equations (1a)-(1c)) are solvable, they do not admit a one-to-one map from the parameters of the GWV diagram to its mechanical state. Said another way, distinct mechanical configurations correspond to the same geometric configuration. In physics these are referred to as zero modes of a system, and as such aren't constrained by the data at hand. Quantitatively speaking, we identify 3 global scalar zero modes in our theory that require information beyond the observed geometry of a mutlicellular aggregate.

Generalized weighted Voronoi constructions parameterize SCP tessellations: SCP tessellations are comprised of vertices connected by lines that are sections of circles, which are themselves connected by curved membranes that are sections of spheres. We define each membranes centroid to be the 3D coordinate of the center of the sphere that best approximates it. Manifestly, vertex and membrane centroids are inter-dependent. For instance, the vertices of a given membrane must be equidistant from its centroid. As for polyhedral tessellations with planar membranes, vertex locations of one membrane are constrained to lie on a plane -ensuring that they cannot be independent of each other. Additionally, because each junction is the intersection of three membranes, the three corresponding membrane centroids must be compatible with each other. In the SI (Part I.B), we demonstrate that the centroids of three adjacent membranes must be co-linear. Summarizing, the parameters of a SCP tessellation are constrained to avoid geometric conflicts. Accounting for these constraints and conflicts we prove in the SI (Part I) that, given, C, the number of cells in an aggregate , an SCP tessellation has 5C independent degrees of freedom (4C for polyhedral tessellation -An SCP tessellation with flat surfaces of contact and straight lines of triple contact). Following this, we demonstrate that we can parameterize an SCP tessellation using a generalized weighted Voronoi (GWV) diagram (Figure 2), which is defined by 5C independent parameters Ψ = {q q q α , θ α , p α } α=1,2,...,C . q q q α is the original site that generates a three-dimensional cell region, θ α and p α are the weight and the power of the distance definition (see SI Part I.C for a description of the generalized weighted Voronoi construction). Loosely speaking, the generators of the diagram identify the location of cells, while the weights and powers determine their size and the curvature of cell-cell contacts. Polyhedral tessellations are a special case of SCP tessellations where all the powers are identical, p α = p 0 .

An analytical solution to the inverse mechanical problem:

The GWV parametrization provides a minimal, but exhaustive characterization of SCP tessellations. To solve the inverse mechanical problems requires deriving an analytical solution of the equations of motion, (1a)-(1c), that relate the parameters of a GWV construction to its effective pressures and tensions.

Before we present our results in three-dimension, we first describe how our construction works in two-dimensions. The derivation of the two-dimensional results is distinct to Noll et al's approach. [2]. The novel approach to deriving the two-dimensional results makes clear the path to the three-dimensional generalization. We will heavily rely on a few elementary geometric facts. Central amongst them being the concept of a "vector area", which is a vector S S S = S n n n with the face area S as the magnitude and the face normal n n n as the direction -see Figure 3A. For a curved face, the vector area is the integral of face elements S S S = n n ndS. A consequence of the divergence theorem, derived in the SI (Part II.A), is that the sum of all face vector areas of a compact closed shape must be zero: i S S S i = 0 0 0.

Let us now consider the two-dimensional scenario with non-uniform pressures. We remind the reader that a two-dimensional system with cell-pressures and membrane tensions at static equilibrium must correspond to a Circular Arced Polygon (CAP) tessellation (the two-dimensional analogue of SCPs), as shown in black in Figure 3B. These can be parameterized by a 2d GWV tessellation, taking {q q q α , θ α , p α } α=1,2,...,C as parameters. Before introducing any mechanics we must introduce the global geometry of a GWV tessellation in two-dimensions. Following the schematic shown in Figure 3B, for any point r r r on the membrane or vertex of cell α, we define the local dual point as q q q * α (r r r) = r r r + p α (q q q α -r r r). A corresponding dual point exists in the adjacent cell, q q q * β (r r r). The two local dual points corresponding to the point r r r on the membrane form the local dual line q * αβ corresponding to the membrane between cell α and β, with length:

q * αβ = q q q * α -q q q * β = (p α -p β )(r r r -ρ ρ ρ αβ ) = |p α -p β |R αβ .
Since the length is independent of The blue dots are voronoi sites q q q, and the red dots are local dual points q q q * with respect to points r r r on membranes or vertices. The dark red triangles are dual to the vertices, the light red curved quadrilaterals are dual to edges. (C) The divergence theorem of local dual graph reveals mechanical balance in 2D. Membrane tension balance T T T αβ + T T T βγ + T T T γα = 0 corresponds to triangle divergence theorem q * αβ n n n 1 + q * βγ n n n 2 + q * γα n n n 3 = 0 at vertex r r r i . Pressure-tension balance

T αβ n n n 1 + T αβ n n n 2 -P α dr n n n 3 -P β drn n n 4 = 0 corresponds to quadrilateral divergence theorem q * αβ n n n 1 + q * αβ n n n 2 + dq * α n n n 3 + dq * β n
n n 4 = 0 of edge element dr. (D) A two-cell example of the mechanical zero mode. In the upper case, no line tension results in three equal membrane tensions. In the lower case, line tension results in smaller middle membrane tension than two apical tensions.

the location along the membrane it implies that the length is a constant along the membrane. Additionally, the local dual line is parallel to the radial direction and thus perpendicular to membrane at r r r. The three local dual lines defined at a vertex form a local dual triangle, centered at a vertex. Furthermore, along the membrane, the local dual lines sweep out a curved quadrilateral. These constructions are shown in red in Figure 3B. This is the geometry of a generalized Voronoi construction for a CAP tessellation.

We are now positioned to apply the divergence theorem and connect the parameters of the Voronoi construction to mechanics in the 2D scenario. To help follow the mathematical procedure we have schematically excised the triangle and quadrilateral and displayed them in Figure 3C, along with the force balance relations they correspond to. Applying the two-dimensional divergence theorem to a local dual triangle, gives us q * αβ n n n 1 + q * βγ n n n 2 + q * γα n n n 3 = 0. We notice that the three vectors have same directions as the corresponding membrane tension balance T αβ n n n 1 + T βγ n n n 2 + T γα n n n 3 = 0. Therefore, the length of local dual line is proportional to the membrane tension: T αβ = cq * αβ , where c is a single constant for all membranes in a tessellation.

We now leverage the divergence theorem yet again, but now along the quadrilateral spanning membranes to help identify the non-uniform pressures from the underlying Voronoi construction (see Figure 3C). Considering the edge element dr r r along the membrane adjacent to cell's α & β, the dual graph is a curved quadrilateral as shown in Figure 3C. Applying again the two-dimensional divergence theorem to the quadrilateral, we obtain

q * αβ n n n 1 + q * αβ n n n 2 + dq * α n n n 3 + dq * β n n n 4 = 0.
Observing that the vectors have same directions as the corresponding tension-

-P β = -c dq * α -dq * β dr = c(p α -p β ).
Finally giving as the pressure in cell α, P α = cp α + b, where c is the same constant as before, and b is a single global constant setting the overall scale of pressures.

The extension and derivation of this machinery to three-dimensional situations can be found in the SI (Part II.C). As the reader no doubt anticipates, the additional complexity associated with three-dimensions is now the presence of vertices, lines, membranes, and cells, as opposed only the vertices, lines, and cells in two-dimensions. This additional topological and mechanical complexity results in one additional application of the divergence theorem. As expected, this additional using of the divergence theorem results in one additional bridge between the GWV representation and its mechanical state, introducing one additional unknown global constant. We only present the final results here, and its detailed derivation and intuition can be found in the SI (Part II.C),

F αβγ = cS * αβγ , T αβ = 1 2 (b -c(p α + p β ))q * αβ , P α = cp 2 α -bp α + a.
(2)

These three equations are the most central result of the paper and a solution to the inverse mechanics problem. They are what make possible revealing and embryo's mechanical atlas at single cell resolution over time. We highlight that the map is nonlinear and in no way could one intuitively derive this based on simplistic mathematical and physical considerations.

Zero modes: We notice that there are three undefined constants a, b and c in the expressions of forces. This indicates that for a given SCP tessellation, there are three zero modes in the correspondence from the mechanical degrees of freedom admitted in our model to a given geometry [3]. These are three global scalar quantities that determine the overall relative contributions of pressure, membrane tension, and line tension to the observed geometry. Said conversely, there exists a mechanical degeneracy wherein a single SCP tessellation corresponds to a three parameter family of mechanical states. The constant a refers to the background pressure, which can be set to 0 without loss of generality. Physically speaking, this freedom is permitted by the fact that only changes in pressures enter our mechanical model, and not the absolute scale of pressures themselves. The constant b is most simply understood through considering a simpler system devoid of line tensions. In this case, the overall constant b determines the overall scale of forces in physical units. Manifestly, an inference based purely on geometry can never inform as to whether the tension on a membrane is on the pico vs nano Newton scales. This constant is the same as in the two dimensional theory that precludes us from determining the physical units mechanical stresses in the system. The constant c is truly a three-dimensional feature that the presence of line and membrane tensions gives rise to.

What is the physical interpretation of the third zero mode, embodied in the value of the constant c? A physical demonstration of this zero mode can be seen in a simple two cell system shown in Figure 3D where the two cells have equal pressure, and the angles between three membranes are equal at the junction. In the absence of a line tension, the three membrane tensions are required to be in balance, and, thus, equal to each other. In the presence of a line tension, which provides an additional force towards the center of the middle membrane, to maintain the same geometry, the middle membrane tension must adopt a smaller value to accommodate the presence of the line tension and maintain static force balance.

Summarizing, while there are formally 3 zero modes only the parameter µ = c/b breaks the mapping from geometry to mechanics. Generalizing to systems more complex than a symmetric two-cell system, once we admit the possibility of an independent set of line tensions, in addition to membrane tensions, the same complex SCP geometry can be sustained by a one parameter family of mechanical states, that, loosely speaking, put more or less weight in the line tensions in the system. A detailed analysis of this can be found in the SI (Part II.D). We can deduce bounds on the possible values that µ can adopt by simply requiring that all the membrane and line tensions in the system are pulling/contractile forces. This requirement is biologically motivated by the fact that the membrane and junction-localized cytoskeletal elements generate contractile forces, and generically do not sustain loads. In particular, the value of µ is bounded by 0 ≤ µ ≤ µ 0 , where µ 0 = min α,β (p α + p β ) -1 . Derivation of this result can be found in the SI (Part II.D). Manifestly, to determine the value of the parameter µ requires information beyond the geometric data this current investigation is based on. We are thus limited to have an inference of mechanical stresses in an embryo up to a single global constant.

Fitting the data

Overview:

The above relations provide the map from parameters of a weighted Voronoi diagram of a 3D multicellular aggregate to its mechanical state -up to zero modes that cannot be identified by any image-based approach. What is thus required is to assess whether the observed, almost certainly more complex, geometry of the embryo can be well approximated by the simpler geometry of a GWV diagram. This global projection of the data onto a simpler class of tessellations is done numerically, giving us a quantitative handle on the degree of mismatch between the observed geometry of the embryo and its closest fitting GWV diagram. Thus, not only does the scheme proposed fit the model to data, it also provides a quantitative assessment of the quality of the fit both in space and time. Below we discuss the quality of the fitting procedure, over space and time, which gives strong evidence in favors of our theory's assumptions and the quality of our fitting scheme.

A minimization statement: Quantitatively speaking, the segmented image obtained from live embryo imaging provides the location of membrane pixels, r r r, from which the geometry needs to be reconstructed. Denoting the centroid and radius of membrane sphere by ρ ρ ρ i and R i , the deviation of the pixel from sphere is i = |r r r i -ρ ρ ρ i | -R i (Figure 4D). Since ρ ρ ρ i and R i are both functions of Ψ, we can recover the geometric parameters by least-square fitting to a SCP tessellation. Specifically, we minimize the mean-squared-deviation (MSD) function,

E(Ψ) = 1 N N i=1 [|r r r i -ρ ρ ρ i (Ψ)| -R i (Ψ)] 2 , ( 3
)
where N denotes the total number of pixels, to find the 'closest' SCP tessellation to the observed data. Our minimization statement falls into the class of nonlinear optimization problems, requiring an educated initial guess for the desired parameters that takes advantage of the properties of GWV diagrams -see SI (Part III.B) for further details. The minimized value of E(Ψ) is a global measure of the average deviation of the empirically observed geometry and the fit SCP tessellation. However, as we will demonstrate below, our approach provides a finer-grained spatial information pertaining to the errors in our approximation.

Precision and robustness of data fitting: Prior to discussing results from the ascidian embryo we wish to validate and assess our entire schema in an in silico scenario where we have access to the ground-truth. Given 9/32 (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission.

random parameters Ψ 0 of C cells in a 3D region, a GVW diagram can be constructed from which the membrane pixels in the 3D image can be identified. After introducing Gaussian noise into all pixel locations, synthetic noisy images of cellular aggregates are generated, as shown in Figure 4A. Our fitting scheme recovers a Ψ and we can compare it with Ψ 0 , the ground truth. The results (Figure 4B) report the degree of mismatch, or error, between the parameters of the weighted Voronoi diagram inferred by our scheme, relative to ground truth, as a function of the degree of noise injected into the geometry. Inferring the different parameters of the weighted Voronoi display varying degrees of susceptibility to noise, however, the schema recovers parameters with up to 99% accuracy with as much as 20% noise injected into the location of membrane pixels. See SI (Part III.C) for more details.

To provide further evidence for the robustness of our numerical fitting scheme we perform a sensitivity analysis. Physically speaking, the sensitivity analysis we perform probes to what extent the inferred parameters of a Voronoi diagram jiggle when we jiggle the data around some configuration. Mathematically this involves performing an analysis of the eigenvalue distribution of the system's linear response [1]. A robust inference scheme is one with all its eigenvalue smaller in magnitude than 1, ensuring that no perturbations to the system can generate disproportionately large deviations in the values of inferred parameters. In particular, at the minimum of the MSD function, δE(Ψ, r r r) = 0, thus, for any r r r i , there must be

δ i = δ [|r r r i -ρ ρ ρ i (Ψ)| -R i (Ψ)] = 0.
This permits us to derive a local linear approximation, Kδr r r + M δΨ = 0, where K = ∂ /∂r r r and M = ∂ /∂Ψ are two matrices. As such, the relation between observed pixels perturbation δr r r and geometric parameter deviation δΨ is given by δΨ = Lδr r r = (-M -1 K)δr r r, where M -1 is the pseudoinverse matrix of M . For synthetic data, we compute the eigenvalues λ i of L and find they are all less than 1 in magnitude (Figure 4C). This indicates that the fitting scheme is robust to noise.

Results

Ascidian gastrulation is an approximately adiabatic process

As demonstrated above, our numerical scheme permits a quantitative evaluation of the quality or accuracy of the approximations inherent to our physical model and numerical scheme. The most severe of the three assumptions made in our theory is the dynamical one, suggesting that the multicellular aggregate is close to static equilibrium. Deviations from this would introduce strong velocity-dependent dissipative components to the stresses, resulting in geometries that would be poorly approximated by SCP tesselations. Additionally, any elastic contribution to the most salient constitutive properties of a membrane would produce significant deviations in membrane shapes away from the spherical expectations of our model.

To test the validity of our approximations, and confidence in our numerical procedure, we apply the fitting scheme on imaging data of ascidian gastrulation. In Figure 4E we report the distribution of errors, i , across all pixels of an embryo over an entire time course. The distribution can be well approximated as a Gaussian with a zero mean and a standard deviation of approximately 4 -5%.

We can also take a view of the spatial distribution of errors in the embryo (Figure 4F). There are clear deterministic deviations in the quality of the approximation. In particular, the posterior mesodermal lineages at the 64 cell stage (Figure 4F) display larger, and approximately left-right symmetric, errors across multiple embryos. As expected, the mechanical model and inference scheme is unsuitable for cells about to divide wherein the cell shape adopts a dumbbell-like shape just prior to the computational identification of a new membrane -See SI (Part III.E) for details. While present, the magnitude of errors is only around 5%. Much of the remainder of the spatial distribution of errors map in Figure 4F is random with only slight deviations from an SCP tessellation. The spatial location of errors is uncorrelated with faster moving regions of the embryo during the process of gastrulation, consistent with our adiabatic approximation.

As such, the quality of the fits that we are able to generate to live-imaging data of ascidian gastrulation provides quantitative evidence for the dynamical assumptions underlying our static mechanical theory. In particular, our results suggest that the embryo is very close to a static equilibrium configuration of forces. The observed dynamics can thus be considered the result of "adiabatic" changes in the configuration of the mechanical state, the morphogenetic timescale being much larger than the mechanical relaxational timescales of the system. More geometrically speaking, the dynamics of ascidian gastrulation can be very well approximated by dynamical trajectories that connect close by SCP tessellations. Despite the unfathomable complexities of sub-cellular, cellular, and tissue physics and biology manifest in this embryo, a highly general coarse-graining into just a handful of effective mechanical degrees of freedom can sufficiently account for the observed spatial and temporal variation in observed geometries. From a biological perspective, the observation that the entirety of the embryo dynamically maintains such proximity to a static equilibrium points to mechanisms that can globally coordinate the generation and regulation of forces.

A mechanical atlas for ascidian gastrulation

In light of the accuracy and robustness of our inference scheme, and of the quality of the model's fit to real images of the ascidian embryo we present our first mechanical atlas (Figure 5). Once we reconstruct the geometry and get the best fitting geometric parameters Ψ, the force maps can be visualized for any choice of µ -the undetermined zero mode in the system. While the quantitative details of the mechanical patterns depend on this undetermined parameter the courser-grained patterns remain invariant -see SI (Part III.E) for details. The mechanical patterns remain relatively constant in between rounds of cell divisions, so we take one time point in each stage for representation. Additionally, as shown in the SI (Part III.E), the reproducibility of the mechanical patterns across three separately imaged embryos is striking.

Our mechanical atlas makes manifest many of the known symmetries of the embryo. At this stage of ascidian development the embryo is symmetric about its left-right axis, and our mechanical atlases reflect this symmetry. Cross-sections of the embryo displaying cell pressures and membrane tensions at the 112-cell stage in Figure 5 clearly display large mechanical disparities along the animal-vegetal axis of the embryo. Beyond the morphological symmetries that are recapitulated in our mechanical atlases, it also reveals the mechanical differences between the three emerging germ layers -ectoderm, mesoderm, and endoderm. For example, a ring of high pressure (and high line tension) in the mesodermal lineages (starting along the posterior and sweeping anterior) assembles by the 76-cell stage, encompassing the invaginating endodermal layer. Additionally, the endodermal layer displays high surface and line tensions at the 112-cell stage as invagination ensues. The ectodermal lineage stays largely mechanically quiescent throughout gastrulation. Taken together, the quality of the model's fit to the data, its temporal robustness, reproducibility across embryos, the manifest morphological symmetries, and the distinct mechanical states of germ layers gives us confidence in the biological relevance, accuracy and robustness of our mechanical atlas of ascidian gastrulation.

Heterogeneities in cellular pressures during ascidian gastrulation: Here, we highlight patterns in cellular pressures, membrane surface-tensions, and junctional line-tension through the early phases of ascidian gastrulation that remain invariant for all values of µ. In particular, we show its force maps for µ = µ 0 /2 in Figure 5, and the force maps for other values of µ can be found in our SI (Part III.E).

In reference to Figure 5, at the very onset of gastrulation (64-cell stage) we consistently observe raised levels of pressures, apical surface tensions, and apical line tensions in the posterior mesoderm cells located in the vegetal hemisphere of the embryo.By contrast, the endodermal and ectodermal lineages lack significant mechanical heterogeneities at the 64-cell stage. While never observed before, our mechanical atlases make clear that the mechanical process of gastrulation initiates through mechanical activity at in the posterior mesodermal lineage of cells. By the 76-cell stage, all three mechanical readouts have increased in the posterior mesodermal lineage and spread to almost the entirety of the mesodermal lineage, which surrounds the mechanically less active endodermal cells. The ectodermal lineage displays few mechanical heterogeneities. By the 112-cell stage, the endodermal cells, which are in the process of invaginating, have raised mechanical activity. In particular, the cellular pressures, apical surface tensions, and line tensions, are now amongst the highest in the invaginating population of endodermal cells. The posterior population of mesodermal cells also displays raised levels of mechanical activity in all three mechanical measurements. The patterns reported here are observed across all permissible values of µ, and are robustly observed in the 3 embryos that we analyzed (see SI Part III).

Figure 5J-L shows a cross-sectional view of the embryo at the same 3 representative developmental stages. In line with studies focused on analysis the spatio-temporal distribution of myosin we too observe an increase in the baso-lateral membrane tensions of the endodermal population of cells, which is visualized in a distinct color map in Figure 5D-F. Additionally, we recover the progressive increase in apical junctional line tensions in the endodermal lineage through the process of its invagination, which is believed to be generated by a doubly phosphorylated population of the myosin II regulatory light chain.

However, we recover important and novel differences from the currently held view of ascidian gastrulation and embryonic morphogenesis more broadly. While the apical view of the pressure field in Figure 5A-C gives some sense of the heterogeneities in cellular pressures, it is the cross sectional images in Figure 5J-K that gives a the correspondence between the embryonic shape change and emergent pressure patterns. In particular, after the stabilization of a mesodermal pressure ring, all the non-ectodermal cells stabilize a pressure that is twice as much as what is measured in the ectodermal layer of cells. Even more striking is the emergence of a cage of membrane surface tension that encapsulates this high pressure region of the embryo. Together, the raised pressure and surface tension cage stabilizes a stiff layer of mesodermal and endodermal cells as they collectively invaginate, displacing a softer layer of ectodermal cells. An additional deviation from the current view of ascidian gastrulation is the increasing level of apical membrane surface tensions observed in the endodermal layer as it invaginates.

12/32 Figure 6. Lineage maps for gastrulation process in Phallusia mammillata embryo. A) The lineage map of inferred pressures of all cells from 64-cell stage through 112-cell stage. Each row of circles represent all cells at one time point. Cells are sorted by ecotoderm, mesoderm and endoderm, from anterior to posterior. B)The lineage map of inferred apical tensions of all cells from 64-cell stage through 112-cell stage. The undetermined zero-mode parameter µ is set to be µ = µ 0 /2. Inferred forces are normalized by the median of ectodermal pressures at each timepoint.

Our observations highlight the importance of developing new approaches to measuring the mechanical stress state of an embryo that go beyond a focus on the forces generated by the actomyosin cytoskeleton. Our global geometric approach reveals that you must account for significant patterned heterogeneities in cell pressures to explain the observed geometries through the process of gastrulation. It is worth emphasizing that our inference approach only takes in observed geometry and knows nothing about the temporal nature of the data, lineages, germ-layers, or the morphogenetic axes and symmetries of the embryo. Despite this, the measurement highlights smooth changes in mechanical patterns over time, the patterns themselves appear to respect lineage boundaries providing evidence of their biological significance, and displays the expected left-right symmetry. The observations as a whole, their spatio-temporal lineage-specific and symmetric patterns, give confidence in the accuracy and robustness of our mechanical atlas.

Lineage and cell-cycle dependent and independent features of the mechanical atlas: A virtue of the stereotyped lineage of ascidian embryos is that we can display the mechanical data on the invariant lineage maps. In Figure 6A we can see the time-dependent trajectories of pressures in cells. These maps have been sorted according to germ layers, and by the relative locations of cells along the anterior-posterior axis, within each germ layer. The lineage view of the mechanical data makes apparent the very strong lineage-dependent patterns, with only the mesoderm displaying an increase in pressures prior to the 112-cell stage, and the ectoderm never showing a significant change in its pressure distributions. Furthermore, at the 112-cell stage when the embryo begins to invaginate, there is a clear anterior-posterior gradient in cell pressures within the mesodermal layer of cells. Apical membrane tensions, shown in Figure 6B, show similar overall patterns, though the endodermal population does reveal more dynamic 13/32 patterns of apical surfaces than in its pressure field.

The lineage view of the mechanical data does also make manifest both cell-cycle dependent and independent patterns. Both the pressure and apical stress data reveals two phases of mechanical activity. A first, coupled to the wave of divisions in the mesoderm that transforms the embryo from the 64-to 76-cell stage, and a second wave within the 112-cell stage that appears decoupled from the cell cycle. The first phase of mechanical activity is initiated by the division of a mother cell, transiently retaining high pressures and apical/lateral/basal stresses (see SI) in the daughters, eventually dissipating just before the division of ectodermal cells. Both single-cell and tissue-based studies of cell divisions suggest a transient alteration in the local mechanical environment, and it is reassuring that our mechanical atlas makes manifest this anticipated feature.

The second phase of mechanical activity appears independent of any wave of cell divisions in the embryo, and is when the embryo begins to significantly alter its morphology. It is at this phase that clear anterior-posterior and animal-vegetal asymmetries in the mechanical atlas are manifest. Integrating all this data and distilling a mechanical logic is made challenging by the empirical observation that since the system is arguably close to a static equilibrium at all times the outward pushing cell pressures and contracting membrane and line elements must be in balance. Thus, the correlations apparent within the lineage maps between, say, cell pressures and apical stresses are both anticipated and confounding. Fortunately, our present theory of mechanics suggests a principled approach to integrating this complex data so as to reveal the independent mechanical features at a cell-by-cell basis.

The mechanical logic of ascidian gastrulation: While the above apical-centric view of the mechanical patterns demonstrates rich phenomena, it occludes a more holistic and 3D view of the gastrulating ascidian embryo. In particular, we require a means of integrating the, geometry-dependent, cellular pressures, surface tensions, and line tensions into cell-centric quantities. To achieve this, quantitatively, the physical concept of a mechanical stress tensor is useful. Inside cells, the cell volume is under an isotropic 3D stress σ σ σ iso = P α I I I (where I I I denotes the identical matrix), which is sourced by the pressure we measure. Over a membrane, surface tension cause an in-plane isotropic stress σ σ σ = -T αβ D D D (where D D D = t t t 1 ⊗ t t t 1 + t t t 2 ⊗ t t t 2 , t t t 1 and t t t 2 are tangential vectors of surface). Along each edge, line tension generate a 1D stress σ σ σ = -F αβγ E E E (where E E E = t t t ⊗ t t t and t t t is the tangential direction). Accounting for all of these, the average stress over a region of volume V is

σ σ σ =   α P α I I IdV - αβ T αβ D D DdA - αβγ F αβγ E E EdL   /V. ( 4
)
For convenience, we pick the animal-vegetal (AV), left-right (LR) and anterior-posterior (AP) axes as our orthonormal 3D coordinate system. As a result, we integrate all the data within a cell, and on its membrane and junctional surfaces and boundaries into a single 3 × 3 cell-based tensor. This tensor is symmetric by construction, and the scaling by the volume of the cell ensures that the tensor itself is an intensive quantity, whose magnitude is independent of the size of the system. Perhaps more importantly, constructing a principled quantity such as the cell-based stress tensor allows for the correct conglomeration of the, geometry-dependent, membrane and junctional tensions that can produce shear stresses on the cell. A simple way to see these patterns is to decompose the stress tensor, per cell, into its hydrostatic and deviatoric (its traceless part, which causes shape change) components. The second characteristic of the deviatoric matrix, J 2 , is related to a scalar quantity termed the von Mises Stress, σ V M = √ 3J 2 , which is a cell-based scalar measure of the extent of shear stresses in the system. We heatmap the variation of σ V M on the lineage diagram of the embryo in Figure 7A, and present its spatial patterns from an animal and vegetal view at the 112-cell stage stage in Figure 7B. The pattern of the von Mises stress in Figure 7A are manifestly distinct from either of the patterns seen in Figure 6. This is naturally expected since the pressure (hydrostatic) and deviatoric contributions to overall stress are independent. Contrasting previously observed patterns, ectodermal cells experience shear stress. To gain a more physical sense of the pattern we visualize the von Mises stress on the embryo in Figure 7B. One can clearly observe patterns, including the high levels of von Mises stress in the posterior mesodermal lineage.

Discussion

Recap of main findings

Increasing the precision, resolution, and robustness of our measurements of mechanical stresses in embryos is necessary to understand of the principles of morphogenesis and pattern formation. Here, we have developed a three-dimensional mathematical theory for static mechanical equilibrium in multicellular aggregates. The theory facilitates the solution of an inverse mechanics problem, and a data-processing and analysis scheme. Applying our 14/32 scheme to gastrulation in the ascidian embryo, we construct the first mechanical atlas at single-cell resolution in a whole intact embryo through the time course of a large scale change in shape. The quality of the model's fit to the data provides strong evidence that the cellular, membrane, and line forces are close to a static equilibrium with each other throughout the course of the dynamics. As such, a central result of the paper is that the mechanical dynamics are adiabatic, always maintaining close proximity to the manifold of static geometries, suggesting the presence of yet to be discovered mechanisms that ensure this. We find that cell pressures are patterned in space and time, appearing to correlate strongly with cell lineage and germ layers. We finally demonstrate our ability to measure spatio-temporal patterns of complex mechanical objects such as the von Mises stress that provides insight into the shear stresses in the system.

What we observe agrees with many of the known dynamical trends of the actomyosin cytoskeleton within this context [4]. A striking qualitative deviation is the possible role that spatio-temporal heterogeneities in cell pressures could be playing in the observed shape change. Thus far, in the context of ascidian gastrulation and most case studies of multicellular animal morphogenesis, the model is solely focused on the dynamics of the actomyosin cytoskeleton that imbue membranes and junctions with contractile forces. Here, running counter to a this very prevalent paradigm, we demonstrate that the observed shapes are concomitant with variations in cell pressures across the embryo.

New avenues of research

Our measurements of mechanical stress open up many exciting avenues for research. First, despite the statistical evidence for the accuracy of our measurements, falsification with alternate measurements of mechanical stresses must be pursued. Second, our computational approach to a mechanical atlas cannot identify the origins of the pressure heterogeneities we observe. While a 3D theory constructs an intensive conjugate variable to volume, unlike analogous measurements in 2D schemes, we cannot assess whether the pressure variations are generated by changes in its hydrostatic components or its more active components. Addressing what the mechanistic origins of the dynamical changes in cell pressures must be addressed in the future. Third, our formulation of a theory for the static equilibrium 15/32 (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission.

of 3D cellular aggregates is necessarily an approximation of a more dynamical phenomena. While velocity-dependent components to stresses have been ignored, it points to the importance of incorporating such effects in accurate and robust numerical schemes.

A longer term scientific goal that is made possible with the advance laid out in this study is the possibility of investigation the regulatory origins of mechanical stresses in embryos. In particular, armed with a mechanical atlas and a RNA-based cell atlas at single-cell resolution over time opens the possibility of identifying the regulatory changes that accompany observed mechanical trends in the data. Manifestly, the mechanical degrees of freedom are far fewer than the possible regulatory states thus making it possible to study the statistical properties of a highly degenerate genotype-to-phenotype map.

16/32

This derives

(f f f αβ • l l l βγδ )(f f f αγ • l l l βγ ) -(f f f αβ • l l l βγ )(f f f αγ • l l l βγδ ) = 0, (S5)
and

f f f αβ • l l l βγδ f f f αβ • l l l βγ f f f αγ • l l l βγ f f f αγ • l l l βγδ = 1. (S6)
There are e such constraints in total. However, when we multiply all these equations together,

<β,γ> f f f αβ • l l l βγδ f f f αβ • l l l βγ f f f αγ • l l l βγ f f f αγ • l l l βγδ = β γ f f f αβ • l l l βγδ f f f αβ • l l l βγ = 1 (S7)
The product in the brackets forms a loop, giving us

β (1) = 1 (S8)
an identity. So these e constraints have 1 degree of dependence. On the other hand, when leg orientations are constrained, the v leg lengths are free. In summary, the single dependence of orientation constraints is recounted, and there are thus C recounted constraints in total.

Figure S1 B. Dimensionality of a SCP tessellation

The geometry of a SCP tessellations can be defined by V vertex locations {r r r i } and F membrane centroids {ρ ρ ρ αβ }, which are not independent parameters and thus are not a useful parameterization of the system.

To achieve one let us first consider the constraints on 3F centroid parameters. Each edge, a circular arc, is an intersection of three spheres. So the central axis of the circle goes through all three corresponding centroids, as shown in Figure S1b. Thus, the three centroids are co-linear. For ρ ρ ρ αβ , ρ ρ ρ βγ and ρ ρ ρ γα , the co-linear relation is given by

ρ ρ ρ αβ -ρ ρ ρ βγ = χ αβγ (ρ ρ ρ αβ -ρ ρ ρ γα ), ( S9
)
where χ αβγ is ratio of lengths connecting the three centroid, and the index sequence represents the sequence of centroids in equation. This equation gives 2 constraints since the ratio is free. As such, there are 2E constraints in total, but manifestly they are coupled. At each vertex r r r αβγδ , the six centroids are coplanar, as shown in Figure S1c. Let us suppose the first three co-linear constraints, corresponding to χ βγα , χ γδα and χ δβα , are satisfied. These three axes determine a plane, so ρ ρ ρ βγ , ρ ρ ρ γδ and ρ ρ ρ δβ can only move in plane. This suggests that there ought to be one constraint instead of two for the last co-linear constraint. The last constraint can be given by

χ βγα χ γδα χ δβα = 1, ( S10
)
18/32 according to Menelaus' theorem. There is one dependence of constraints at each vertex, and thus we recounted V constraints in total. We also observe that for each cell, if we construct a product of all the Menelaus equations of its vertices, r r r αβγδ ∈α χ βγα χ γδα χ δβα = e e e αβγ ∈α χ βγα χ γβα = 1, (S11) it gives us an identity. That is because χ βγα χ γβα = 1, and thus all left-hand-side terms cancel each other. In summary, 2E constraints of centroids have V + C dependence in SCP tessellation, so the centroid parameters have

3F -(2E -V -C) = 3F -3V + C = 4C (S12)
degrees of freedom.

We then consider the constraints on 3V vertex parameters, once centroids are given. Vertices of one membrane are equidistant from its centroid,

r r r αβ,1 -ρ ρ ρ αβ = r r r αβ,2 -ρ ρ ρ αβ = ... = r r r αβ,v αβ -ρ ρ ρ αβ , ( S13
)
which gives v αβ -1 constraints. Thus there are 6V -F constraints in total. However, these constraints are not independent. Consider an edge e αβγ and its two end vertices, r r r i and r r r j . When |r r r

i -ρ ρ ρ αβ | = |r r r j -ρ ρ ρ αβ | and |r r r i -ρ ρ ρ βγ | = |r r r j -ρ ρ ρ βγ | are satisfied, |r r r i -ρ ρ ρ γα | = |r r r j -ρ ρ ρ γα |
is satisfied automatically, since the three centroids are co-linear. This is shown in Figure S1b. Therefore, we have 6V -F -E independent constraints on vertices in SCP tessellation, which means

3V -(6V -F -E) = F -V = C (S14)
degrees of freedom of vertex parameters. All in all, SCP tessellation has 5C degrees of freedom.

C. Weighted Voronoi tessellation

A Voronoi tessellation of C cells in 3D space is defined by C sites, and therefore has 3C degrees of freedom. {q q q α } α=1,2,...,C . A cell region R α is a set of points which are closer to q q q α than to other sites. R α = {r r r|d α (r r r) < d i (r r r), ∀i = α}, and here d α = |r r r -q q q α | is the Euclidean distance. So the boundary B αβ between any two neighboring cells is the perpendicular biosector of two corresponding sites, αβ = {r r r|d α (r r r) = d β (r r r)}, which is a flat plane. So the shape of a cell region is a polyhedron. Further more, the edge E αβγ is the intersection of three bisectors and has equal distance to three neighbor sites, which satisfies E αβγ = {r r r|d α (r r r) = d β (r r r) = d γ (r r r)}. The vertex r r r = r r r αβγδ satisfies d α (r r r) = d β (r r r) = d γ (r r r) = d δ (r r r). Now, we can modify the definition of distance by subtracting (or adding) a weight θ 2 α to the squared Euclidean distance, d 2 α (r r r) = |r r r -q q q α | 2 -θ 2 α . In this situation, any two points r r r 1 and r r r 2 at the boundary B αβ satisfy |r r r 1 -q q q α | 2 -|r r r 1 -q q q β | 2 = θ 2 α -θ 2 β = |r r r 2 -q q q α | 2 -|r r r 2 -q q q β | 2 . (S15)

It can be simplified as r r r 1 • (q q q α -q q q β ) = r r r 2 • (q q q α -q q q β ). ( S16) or (r r r 1 -r r r 2 ) • (q q q α -q q q β ) = 0 (S17) r r r 1 -r r r 2 is thus perpendicular to q q q α -q q q β , and the boundary face is still a perpendicular plane to the line between two sites. The shape of cell region is still a polyhedron. This weighted Voronoi tessellation has 4C degrees of freedom (the three coordinates of q q q α and the weights θ α , same as the dimensionality of polyhedron tessellations. Therefore, any polyhedron tessellation can be described by a weighted Voronoi tessellation. Said another way, there exists a one-on-one mapping between these two families of tesselations. Now, we can further modify the definition of distance by multiplying a power p α to the definition of distance, d 2 α (r r r) = p α (|r r r -q q q α | 2 -θ 2 α ). In this situation, the point r r r at boundary B αβ satisfies p α (|r r r -q q q α | 2 -θ 2 α ) = p β (|r r rq q q β | 2 -θ 2 β ). Multiply this equation by (p α -p β ) and simplify, there is p 2 α (r r r -q q q α ) 2 + p 2 β (r r r -q q q β ) 2 -2p α p β (r r r -q q q α ) • (r r r -q q q β ) = p α p β (q q q α -q q q β ) 2 + (p α -p β )(p α θ 2 α -p β θ 2 β ). (S18)

19/32

We realize that the left hand side is a square and the right hand side is a constant independent of r r r. So this equation can be further simplified as (p α -p β ) 2 (r r r -p α q q q α -p β q q q β p α -p β ) 2 = const. (S19)

Let us define ρ ρ ρ αβ = p α q q q α -p β q q q β p α -p β (S20)

and

R 2 αβ = p α p β (q q q α -q q q β ) 2 + (p α -p β )(p α θ 2 α -p β θ 2 β ) (p α -p β ) 2 . (S21)
Then the equation of r r r is |r r r -ρ ρ ρ αβ | 2 = R αβ , which is a sphere with centroid ρ ρ ρ αβ and radius R αβ . So all boundaries are spherical sections. Therefore, this generalized weighted Voronoi tessellation is a SCP tessellation and has 5C degrees of freedom. This indicates that there is a one-on-one mapping between SCP tessellation and generalized weighted Voronoi. Any SCP tessellation can be parameterized by generalized weighted Voronoi.

Part II: Mechanical Dual Graph

A. Divergence theorem of vector area

In this part, we'll solve the force balance equations ( 1a)-( 1c) by introducing the mechanical dual graph. Every balance relation corresponds to a dual shape, and every force corresponds to the vector area of a face in this dual shape. The concept of a "vector area" is a vector S S S = S n n n with the face area S as the magnitude and the face normal n n n as the direction. For a curved face, the vector area is the integral of face elements S S S = n n ndS. first claim that the sum of all face vector areas of a compact closed shape is zero:

S S S i = 0 0 0. ( S22
)
This result is a simple consequence of the Divergence theorem. For a vector or tensor field F F F , the divergence theorem on a compact closed volume is given by

S (F F F • n n n)dS = V (∇ • F F F )dV. (S23)
Taking the identity tensor field F F F = I I I, we get

S n n ndS = 0 0 0. (S24)
So the integral of surface vector area for a compact shape is zero.

B. Local dual graph for 2D cellular lattice

In a 2d cellular lattice, the equilibrium geometry is a Circular Arced Polygon (CAP) tessellation [2]. This can be parameterized by 2d generalized Voronoi tessellation, taking Ψ = {q q q α , θ α , p α : α = 1, 2, ..., C} as parameters.

For any point r r r on the membrane of cell α, we define the local dual point as q q q * α (r r r) = r r r + p α (q q q α -r r r).

(S25)

The two local dual points form the local dual line q * αβ corresponding to the membrane between cell α and β, with the length:

q * αβ = q q q * α -q q q * β = (p α -p β )(r r r -ρ ρ ρ αβ ) = |p α -p β |R αβ . (S26)
This means the length is a constant along the membrane. Also the local dual line is parallel to the radius and thus perpendicular to membrane at r r r. At a vertex, the three local dual lines form a local dual triangle. Along the membrane, the local dual line is swept to form a curved quadrilateral. See Figure S2b.

Figure S2

We now consider the local dual triangle q q q * α q q q * β q q q * γ constructed from vertex r r r i . Applying the 2d divergence theorem to this local dual triangle gives us

q * αβ n n n 1 + q * βγ n n n 2 + q * γα n n n 3 = 0. (S27)
We notice that the three vectors have same directions as the corresponding membrane tension balance at r r r i :

T αβ n n n 1 + T βγ n n n 2 + T γα n n n 3 = 0. (S28)
Therefore, the length of local dual line is proportional to the membrane tension:

T αβ = cq * αβ , ( S29
)
where c is a constant for all membranes, representing the scale factor. We consider next the edge element dr r r of membrane αβ, the dual graph is a curved quadrilateral as shown in Figure S2c. Applying the 2d divergence theorem to the quadrilateral gives us

q * αβ n n n 1 + q * αβ n n n 2 + dq * α n n n 3 + dq * β n n n 4 = 0. (S30)
We also notice the vectors have same directions as the corresponding tension-pressure balance on dr:

T αβ n n n 1 + T αβ n n n 2 -P α drn n n 3 -P β dr n n n 4 = 0, (S31)
which is the element-wise state of the Young-Laplace relation. Since we have equation ( S29) and n3 = -n 4 , the pressure difference would be:

P α -P β = -c dq * α -dq * β dr = c(p α -p β ). (S32)
Therefore,

P α = cp α + b, ( S33
)
where c is the same constant as in equation (S29), and b is another constant representing the background pressure.

C. Local dual graph for 3D cellular aggregate

As in 2d, we define the local dual point of r r r on the membrane of cell α as equation (S25). The local dual line of membrane αβ has constant length q * α given by equation (S26), and it is also perpendicular to the membrane at r r r. For a point r r r at the edge αβγ, the three local dual line form a local dual triangle. The triangle area along the edge is a constant: S * αβγ = 1 2 q q q * αβ × q q q * βγ , (S34)

where q q q * αβ = q q q * α -q q q * β . Additionally, the local dual triangle is orthogonal to the tangential direction of the edge at r r r. These properties allow us to infer forces from local dual networks.

At a vertex r r r αβγδ where four cells meet, the four local dual points q q q * α , q q q * β , q q q * γ and q q q * δ form a local dual tetrahedron, with four local dual triangle faces. According to the Divergence theorem, the four vector areas of triangles satisfy where n n n 1,2,3,4 are corresponding normal vectors. Consider the line tension balance at this vertex:

F αβγ n n n 1 + F αβδ n n n 2 + F αγδ n n n 3 + F βγδ n n n 4 = 0 0 0, (S36)
where n n n 1,2,3,4 are tangential vectors of edges. According to the local dual properties, these vectors match up to the normal vectors in equation (S35). the four unknown line tensions satisfy three linear equations with a specific solution {S * αβγ , S * αβδ , S * αγδ , S * βγδ }. As such, the value of line tension is proportional to the triangle area:

F αβγ = cS * αβγ . (S37)
Since the triangle area is constant along the edge, the coefficient c here should be a global constant so that the solution of tensions matches with the equations of the other vertices.

First, as we mentioned previously, the line tension is proportional to S * αβγ , ensuring that the line tensions balance at a vertex. Second, the local dual triangle sweep out a triangular tube along the edge, where the three lateral faces are orthogonal to membrane tensions, as shown in Figure S3e. According to the Divergence theorem, for a small tube element corresponding to dr r r, the five area vectors satisfy

S S S

(1) αβγ + S S S

(2) αβγ + dS S S αβ + dS S S βγ + dS S S γα = 0.

(S38)

This corresponds to the force balance on an edge (balance of 2 line tensions and 3 membrane tensions). Since the triangle area is proportional to line tensions and is orthogonal to line tension, cS S S

(1) αβγ + cS S S

(2) αβγ = F αβγ κ κ κdr. Therefore, the membrane tension is proportional to the lateral face area cdS αβ = T αβ dr. Hence, the membrane tension is

T αβ = c dS αβ dr = cq * αβ 1 2 dq * α dr + dq * α dr = cq * αβ 1 2 (2 -p α -p β ). (S39)
Third, for a small surface element dr 1 dr 2 on a membrane, the dual line sweeps out a quadrilateral prism, as shown in Figure S3f. In the same manner, we have

P α = cS * α dr 1 dr 2 = c dq * α dr 1 dq * α dr 2 = c(1 -p α ) 2 . (S41)
We have now solved for all the forces using the local dual network. However, the local dual network is not unique. If we set q q q * α (r r r, λ) = r r r -λp α (r r r -q q q α ), which implies a rescaling of the local dual by λ. In this rescaled dual network, all the orthogonality remains and thus we can use it to infer force. The solutions would be

F αβγ = cλ 2 S * αβγ , T αβ = 1 2 c(2 -λ(p α + p β ))λq * αβ , P α = c(1 -λp α ) 2 . (S42)
Relabel the coefficients by c = cλ 2 and b = 2cλ. Besides this, we can add another constant a to pressure that will leave the balances unaffected. Then the solutions of force are

F αβγ = cS * αβγ , T αβ = 1 2 (b -c(p α + p β ))q * αβ , P α = cp 2 α -bp α + a.
(S43)

D. Mechanical zero mode

Let us take a look at a simple two-cell example as a case study for the application of force inference and to facilitate a physical understanding of the three zero modes in the solution.

To set the scene, the two-cell system is defined by the following geometry parameters as: q q q 1 = [1, 0, 0], q q q 2 = [-1, 0, 0], p 1 = p 2 = 1, θ 1 = θ 2 = 2. And also the background parameters: q q q 0 = [0, 0, 0], p 0 = 0, θ 0 = 0. This results in two spherical membranes and a flat membrane in the middle (Figure S4 a andb). The centroids and radius are: ρ ρ ρ 10 = q q q 1 , ρ ρ ρ 20 = q q q 2 , R 10 = R 20 = 2; ρ ρ ρ 12 = ∞, R 12 = ∞. Now the forces can be inferred from these parameters. Take any point r r r on the edge (e.g. r r r = [0, 0, √ 3]), the local dual points are: q q q * 1 = q q q 1 , q q q * 2 = q q q 2 , q q q * 0 = r r r. Then the local dual line and triangle are: q * 10 = q * 20 = q * 12 = 2, S * 123 = √ 3. So the force solution is:

P 0 = a, P 1 = P 2 = a -b + c, T 12 = b -2c, T 10 = T 20 = b -c, F 012 = √ 3c. (S44)
Here the constant a represents the background pressures in the simple system considered here. Thus, increasing the pressures everywhere by a constant value leaves the pressure differences unchanged and the system maintains its balance. The magnitude of constants b and c represent the unit of all forces. Multiplying them by same scalar factor will not affect the relative force configurations. But the ratio µ = c/b does. In this example, take c = 0 (non-line-tension case), then the three membrane tension are equal T 10 = T 20 = T 12 = b. When considering a small amount of line tension (e.g c = 0.1), we see the middle tension is lower than the other two membranes T 10 = T 20 > T 12 so that the forces are balanced at edge. Additionally, we assume that the cytoskeletal machinery at membranes and lines is contractile, which suggests that solutions of F and T must be positive, which requires µ > 0 and µ < max αβ (p α + p β ) -1 .

Figure S4 E. 2D reduction of 3D mechanical constraints

The existence of mechanical dual graph results in the solvability of force balance equations. The consequence is that every polyhedral or SCP tessellation is a mechanical equilibrium geometry. However, in 2D, as noted by Noll et al. [3,2], all polygonal tessellations are not in equilibrium. In particular, a dual network only exists when the outer angles of each cell satisfy

i∈V α sin φ iγ sin φ iβ = 1, ( S45
)
where φ iγ and φ iβ are shown in Figure S4c. What underlies this distinction? In particular, why does any polyhedral tessellations correspond to a state at mechanical equilibrium, while polygonal tessellations do not?

Here we demonstrate that the 3D constraints reduce smoothly to the 2D constraints, and that the apparent additional mechanical constraints in 2D are geometric in origin. We can investigate membrane F 0α and its adjacent faces in 3D, shown in Figure S4c. Coplanarity of face F 0β gives the equation (r r r i -r r r j ) • (e e e i ∧ e e e j ) = 0, (S46)

where e e e i and e e e j is the direction of edges from r r r i and r r r j . We use z z z to denote the normal vectors of F 0α , and use f f f β to denote the unit vector along F 0β and perpendicular to E 0αβ . There is a parallel relation that (z z z ∧ f f f β ) (r r r i -r r r j ), from which we can derive The first product is the same as the left hand side of equation, recapitulating the 2D constraint. Thus, the 2D constraints can be derived from 3D coplanarity constraints. In the 2D limit, when the angles between faces tend to zero, the coplanarity constraint is released, thus giving polygonal tessellations additional degrees of freedom.

(z z z ∧ f f f β ) • (e

Part III: Numerical Method for Geometric Reconstruction

A. Mean-squared-deviation of membrane pixels

In this part, we'll recover the geometric parameters Ψ = {q q q α , p α , θ α : α = 1, 2, ..., C} from a three-dimensional segmented membrane image. In the image, we know all the membrane pixel locations {r r r i : i = 1, 2, ..., N } and the two cells α i , β i adjacent to the membrane in question. The model assumes that a membrane pixel lies on a section of sphere, denoting the centroid ρ ρ ρ i = ρ ρ ρ α i β i and the radius R i = R α i β i . However, we anticipate a deviation of the pixel from the sphere i = |r r r i -ρ ρ ρ i | -R i in empirical data. Since ρ ρ ρ i and R i are both functions of Ψ by equation, we can recover the geometric parameters by least-square fitting to a SCP tessellation. Specifically, we minimize the mean-squared-deviation (MSD) function:

E(Ψ) = 1 N N i=1 2 i = 1 N N i=1 [|r r r i -ρ ρ ρ i (Ψ)| -R i (Ψ)] 2 . (S50)

B. Initial guess of least-squared fitting

In our model, the shape of membrane must be sphere (or plane -a sphere with 0 curvature). So we first use a sphere to fit each membrane, and record the best fit sphere's centroid and radius {ρ ρ ρ * αβ , R * αβ }. Manifestly, these parameters do not satisfy the constraints of equilibrium. We try to define {q q q 0 α }, {p 0 α } and {θ 0 α } successively. For an interior membrane F αβ , due to the expression ρ ρ ρ αβ = p α q q q α -p β q q q β p α -p β , three points q q q α , q q q β and ρ ρ ρ αβ should be co-linear. Therefore, we need to minimize the error E 1 (q q q α ) = mem in 1 -|(q q q α -q q q β ) • (q q q α -ρ ρ ρ * αβ )| |q q q α -q q q β ||q q q α -ρ ρ ρ * αβ | 2 (S51)

But for an outer membrane F 0α , such a constraint does not work since q q q α = ρ ρ ρ 0α . But we still can add terms on E 1 that keep q q q α close to ρ ρ ρ * 0α .

E 1 (q q q α ) = mem in 1 -|(q q q α -q q q β ) • (q q q α -ρ ρ ρ * αβ )| |q q q α -q q q β ||q q q α -ρ ρ ρ

* αβ | 2 W αβ + mem out |q q q α -ρ ρ ρ * 0α | R * 0α 2 W 0α (S52)
where W αβ and W 0α are weights for each membrane. Due to the form of main error function E (equation [WHICH EQUATION]), we notice that larger membrane should be weighted more than smaller membrane. Therefore, we set W αβ = (n αβ /n 0 ) 2 , where n 0 is the characteristic membrane size. For outer membranes, because the error function has a different form, we set W 0α = ω(n 0α /n 0 ) 2 . Here we set ω small since the geometry of an outer membrane may have a segmentation error larger than their interior counterparts (they are darker than interior membranes in an typical image).

After this minimization step (equation (S52)), we obtain an initial guess for dual points {q q q 0 α }. Now we can use sphere to fit each membrane with the constraint that ρ ρ ρ * * αβ is on the line defined by q q q 0 α and q q q 0 β . And we get a distance ratio k αβ = |q q q 0 α -ρ ρ ρ * * αβ | |q q q 0 β -ρ ρ ρ * * αβ | , which should be ratio of pressure k αβ = p β /p α if equilibrium. Therefore, we try to minimize the error

E 2 (p α ) = mem in (p β -p α k αβ ) 2 W αβ (S53)
where

W αβ = n αβ /n 0 |q α -q β | 2
is the weight for each membrane. Note that we don't consider outer membrane in this step because we have already set p 0 = 0. Also note that k αβ > 1, if not, we can replace p β -p α k αβ by p α -p β (1/k αβ ).

After this minimization step (equation S53), we get the initial guess of pressures {p 0 α }. Now we can have ρ ρ ρ * * * αβ = p 0 α q q q 0 α -p 0 β q q q 0 β p 0 α -p 0 β as the centroid for a membrane. These centroids can exactly satisfy the constraint in the model (co-linearity of three centroids of an edge). Using a sphere with a centroid on ρ ρ * * * αβ to fit each membrane, the radius R * * * αβ can be calculated. In order to make the radius satisfy constraints in the model we minimize the error

E 3 (θ α ) = mem |R * * * αβ -R αβ (p 0 α , q q q 0 α , θ α )| 2 W 2 αβ (S54)
and obtain an initial guess of {θ 0 α }.

Finally, we use these initial guesses {p 0 α , q q q 0 α , θ 0 α } to minimize the primary error function. This step can be done using MATLAB in built function fmincon or lsqnonlin. Perturbing the initial guess a little, and repeating the final minimization step, improves the result further. Repeated cycles of of perturbation-minimization steps can produce a strong convergence to a global minimum.

C. Synthetic verification of numerical scheme

In order to verify our numerical scheme for recovery of the empirically observed geometry , we generate in-silico 3d image data. We first initialize the geometric parameters Ψ 0 = q q q α , p α , θ α in a given 3d region. We cut the 3d region into cubic grids of a unit size. In each grid we randomly pick a site q q q α . The weights θ α and powers p α are sampled from a Gaussian distribution centered at the unit scale and 1, respectively. The background parameters are set to be zero.

Then we set the resolutions of the 3d image, i.e. set the pixel size to be 1/100 of the unit scale, which is close to the resolution of real imaging data. Now we are able to generate Voronoi tessellation of this 3d image. For each pixel, we measure the distance of weighted Voronoi definition to the sites and label the pixel by the cell index with smallest distance. After labeling every pixel, we mark the pixels at cell boundary as membrane pixels. We then add Gaussian noise on the membrane pixels locations and apply the MSD minimization to infer the best fitting parameters Ψ.

Given different level of noise, from 0% to 20%, normalized by the unit scale, we find out that the inferred Ψ is 99% close to the ground truth Ψ 0 . This indicates the MSD minimization is a robust way to recover the geometry from noisy imaging data. Taking a close look at the parameters, we observe that the sites q q q α and the weight θ α is more robust than the powers p α .

D. Sensitivity analysis

At the minimum of the MSD function (equation 5), the deviation of E(Ψ, r r r) ought to be zero: where n n n i = (r r r i -ρ ρ ρ i )/|r r r i -ρ ρ ρ i | is the normal direction of pixel to the sphere. Therefore, the deviation at r r r i is

δE(Ψ, r r r) = 1 N N i=1 2 [|r r r i -ρ ρ ρ i | -R i ] δ [|r r r i -ρ ρ ρ i | -R i ] = 0, ( S55
n n n i • δr r r i -(n n n i • ∂ρ ρ ρ i ∂Ψ + ∂R i ∂Ψ )δΨ = 0. (S57)
Writing all these equations together, we get a linear approximation at the minimum Kδr r r + M δΨ = 0. The vector δr r r = [δr 1x , δr 1y , δr 1z , ..., δr ix , δr iy , δr iz , ...] T is the deviation of pixel locations. And the elements of matrix K are M i,5α-3 = -n ix ∂ρ ix q αx -n iy ∂ρ iy q αx -n iz ∂ρ iz q αx -∂R i q αx M i,5α-3 = -n ix ∂ρ ix q αy -n iy ∂ρ iy q αy -n iz ∂ρ iz q αy -∂R i q αy M i,5α-3 = -n ix ∂ρ ix q αz -n iy ∂ρ iy q αz -n iz ∂ρ iz q αz -∂R i q αz M i,5α = -∂R i ∂θ α (S59)

[K i,3i-2 , K i,3i-1 , K i,3i ] = [n ix ,
where i = 1, 2, ...N ; α = 1, 2, ..., C.

Figure 1 .

 1 Figure 1. The mechanical ingredients of the theory. A) Surfaces of contacts between cells (labeled α and β are defined by two numbers corresponding to their mean curvature, H αβ , and a homogeneous and isotropic effective surface tension, T αβ . The Young-Laplace Law relates the difference between the pressures in the two adjacent cells to the surface tension and mean curvature of the surface of contact, |P α -P β | = 2T αβ H αβ . B) Lines of contact are formed when three surfaces of contact coincide, and admit their own effective line tension. The 1D analogy of the usual Young-Laplace Law relates the three surface tensions, which attempt to curve the line, to its curvature and line tension through T T T αβ + T T T βγ + T T T γα = -F αβγ κ κ κ. C) Multiple lines of contact meet at vertices within multicellular aggregates. The static theory developed assumes that the line tension balance each other at vertices F F F αβγ + F F F βγδ + F F F αβδ + F F F αγδ = 0.

Figure 2 .

 2 Figure2. Increasingly complex voronoi diagrams. A) 2D voronoi diagram, the blue dots are voronoi sites q q q, the black lattice are cell regions {r r r|d 2 α (r r r) ≤ d 2 j (r r r), ∀j = α}, under Euclidean distance d 2 α (r r r) = |r r r -q q q α | 2 . B) 2D weighted voronoi diagram, with same sites as A), cell regions under the weighted Euclidean distance d 2 α (r r r) = |r r r -q q q α | 2 -θ 2 α . Size of blue dots represents the weight θ. C) 2D generalized weighted voronoi tessellation, with same sites and weights as B), cell regions under the distance d 2 α (r r r) = p α (|r r r -q q q α | 2 -θ 2 α ). The color scale of the blue dot represents the power p.

Figure 3 .

 3 Figure 3. Analytical solution to the inverse mechanical problem. (A) Divergence theorem of vector area. The dark gray square represents the vector area, a face element dS with the normal vector nn n. The integration of vector are over the surface of a closed 3D object is zero. (B) Local dual graph of 2D GWV tessellation (black line). The blue dots are voronoi sites q q q, and the red dots are local dual points q q q * with respect to points r r r on membranes or vertices. The dark red triangles are dual to the vertices, the light red curved quadrilaterals are dual to edges. (C) The divergence theorem of local dual graph reveals mechanical balance in 2D. Membrane tension balance T T T αβ + T T T βγ + T T T γα = 0 corresponds to triangle divergence theorem q * αβ n n n 1 + q * βγ n n n 2 + q * γα n n n 3 = 0 at vertex r r r i . Pressure-tension balance T αβ n n n 1 + T αβ n n n 2 -P α dr n n n 3 -P β drn n n 4 = 0 corresponds to quadrilateral divergence theorem q * αβ n n n 1 + q * αβ n n n 2 + dq * α n n n 3 + dq * β n n n 4 = 0 of edge element dr. (D) A two-cell example of the mechanical zero mode. In the upper case, no line tension results in three equal membrane tensions. In the lower case, line tension results in smaller middle membrane tension than two apical tensions.

Figure 4 .

 4 Figure 4. Geometry reconstruction from image. (A)Synthetic image of cellular aggregate. (B)Accuracy of fitting scheme on synthetic data with noise. (C)Distribution of eigenvalues at the minimum of error function on a synthetic data. (D)Error of pixel locations from a sphere in SCP tessellation. (E)Error distribution in an image of a 64-cell stage Phallusia mammillata embryo. (F)Spatial distribution of error map for the same image as E).

Figure 5 .

 5 Figure 5. The mechanical atlas for gastrulation process in Phallusia mammillata embryo. The top three panels show the inferred pressures, apical surface tensions and apical line tensions in a Phallusia mammillata embryo at the 64-cell, 76-cell and 112-cell stages, in vegetal view. The lower panels show the pressures and membrane tensions in the frontal cross section of the middle line in vegetal view. The undetermined zero-mode parameter µ is set to be µ = µ 0 /2. The color maps give the relative values of inferred forces.

Figure 7 .

 7 Figure 7. Lineage map and spatial visualization of von Mises stress for gastrulation process in Phallusia mammillata embryo. A) The lineage maps of von Mises stress for all cells from 64-cell stage through 112-cell stage. B) The vegetal, sagittal and animal view of von Mises stress at 112-cell stage.

  e e i ∧ e e e j ) = 0. (S47) This equation can be written asf f f β • e e e i f f f β • e e e j = z z z • e e e i z z z • e e e j , (S48)where the left hand side is exactlysin φ iβ sin φ jβ. Taking a product around a given face, F 0α , gives β

)

  So for anyi, there is δ [|r r r i -ρ ρ ρ i (Ψ)| -R i (Ψ)] = 0. In this equation, δ|r r r i -ρ ρ ρ i (Ψ)| = n n n i • (δr r r i -δρ ρ ρ i ) = n n n i • (δr r r i -

7/32 pressure

  balance on dr, T αβ n n n 1 + T αβ n n n 2 -P α drn n n 3 -P β dr n n n 4 = 0, which is the differential statement of the Young-Laplace equation. Since T αβ = cq * αβ and n3 = -n 4 , the pressure difference can be rewritten as P α

  four lateral area vectors add up to 2T αβ H αβ dr 1 dr 2 . So the upper area is cS * α = P α dr 1 dr 2 . The pressure is

	S S S * α + S S S * β + S S S	(1) αβ + S S S	(2) αβ + S S S	(3) αβ + S S S	(4) αβ = 0.					(S40)
	This corresponds to a force balance on a membrane (balance of 4 membrane tensions and two pressures). As proved
	previously, cS αβ = cS (1) αβ = T αβ dr 1 , and thus cS S S (3)	(1) αβ + cS S S	(3) αβ = T αβ dr 1	dr 2 R αβ n	. Similarly, cS S S	(2) αβ + cS S S	(4) αβ = T αβ dr 2	dr 1 R αβ n	.
	Therefore, the

  n iy , n iz ], i = 1, 2, ..., N (S58)and K i,j = 0 otherwise. Similarly, δΨ = [δp 1 , δq 1x , δq 1y , δq 1z , δθ 1 , ..., δp α , δq αx , δq αy , δq αz , δθ α , ...] T is the deviation of geometric parameters. Thus the elements of the matrix M is

	M i,5α-4 = -n ix	∂ρ ix p α	-n iy	∂ρ iy p α	-n iz	∂ρ iz p α	-	∂R i p α

8/32(which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission.

Supplementary information Part 0: Formulation of Force Inference Problem

We assume a highly simplified model for the dominant mechanical contributions to a multicellular aggregate: cellular pressures P α , membrane tensions T αβ of each membrane, and the line tension F αβγ of each junctional edge. Indices label cells, where membranes and edges are indexed by the two and three cells that define them, respectively. The mechanical balance equations between them are:

T T T αβ + T T T βγ + T T T γα = -F αβγ κ κ κ, (S1b)

Here H αβ is the mean curvature of surface and κ κ κ is the curvature vector of 3d curve.

The force inference problem is to infer the values of P α , T αβ and F αβγ from a given 3d segmented image of membranes. In order to solve equation ( S1a)-(S1c) (Part II), the geometric parameters in equations must be extract from an image (Part III). This requires a parameterization of the geometry of cellular aggregate.

Part I: Parameterization of Cellular Aggregate Geometry

By assumption, the pressures and tensions are homogeneous and isotropic over each cell, membrane and edge. So the membrane shapes are constrained to be spheres (or planes -spheres with no curvature) and the edges are circular arcs (or straight lines -circular arcs with no curvature). Therefore, the equilibrium geometry of a cellular aggregate is polyhedron tessellation or spherically curved polyhedron (SCP) tessellation. Before parameterizing the two tessellations, we count their respective dimensionalities.

A. Dimensionality of a polyhedron tessellation

We use C, F , E, V to denote the number of cells, membrane faces, edges and vertices, respectively. According to Euler's formula, for each polyhedron we have f + v -e = 2, where f , v and e denoting the numbers for one cell. When we sum the C Euler relation across all cells each face counts twice, each edge counts 3 times, each vertex counts 4 times, giving 2F + 4V -3E = 2C. An additional relation is conferred by the fact that each edge has two vertices and each vertex has four edges, giving 2E = 4V . Using these two relations we derive

The geometry of any polyhedral tessellation can be defined by the positions of its vertices, {r r r i }, which are 3V numbers. However, these parameters are not independent owing to the constraint that the faces must flat. We use v k to denote the vertex number of kth membrane. Then the vertices on this membrane are under v k -3 constraints to enforce coplanarity, since three points can define a plane. Thus there are 6V -3F constraints in total, since each vertex is counted 6 times. We will show in next paragraph that these constraints are not independent. We recounted 1 constraint for each cell. So in all there should be 6V -3F -C independent constraints. Therefore, a polyhedron tessellation has

Where does the recounted constraint come from? For one polyhedron cell α, after the coplanarity constraints of its faces are satisfied, let us consider the v outer vertices that connect to this cell. As shown in Figure S1a, each vertex connects to one polyhedron vertex. They have 3v degrees, including v lengths of corresponding 'legs' and 2v leg orientations, {l l l βγδ }. As we counted before, there ought to be e constraints. Let us first consider the constraints on leg orientations. The coplanar object comprising one edge and two legs at its ends is given by (l l l βγδ ∧ l l l βγ ) • e e e αβγ = 0, (S3)

where e e e αβγ denotes the orientation of edge. If we use f f f αβ to denote the normal vector of face, then the following relation holds true, e e e αβγ = f f f αβ ∧ f f f αγ . Allowing us to write the coplanarity constraint as

17/32

(which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission.

Therefore, the relation between observed pixel's perturbation δr r r and deviation in geometric parameters δΨ is given by δΨ = Lδr r r = (-M -1 K)δr r r, where M -1 is the pseudoinverse matrix of M . As long as the error respond matrix L has eigenvalues smaller than 1, the inverse method is determined to be robust. For synthetic data, we compute the eigenvalues λ i of L and demonstrate that the magnitudes are all smaller than 1.

E. Some results on Ascidian gastrulation data

Figure S5

After verifying the geometry reconstruction method on synthetic data, we apply it onto ascidian gastrulation images. The images are well segmented and labeled. We focus on the stage from 64-cell to 112-cell. During the non-dividing phase, the fitting error stays in low value, and uniformly distributed in space. During the dividing phase, the cells are like dumbbell shapes before the new membrane forms (Figure S5a). This causes the high fitting error on such cells, since we use spherical shape to fit the membrane. Figure S5b shows the error distribution of embryos during the dividing phase from 64-cell to 76-cell stage. We realize that the fitting error is still in low value for other cells, thus we could still trust the force inference results on those cells. The lineage map of errors (Figure S5c) has a clearer presence. It also shows that the errors of daughter cells fall back to low values immediately after dividing.

Then we infer the forces of all cells. As we illustrate in previous part, there is one undefined parameter µ that would affect the configurations of force results. Additionally, this parameter is bounded by 0 and µ 0 to ensure the tensions are contractile. For this embryo data, the typical value of µ 0 is 0.1, which means the effect of the zero mode is up to 10% to the force configurations. While we have shown the embryo force maps with µ = µ 0 /2, i.e. µ = 0.05 in the main text, Figure S6 and Figure S7 shows the force maps with µ = 0 and µ = µ 0 .

When consider the force effect in a course-grained view, the contractile stress applying on cell is under even smaller variations. For instance, given higher value of µ, the line tensions provide higher contractile stress; meanwhile, lower value of the membrane tensions provide lower contractile stress, canceling out the line tension effect to a certain extend. Using dimensional analysis, suppose the length scale is L, the membrane tension scale is q, then the line tension scale is q 2 . Therefore, the line tensions provide stress of scale q 2 L, thus the zero mode effect is µq 2 L; similarly, the effect stress by membrane tension is of scale -µpqL 2 . By the analysis of mechanical dual graph, there is pL ∼ q. Therefore, the two effect stresses have the same scale and opposite sign, resulting in the cancelation of zero mode effects. These invariant patterns of stress are shown in lineage maps (Figure S8, Figure S9 and Figure S10).

27/32

(which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission.