Variational Method for Image-Based Inference of Internal Stress in Epithelial Tissues

Cellular mechanics drives epithelial morphogenesis, the process wherein cells collectively rearrange to produce tissue-scale deformations that determine organismal shape. However, quantitative understanding of tissue mechanics is impaired by the difficulty of direct measurement of stress in vivo. This difficulty has spurred the development of image-based inference algorithms that estimate stress from snapshots of epithelial geometry. Such methods are challenged by sensitivity to measurement error and thus require accurate geometric segmentation for practical use. We overcome this difficulty by introducing a novel approach—the variational method of stress inference (VMSI)—which exploits the fundamental duality between stress and geometry at equilibrium of discrete mechanical networks that model confluent cellular layers. We approximate the apical geometry of an epithelial tissue by a 2D tiling with circular arc polygons in which arcs represent intercellular interfaces defined by the balance of local line tension and pressure differentials between adjacent cells. The mechanical equilibrium of such networks imposes extensive local constraints on circular arc polygon geometry. These constraints provide the foundation of VMSI which, starting with images of epithelial monolayers, simultaneously approximates both tissue geometry and internal forces, subject to the constraint of equilibrium. We find VMSI to be more robust than previous methods. Specifically, the VMSI performance is validated by the comparison of the predicted cellular and mesoscopic scale stress with the measured myosin II patterns during early Drosophila embryogenesis. VMSI prediction of a mesoscopic stress tensor correlates at the 80% level with the measured myosin distribution and reveals that most of the myosin activity in that case is involved in a static internal force balance within the epithelial layer. In addition to insight into cell mechanics, this study provides a practical method for nondestructive estimation of stress in live epithelial tissue.

Popular Summary

Morphogenesis, the process by which a developing embryo shapes itself, requires precise orchestration of growth and cell rearrangements. Mechanical interactions between cells play a prominent role in controlling morphogenesis, but their study has been impeded by the difficulty of directly measuring mechanical stress in live tissues. One promising solution is based on highly developed fluorescent imaging of live tissues and aims to infer mechanical stress from the observed 2D geometry of cells. Current algorithms for this sort of inference are hindered by image noise and the quality of the underlying image analysis. Here, we present a new approach to inferring mechanical stress from tissue images that overcomes these hurdles.

Our approach is based on a novel parametrization of cellular packings (in epithelial tissues) in terms of polygonal tilings with circular arc edges. These polygons represent the equilibrium geometry of an array of cells, allowing us to identify the observable constraints on the tiling geometry that are imposed by the force balance needed for mechanical equilibrium. Implementing these constraints into the image analysis of single layers of epithelial cells, we formulate a robust algorithm that can infer cell-scale stress over the whole surface of the tissue. Testing our algorithm on data for the embryonic development of fruit flies, we find that the inferred stress strongly matches the measured global pattern of myosin II, a molecular motor known to generate mechanical forces within tissue.

Our algorithm should be immediately useful to researchers studying the mechanics of tissue and organ development, providing a practical and nondestructive approach to measuring stress in live tissue.

Figure 1

Circular arc polygonal (CAP) tiling and its dual tension triangulation. (a) Circular arc polygons provide an approximate representation of the equilibrium geometry of a cell array, in which edges are represented as circular arcs and vertices are loci where three cells meet. ${\varphi }_{j,\beta }\sim \pi$ implies that $T_{ij}\ll T_{ik},T_{il}$. (b) Force balance at vertex $i$ requires the three tension vectors tangent to each edge to sum up to zero, which defines a local triangle. As adjacent vertices in (a) share edges, the tension triangles of each vertex form a triangulated surface—the dual representation of force balance—with each triangular face corresponding to each vertex. Pressure differentials between cells result in a nonflat dual triangulation. The height above the $xy$ plane is represented as gray scale. (c) In mechanical equilibrium, the curvature of edge $\alpha \beta$ is controlled by pressure differences between adjacent cells, $p_{\alpha }-p_{\beta }$, and the interfacial tension $T_{\alpha \beta }$ by the Young-Laplace law [40].

Figure 2

Determination of equilibrium CAP tilings. (a) Fitting an equilibrium CAP tiling to imaging data. The best fit to raw data (shown in gray scale) by an equilibrium CAP tiling found by minimizing Eq. (5) is shown in red (with vertex positions in blue). The CAP contour is close, but distinct from the result of standard cell segmentation (see Sec. 7), here shown in green. (b) An example of vertex position ${\mathit{r}}_i$ derived from the generalized Voronoi construction—i.e., generating points ${\mathit{q}}_{\alpha }$ displaced to height $z_{\alpha }$ above the plane and a metric rescaling $p_{\alpha }$. Each CAP edge $\beta \gamma$ is the locus of points with a fixed ratio of distances to a corresponding pair of generating points—e.g., $\sqrt{p_{\beta }}{d}_{\beta }=\sqrt{p_{\gamma }}{d}_{\gamma }$ in the above image. The vertex is defined by the intersection of all three arcs. Parameters of this generalized Voronoi construction (including $p_{\alpha }$) provide a convenient set of reduced variables used in the minimization of Eq. (5).

Figure 3

Validation of VMSI for force inference using synthetic data. (a) Comparison of actual and inferred tensions based on the matrix-inverse-type inference method introduced in Ref. [29] made overdetermined using measured edge curvature data, as in Ref. [35]. 10% noise was added to the synthetic input data (see Appendix pp4 for details). $\rho$ denotes Pearson’s correlation coefficient. (b) Same as (a) but using the VMSI on the same synthetic data. Note the significant increase of the correlation coefficient.

Figure 4

Stress inference results compared to measured myosin distribution for the germ-band extension in Drosophila embryo. (a),(b) Raw in toto images showing fluorescent-labeled membrane [mCherry, shown in (a)] and myosin [Sqh-GFP, shown in (b)] obtained using light sheet microscopy and unrolled into the plane using imsane [47]. Time corresponds to roughly 26 min after cephalic furrow (CF) formation. (c) The correlation between myosin and inferred tensions comparing the matrix-inverse-based method defined in Ref. [29] and the VMSI algorithm during the first 40 min of GBE. (d) An example of force inference, combining two overlapping image frames used for the analysis, shown directly on the embryo and on a cylindrical projection with the ventral line cut and mapped onto the top and bottom edges of the image with the dorsal side along the midline. Inset: Color-coded relative edge tension (in arbitrary units, arb. units) inferred by the VMSI algorithm. The inferred mechanics displays stress cables as expected for the lateral region of the embryo at the time considered. (e) The mesoscopic anisotropy of myosin, 26 min post CF formation, shown on the embryo and its cylindrical projection. Anisotropy is largest in the lateral region. The principal axis of myosin in this region points along the DV axis. (f) The spatial distribution of the normalized rms difference [$\mathrm{\Delta }(r)$] between the VMSI inferred stress tensor and the measured total myosin tensor (at 26 min post CF). Both tensor fields are represented as ellipses for direct comparison. The large discrepancy ($\mathrm{\Delta }\sim 0.8$) found in the vicinity of the anterior and posterior poles (mapped to the left and right edges of the cylindrical projection) is due to the poor imaging of the poles. The large discrepancy in the center of the dorsal region is real and can be explained by the difference between the total and “balanced” myosin distributions as explained in the text.

Figure 5

Balanced myosin versus total myosin. (a) The magnitude of the myosin tensor measured in Ref. [36] 33 min after the formation of the ventral furrow. (b) Same for the balanced component myosin tensor, defined by Eq. (13), at the same time. Note that the balanced myosin is approximately constant along the DV axis in contrast with the total myosin displayed in (a). (c) The magnitude of the inferred stress tensor at the same time point. Note greater similarity with (b) than with (a). The dotted gray box shows the region that was amenable to image segmentation. (d) The average misalignment Eq. (10) between inferred stress and myosin tensors [balanced (total) displayed as solid (dashed) orange lines]. The balanced fraction of myosin is shown in blue. As expected, the inferred stress tensor provides a better approximation for the balanced component of myosin tensor then for the total.

Figure 6

Correlation of inferred stress and cell-division axis. (a) Confocal image of the Drosophila lateral ectoderm, near the cephalic furrow, during the late phase of germ-band extension (GBE) at the onset of cell divisions, approximately 45 min after the formation of the cephalic furrow. (b) An overlay of inferred tensions on the CAP array. Note that a low-tension edge (dark red) meets two high-tension edges; the latter meet at an angle close to $\pi$ while low-tension edge is nearly perpendicular. (c) Same as (b) but with the average stress tensor for each cell plotted as an ellipse. The major or minor axis of the ellipse corresponds to the principal axes of stress. (d) Comparison of the cell cleavage axis for cell division events (within mitotic domains 6 and 11 [49]) and the principal (extension) axis of inferred local stress (using different methods). Stress for mitotic cells was estimated 5 min before the registered time of division. While VMSI-based analysis finds significant correlation, predictions based on the extended “matrix inverse” [29, 35] or cell elongation axis are not distinguishable from random (shown as dashed gray line).

Figure 7

An illustration of the definitions used in the procedure to find an initial estimate of $\{p_{\alpha },{\mathit{q}}_{\alpha }\}$, plotted on top of simulated data. The segmentation is shown in red while the idealized edge (to be fitted) is shown in light blue. The tangent vectors ${\mathit{\tau }}_{i,j}$ are estimated from the segmentation. The minimization of ${\mathcal{E}}_I$ [Eq. (c1)] finds $\{p_{\alpha },{\mathit{q}}_{\alpha }\}$ minimizing the overlaps ${\widehat{\mathit{t}}}_{j,i}·{\widehat{\mathit{\tau }}}_{j,i}$.

Figure 8

(a) An example plot of the image of weighted distance from each generating point produced using matlab’s bwdist function and the algorithm described. The locations of the generating points ${\mathit{q}}_{\alpha }$ are shown as blue dots. (b) The resultant CAP network. The segmentation output from matlab’s watershed algorithm is used to infer the triangulation topology. This is then passed into our forward equations to generate the exact positions of vertices and curvatures of edges. The result, using the original values for the generating points and the measured triangulation topology, is overlayed in blue over the original distance map.

Figure 9

Correlation between inferred tensions and pressures against the known synthetic values for three mechanical inference schemas—CVMI (blue), our implementation of the cellfit method [35] (orange), and the inverse used within Ref. [32] (green)—as a function of the contribution of pressure to mechanical balance (a) and the strength of random perturbations to vertex position (b). Panel (a) represents a simulated lattice with constant noise of 0.03 (strength of perturbation is normalized to the average edge length). Similarly, panel (b) is taken at a finite average pressure differential of 0.03 (average ratio of edge length to radius of curvature). All algorithms were run on the same synthetic datasets. Each curve depicts the mean and variance of 100 networks of size $\sim 120$ cells generated with equivalent parameters as described in the text.

Figure 10

Two examples of simulated CAP networks used in the in silico verification presented in Appendix pp4 with varying average pressure differences. Networks (a) and (b) were generated with varying degrees of intracellular pressure differentials. The average edge curvatures, normalized to edge length—i.e., $R_{\alpha \beta }/r_{\alpha \beta }$—for (a) and (b) are 0.01 and 0.20, respectively.

Figure 11

(a) A synthetic spherical embryo of $\sim 3000$ cells, plotted both in the embedded space, as well as the cylindrical unwrapping analogous to the embryonic data shown in the main text. Simulated interfacial tensions are plotted as a heat map in both domains. (b) Graphical depiction of the process used to infer stresses in the tangent plane of a small region of cells. All cells are segmented in 3D using imsane [47]. Each local cellular patch is projected from the 3D embedding space onto the 2D plane that minimizes the sum of squares of deviation. The resulting projection is used as an input into the VMSI method. (c) Scatter plot between the inferred tensions using the work flow outlined and the known tensions shown in (a). The inset displays the dependence on the number of patches used to cover the sphere. As was expected, correlation is monotonic with sampling resolution.

Figure 12

(a) A confocal image of the patch of lateral ectoderm during GBE comparing vertex positions obtained by segmentation (in blue) and the position obtained by VMSI (in red). (b) The plot of residual [Eq. (5)] normalized by the average bond length over the time course of the movie analyzed. The inset displays the distribution of values. The VMSI begins to make qualitative errors during the late stage GBE, in agreement with our observation that myosin becomes unbalanced.

Figure 13

Distributions of cell area (a) and cell eccentricity (b). Solid lines indicate means while surrounding outline indicates one standard deviation. We observe that the onset of the “rounding up” of the cell and the start of area growth occurs at 5 min before division.