## Abstract

Although the importance of cellular forces to a wide range of embryogenesis and disease processes is widely recognized, measuring these forces is challenging, especially in three dimensions. Here, we introduce CellFIT-3D, a force inference technique that allows tension maps for three-dimensional cellular systems to be estimated from image stacks. Like its predecessors, video force microscopy and CellFIT, this cell mechanics technique assumes boundary-specific interfacial tensions to be the primary drivers, and it constructs force-balance equations based on triple junction (TJ) dihedral angles. The technique involves image processing, segmenting of cells, grouping of cell outlines, calculation of dihedral planes, averaging along three-dimensional TJs, and matrix equation assembly and solution. The equations tend to be strongly overdetermined, allowing indistinct TJs to be ignored and solution error estimates to be determined. Application to clean and noisy synthetic data generated using Surface Evolver gave tension errors of 1.6–7%, and analyses of eight-cell murine embryos gave estimated errors smaller than the 10% uncertainty of companion aspiration experiments. Other possible areas of application include morphogenesis, cancer metastasis and tissue engineering.

This article is part of the themed issue ‘Systems morphodynamics: understanding the development of tissue hardware’.

## 1. Cellular forces

The importance of cellular forces to morphogenesis, wound healing and disease is now widely recognized [1–5]. Experimental studies have shown that irregularities in these forces can give rise to developmental defects and other malformations [6], and computational models have shown that alterations as small as 20% can affect clinical outcomes [7].

In many settings, these driving forces can be treated as equivalent interfacial tensions *γ _{i}* along cell–cell and cell–medium boundaries [8,9] (electronic supplementary material, figure S1), a concept considered as early as the 1960s [10,11] but quickly dropped in favour of the differential adhesion hypothesis [12]. Computer simulations of sorting and other cell movements carried out by the authors, however, showed that a wide range of cell and tissue movements were in fact driven largely by interfacial tensions [13–15] (to which cell–cell adhesions make a counteracting contribution [8]), leading to a new paradigm [2,16–20]. Information about these tensions can be obtained through a variety of experimental techniques, that in generally decreasing length scale include: pipette aspiration, deformable substrates, engineered droplets, laser ablation, atomic force microscopy, optical tweezers, magnetic cytometry, and fluorescence resonance energy transfer (FRET) microscopy [18,21–31]. Unfortunately, these techniques can be expensive, time consuming, invasive, or destructive, and all except substrate deformation, FRET and droplet methods provide information about only a single cell surface at one moment in time.

When force inference techniques for two-dimensional systems [32–36] entered the scene in 2010, they made spatio-temporal maps of cellular forces possible for any available image, including historical ones. Force inference techniques, like CellFIT which is here denoted as CellFIT-2D to distinguish it from the present three-dimensional version, are based on two primary mechanical assumptions: that the tension-carrying boundaries that meet at any triple junction (TJ) [36] produce mechanical equilibrium there, and that these tensions are constant over the span of any given boundary or interface (electronic supplementary material, movie S1). Force inference techniques in this class require only that cell edges be visible, they involve no mechanical intervention, and they produce no additional damage beyond that caused by microscopy. Data from CellFIT-2D revealed tension variations around the perimeters of individual cells, differences between cells in a single population and between populations, elevated tensions along inter-population boundaries [36], and purposeful temporal variations [17].

On the strength of this and other evidence, we propose that cells move and organize primarily by gradual and carefully orchestrated changes in interfacial tensions and protrusion contractions, with the latter acting along the interface between cells and serving as a special case of these tensions. Finite-element studies have shown that force imbalances at TJs cause those junctions to move until TJ angles produce an equilibrium configuration (electronic supplementary material, movie S1) or until a high-tension boundary shortens enough to produce a neighbour change [15]. If these forces were to change systematically and gradually with time, a series of carefully controlled cellular and tissue morphologies could be produced. Modelling studies of whole embryos [37] found that physiological motions were best matched by time-varying driving forces, but the reason for this finding was unclear when those studies were carried out. The finding was contrary to the terraced nature of these forces tacitly assumed by modellers and biologists. In the context of embryogenesis, the context in which the most active research on cell mechanics and morphogenesis was occurring, the underlying mindset was that gene expression or some other process set up conditions for the next developmental step and then it moved ahead on the basis of those conditions. The revised understanding for which we here argue, is that regulatory networks and mechanobiological circuits play a much more active role in controlling the driving forces from moment to moment. If cellular forces do indeed change with time as these studies suggest, then techniques like CellFIT—which can take snapshots of the forces in all of the cells visible in a field of view and do so at multiple sequential times—become particularly useful.

Other kinds of forces—such as intracellular pressures and deviatoric stresses associated with cytoplasm deformation—could be included, as was done in video force microscopy (VFM) [17]. Line tensions along the curvilinear interfaces between trios of cells could also, in principle, be included, but we suspect that their contribution is minimal. The picture of gradual force changes we present should not be confused with models in which increasing forces are applied to elastic systems. Elastic systems would spring back to their initial state should the driving forces be eliminated, though real embryos and their tissues do not. Instead, we contend that embryos and their tissues behave primarily in a plastic manner [38], with their morphogenetic deformations being essentially irreversible due to cell neighbour changes. An apparent exception to this statement is the small elastic component of deformation that can occasionally be observed, as when momentary forces are applied, and a reversible response occurs only because cells do not have time to change neighbours and lock the temporary geometry in place. We have argued elsewhere [28] that viscous, visco-elastic and elastic tissue behaviour can all be produced by cells that have constant edge tensions and viscous cytoplasm.

Interestingly, the TJ force balances on which CellFIT is based do not depend on the mechanical characteristics of the cell membrane system, including whether it is elastic, viscous, governed by rate constants, affected by endocytosis, influenced by adhesion systems or altered by cortical components. In a sense, CellFIT operates one level up from these important details, and simply provides the total relative tension acting along any particular interface, without regard to how it is generated. In addition, intracellular pressures and any stresses from volumetric viscous or elastic components seem to act primarily normal to the interface (except perhaps during laser ablation experiments, which violate other CellFIT assumptions) and likewise, do not affect the TJ force balances [36]. These fortuitous circumstances make CellFIT applicable to a very wide range of biological applications. CellFIT would not be applicable, however, to interfaces where cells are adhered to a substrate, because the forces carried by that substrate would generally not be known. If such forces were known, however, governing equations that include them could be constructed [5] and a suitably adapted version of CellFIT applied. It would also not be applicable to interfaces that contain significant spatial force variations along individual cell faces. One might hope that the quality measures that form part of CellFIT would aid in identifying situations of this kind.

In this article, we outline a series of steps (figure 1) that allow CellFIT-2D to be extended so that cellular forces in three-dimensional aggregates and tissues can be inferred from serial sections. When applied to synthetic sections generated using Surface Evolver (SE) [41], the tensions calculated by CellFIT-3D had errors as small as 1.6%. CellFIT-3D was then applied to eight-cell compaction-stage mouse embryos, and the tensions found had estimated errors lower than the 10% uncertainty associated with accompanying aspiration experiments.

## 2. CellFIT-3D

As in the finite-element models we have used to study embryonic cells for many years and in VFM and CellFIT-2D, we here assume that the subcellular forces generated by various structural protein and adhesion systems can be deemed to generate an equivalent tension tangent to the cell membrane [13,15,17,36,42] (electronic supplementary material, figure S1). The total tension along any given cell–cell or cell–membrane interface i is denoted *γ _{i}*, and it is assumed to be specific to that interface, whether an edge in two dimensions or a face in three dimensions [43].

Furthermore, whenever three cells meet at a particular point and a cutting plane is constructed normal to the TJ between those cells at that point, the vector sum of the membrane tensions in that plane must add to zero. Studies of two-dimensional cellular systems showed that the membrane angles at TJs are unaffected by intracellular pressures or deviatoric stresses associated with reshaping of the viscous cellular cytoplasm [36]. Consequently, TJs can be analysed and used to calculate membrane tensions *γ _{i}* without reference to pressure or viscous forces.

Cellular pressures were not calculated as part of CellFIT-3D, but one could presumably calculate them, as is done in two dimensions, after the interfacial tensions *γ _{i}* are determined [36]. In three dimensions, the local principal curvatures

*k*

_{1}and

*k*

_{2}of cell–cell or cell–medium interfaces change with position unless they are spherical. However, the mean curvature (

*k*

_{1}+

*k*

_{2})/2 would be constant throughout any one of these surfaces if it carries a constant isotropic tension

*γ*, and the pressure difference Δ

_{i}*P*across it would be given by the Young–Laplace equation 2.1

_{i}To determine the spatially varying principal curvatures or even the mean curvature of a surface based on cuts through it, however, is a problem beyond the scope of the present study. For the present, there is probably no reason to calculate intracellular pressures or their differences, as their relevance to development is still unclear and methods to verify them experimentally are limited.

Because of the geometric complexity of three-dimensional systems, it is valuable to distinguish between the curvilinear TJs or triple edges (TEs) that arise between trios of contacting cells or between pairs of cells and the medium, and the points at which those junctional curves pass through individual confocal or other sections. We use the term ‘triplet’ to describe the latter locations where three, but occasionally more, cell boundaries are seen to converge (as in the boxed area of figure 1*a*). The cell membranes immediately adjacent to the point are considered part of the triplet.

A series of image processing and geometric and mechanical analyses steps are required to implement CellFIT-3D, and they are outlined here.

### (a) Image capture

As we will show, the quality of a CellFIT-3D analysis depends on image resolution and section spacing. Reliable tension information can be obtained from data like figure 1*a*, a portion of a 512 by 512 pixel image from a confocal stack of an eight-cell murine embryo. There were approximately 125 pixels across the diameter of a typical cell, and individual cells were transected an average of 18 times as a result of the 2 µm section spacing.

### (b) Image enhancement

The cell boundaries in these images may seem well defined, but at the pixel level (figure 1*b*) where computational algorithms work, noise, gaps and other anomalies become apparent. An ideal image would have accurate, gap-free white outlines a single pixel wide and a pure black background. A number of image processing algorithms for amending the images toward this ideal were appraised, and an edge enhancing coherence filter [39] that smoothed the boundaries (figure 1*c*) and largely closed the gaps (figure 1*d*) was selected. Relevant computer code for this step and a number of the others, as well for as the validation procedures, is available through the indicated references.

### (c) Image segmentation

Cells in the processed images were segmented (figure 1*e*) using a watershed algorithm called SeedWater [40,44]. The algorithm treats the bright pixels at cell edges as if they defined mountain range heights and it begins to flood each mountain-surrounded region from a relatively dark pixel or low-lying point called a seed point. The software chooses these seed points, but users can manually adjust them. Boundaries are found indirectly, based on where the flooded zones contact each other or where they end at the outer mountain ranges. This approach is less sensitive to intracellular noise—spurious hills in low lying regions—and partial gaps than are standard edge tracing techniques. The resulting boundaries may, however, still contain pixel-scale positional noise (figure 1*f*) and other imperfections.

### (d) Grouping of triple edges

In three dimensions, individual TEs may appear as triplets in multiple images and it is necessary to associate these occurrences with each other. As with boundary identification, it is preferable to work with cell cross-sectional areas, and we use a cell-grouping algorithm that considers area overlap across adjacent images, centroid collinearity and other geometric features [45] to group cross-sections associated with the same cell (figure 1*g*) and in turn to identify triplets belonging to the same TE. If sections are spaced sufficiently close, seed-point proximity can be used to group the outlines. Grouping also makes possible three-dimensional cell models for visualization and verification of mesoscale geometry and topology.

### (e) Mesh generation

The pixelated boundaries of sectioned cells may have complex shapes with reversing curvature and may contain noise. In order to obtain good triplet approach angles, and informed by CellFIT-2D [36], we fit cubic splines to the boundary pixels along each edge and use these splines to generate uniformly spaced mesh points (figure 1*h*).

### (f) In-plane angles

CellFIT-2D showed that accurate angles are crucial to the tension calculations, and we found that membrane endpoint directions could be estimated well by fitting separate circular arcs to the last four or five mesh points on each end (figure 1*i*), as in two dimensions [36]. Standard software exists for arc fitting, the fit is coordinate indifferent, it provides a second-order approximation to the shape of the boundary terminus, and it attenuates noise.

Triplets that appeared in at least three successive images were identified algorithmically, and shown graphically in their corresponding sections (figure 1*j*). This strategy was more efficient than manually vectorizing triplets of interest. The graphic consisted of a circle at the calculated triplet location and three vectors with circles at their ends in the calculated approach directions. The graphics were automatically overlaid on the raw images and manual angular adjustments were made as needed to best capture the final approach angle at the junction (figure 1*k*).

### (g) Local dihedral angles

To convert the in-plane angles defined by these graphic triplets to true dihedral angles, splines (shown as orange curves in figure 1*l*) were constructed through sets of grouped triplets. Construction of a reliable spline required a TE to appear in at least three images. Local dihedral planes (figure 1*m*) were constructed normal to the spline at each triplet, and the boundary vectors were mapped onto them (white arrows) by projection normal to the planes. Finally, a pair of equilibrium equations was constructed for two arbitrary orthogonal directions in each plane. These equations define the ratios that the three boundary tensions *γ _{i}* must have for that section of the TE to be in equilibrium.

### (h) Average dihedral angles

The equation pairs generated at various triplets along any one TE can give different force ratios, and the quality of the equations varied. For example, if an image plane is strongly oblique to a TE, the fluorescent signals from one or more of its boundaries may be spread out and grainy. In addition, errors in calculated boundary directions can amplify when projected onto highly tilted dihedral planes, making information from these planes less trustworthy. When more than three equation pairs are available, they are automatically checked for consistency and equations from the ends of the spine, where the dihedral planes tend to be more tilted, eliminated if appropriate. Increasing the number of sections per cell improves the geometry of the spline and increases the number of available equation pairs. However, even when there are seven slices per cell on average (orange curve in electronic supplementary material, figure S3), approximately two thirds of the TEs still appear in two or fewer sections. To amalgamate the multiple equation pairs associated with any one TE into a single set of ratios *γ _{i}*:

*γ*:

_{j}*γ*, we recommend a least-squares ratio solver [46].

_{k}In contrast to these steps, one could construct point clouds from the meshing points associated with each boundary. One could then fit mathematical surfaces to each cloud, or a portion thereof, and use mathematical descriptions to calculate the TE contact angles. Unfortunately, the shapes that arise are complex, often having reversing concavity and other challenging features, making this approach impractical [47]. One could also take in-plane angles as surrogates for dihedrals and apply CellFIT-2D to individual sections. However, doing so ignores the oblique angles of typical TEs to these sections, and produces tension errors of 50% [47]. Alternatively, for specific TEs, one might use the in-plane angles of the triplet apparently most normal to that TE [21], but numerical tests show that even this approach introduces approximately 5% additional error in the dihedral angle equations, in addition to foregoing the statistical benefits of multi-triplet averages.

### (i) Assembly and solution of tension equations

The previous step produces either a pair of equilibrium equations or two tension ratios for each of the *n* TEs having a sufficient number of useable triplets. In either case, these equations can be assembled into a single matrix equation of the form
2.2where **G** has 2*n* rows and *m* columns, and *m* is the number of tensions *γ _{i}* that appear in one or more of the 2

*n*force-balance equations, and whose values will in due course be found. These equations correspond exactly to the equation pairs that arise in CellFIT-2D [36], and their assembly, solution and evaluation are identical, unaffected by the dimensionality of the host space.

A unique solution to this homogeneous and overdetermined system can be found by constructing and solving the least-squares system
2.3or
2.4where
2.5a choice that selects the particular set of scaled *γ*s having a unit mean.

### (j) Solution evaluation

The toolkit developed for assessing CellFIT-2D solutions was found to apply equally well to CellFIT-3D. The overall quality of the equations is assessed using the condition number of **G***, a value equal to the ratio of its largest and smallest eigenvalues. This ratio portends the degree to which error in **n** is amplified in **γ*** [48], and so indicates the sensitivity of the calculated tensions to angular and other types of error. The condition number will vary with matrix size, and magnitudes substantially higher than those shown here and in other CellFIT analyses [36] may signal structural problems with the equations and reduced solution accuracy.

In contrast, the solution residual
2.6provides a measure of how well the pair of equations associated with each TE is satisfied by the least-squares solution **γ**. The paired components of **r** indicate the degree to which each TE is out of balance and they are useful for identifying triplets that may have been inaccurately placed during the automated meshing steps and may require manual adjustment. Finally, the scaled cofactors of **G***, also known as standard errors [36], indicate the confidence levels associated with individual tensions.

## 3. Validation

To assess the CellFIT-3D algorithms, they were tested on synthetic sections generated using Surface Evolver [50] (figure 2 and electronic supplementary material) for which ground truth tensions were known. The synthetic data were designed to approximately match the 18 slices per cell of the murine embryo data used here, and the image data were processed as outlined in figure 1. The error in the calculated tensions was found to depend on slice spacing and image resolution, but was insensitive to the number of cells analysed, whether those cells were surrounded by medium or cut cells, and the range of the tension values *γ*_{I} within the model.

Increasing the number of slices per cell increases the number of dihedral planes per TE (electronic supplementary material, figure S3) and allows strongly sloped planes near their ends to be ignored. Seven slices per cell was found to give tension errors of 3.2% in isolated aggregates of eight cells, while 14 slices per cell, often not difficult to obtain experimentally, reduced that error to only 1.6% (figure 2*c*). In all cases, the number of equations was more than adequate to solve for the unknown tensions (see electronic supplementary material).

Portions of isolated synthetic aggregates of eight and 50 cells were eroded away in order to generate a range of aggregate sizes and conditions—from those fully surrounded by medium, to those with mixed medium and cut cell edges, to fragments containing no complete cells. For aggregates with 14 slices per cell, the errors depended primarily on the image resolution, with 400 and 200 pixels per cell diameter generating errors of 3% and 8%, respectively.

Error was then added to the averaged dihedral angles, and figure 3 shows how adding angular noise to the average dihedral angles affects tension error, residual and standard error. The angular noise was randomly assigned using a Gaussian distribution. Conveniently, all three track similarly with noise, allowing the latter two—which can be calculated from the governing equations associated with any given dataset—to indicate the likely angular error (input noise) and tension error (output uncertainty). A mean angular noise of 5°, a value consistent with typical manual digitizing errors of 5° [36], for example, caused tension errors of 7.1% and residual and standard errors of 9.4% and 5.2%, respectively.

## 4. Application to murine embryos

CellFIT was then used to investigate eight-cell murine embryos (figure 4*a* and electronic supplementary materials) undergoing compaction, a process considered to be driven by surface tensions [21,51]. Interfaces were typically seen in at least five slices. Some manual adjustments were necessary in places where the automated meshing was inaccurate. For time A, 90 min after the last four-cell stage blastomeres divided, the CellFIT-3D tensions had residuals of 0.031 and standard errors of 5.4% (figure 4*b*), while the values for time B (figure 4*c*), 150 min later, were slightly higher. The corresponding boxes in figure 3 suggest input errors corresponding to noise levels between 2 and 6 and tension errors of 3–9%.

The surface tensions of all eight cells (figure 4*a*) were measured by pipette aspiration [21] and the experimental error estimated from repeated tests was 10%. Measurement of tensions on interior surfaces is more challenging, but could be done using optical tweezers or laser ablation, though the latter is destructive and precludes the taking of force measurements at multiple times or locations. The root mean square (RMS) errors between the CellFIT-3D and aspiration values at times A and B (figure 4*b*,*c*) were 14.6% and 11.6%, respectively, values consistent with aspiration errors of 10% and standard errors of 5.4% and 6.4%, for times A and B, respectively. The surface tensions of individual cells changed between times A and B, and the time delay between sequential aspiration measurements may have contributed additional error. In contrast, CellFIT typically draws its data from the much narrower window of time required to collect a single set of sections. Several other eight-cell murine embryos were analysed with CellFIT, and similar results were obtained.

## 5. Discussion

This article demonstrates that cellular forces (interfacial tensions) in three-dimensional aggregates can be inferred from serial sections. In principle, any image sources including confocal sections or histological sections could be used, provided that cell boundaries are well defined and there are at least 6–8 slices per cell. Furthermore, tension errors can in principle be as low as 1.6%, considerably lower than the variability of typical experimental tensometry techniques. Because the assembled CellFIT equations are overdetermined, input and solution quality can be assessed using residuals and standard errors, and TEs with apparently discordant angles or other anomalies can be identified and removed without jeopardizing the tension inference process. These quality measures could also serve to identify situations where membrane tensions are not uniform as assumed in CellFIT, as could be the case should undulating membrane geometries arise, or localized tension sources give rise to nonuniformities in membrane tension.

CellFIT-3D builds on its two-dimensional counterpart, and its three primary challenges—obtaining cell outlines, connecting triplets in successive images, and calculating true dihedral angles—can be overcome using relatively straightforward algorithms. Like its counterpart, CellFIT-3D provides only tension ratios, and should scaled values be required, direct force measurements must be obtained. Fortunately, these ratios are often all that is required [21], as it is these ratios that govern cell shapes and motion patterns, with cytoplasm viscosity and force strength determining motion rates, only [52].

Force inference methods are possible because the relative angles at TEs in equilibrium depend uniquely on membrane forces, even though other mesoscale features (see supplementary text) do not [53,54]. Should four or more cell edges impinge on a single junction or there be a rosette [55], the number of unknowns exceeds the number of equations by more than one, and a unique tension ratio does not exist. However, if the equations and unknowns associated with such junctions are assembled with those from other suitably connected TEs, a unique and trustworthy solution should still be possible [36]. In principle, CellFIT-3D could be extended to calculate intracellular pressures, as in CellFIT-2D [36], but calculating the required surface curvatures is much more difficult in three dimensions, there is currently little experimental interest in these pressures, and their consequences for cell movement remain unclear. The effects of viscous cytoplasmic forces could also be incorporated easily, with their contributions to TJ force balances calculated in a manner similar to that used in VFM.

If the cellular forces that drive cell movements, mesoscale assembly and bulk tissue motions turn out to be time-varying, as recent evidence suggests, then techniques like CellFIT that can provide force maps over a whole field of view and at closely spaced times become particularly useful. Hopefully, in time, the technique presented here will provide new insights into the movements of single or small groups of cells during morphogenesis and diseases like metastasis, formation of mesoscale structures such as acini in various organs and engineered tissues, and bulk tissue movements.

## Authors' contributions

The force inference concepts were developed and implemented by G.W.B., J.H.V. and A.E.; the experiments were conducted by J.-L.M. in the lab of T.H.; the Surface Evolver models were developed by S.C. and J.H.V.; the paper was written primarily by G.W.B.

## Competing interests

We have no competing interests.

## Funding

Funding was provided by a Natural Sciences and Engineering Research Council of Canada (NSERC) Discovery Grant to G.W.B.

## Acknowledgements

Animal studies were performed at the European Molecular Biology Laboratory in accordance with Federation for Laboratory Animal Science Associations guidelines and recommendations. Surface Evolver was written by Kenneth Brakke, Susquehanna University, and SeedWater was written by David Mashburn under the supervision of Shane Hutson, Vanderbilt University.

## Footnotes

One contribution of 13 to a theme issue ‘Systems morphodynamics: understanding the development of tissue hardware’.

Electronic supplementary material is available online at https://dx.doi.org/10.6084/m9.figshare.c.3691960.

- Accepted October 3, 2016.

- © 2017 The Author(s)

Published by the Royal Society. All rights reserved.