Physical modeling of cell geometric order in an epithelial tissue
See allHide authors and affiliations

Communicated by Phillip A. Sharp, Massachusetts Institute of Technology, Cambridge, MA, November 26, 2007 (received for review September 24, 2007)
Abstract
In multicellular organisms, cells pack together to form tissues of intricate and well defined morphology. How such cellpacking geometries arise is an important open question in biology, because the functionality of many differentiated tissues depends on their reliable formation. We show that combining adhesive forces due to E and Ncadherin with a quantitative description of cell membrane elasticity in an interfacial energy model explains not only the qualitative neighbor relations, but also the detailed geometry of a tissue. The characteristic cellular geometries in the eyes of both wildtype Drosophila and genetic mutants are accurately reproduced by using a fixed set of few, physically motivated parameters. The model predicts adhesion strengths in the eye epithelium, quantifies their role relative to membrane elasticity, and reveals how simple minimization of interfacial energy can give rise to complex geometric patterns of important biological functionality.
Differentiated biological tissues consist of a variety of cells of distinctive morphology, with the cell shapes often reproduced with astonishing accuracy between individuals and across species (1). Examples for the functional importance of well defined cell geometry in multicellular tissues abound: hexagonal packing in the vertebrate lens minimizes light scattering from plasma membranes and increases transparency (2), and regular packing of sensory hair cells in the inner ear allows the precise alignment of stereocilia bundles and maximizes sensitivity to displacement (3). Among the mechanisms that organize cell geometry, differential adhesion, that is, adhesive energies varying from cell to cell, has been shown to play an important role for cell sorting and distribution in tissue (4), although differential adhesion by itself is not sufficient to quantitatively explain all observed morphologies (5). Other studies have suggested that membrane tension terms and/or contractile forces on cells exerted by microfilaments in the cytoskeleton should also be taken into account (6). Depending on how such forces are modeled, they may arise from complex interactions in the entire cell body (7) or be an effective addition to membrane forces (6). Here, we restrict the focus to interfacial energy terms only and go beyond a qualitative explanation of cell sorting by using a stringently defined interfacial energy functional to determine equilibrium shapes. Such a model shows that the morphology can be explained quantitatively as a simple consequence of a small set of predetermined membrane and adhesion parameters.
We have examined the Drosophila eye as an example of a complexpatterned epithelial tissue. The retina of the eye is a monolayer epithelium composed of a sheet of adherent cells. The sheet is organized as a hexagonal array of repeating multicellular units called ommatidia (Fig. 1). Each ommatidium consists of 20 cells arranged in a precisely defined pattern, including 8 photoreceptor neurons and 12 accessory cells (8). Four cone cells adhere together to form a transparent plate that acts as both the floor of the simple lens and the roof of the underlying chamber of photoreceptors (Fig. 1). Two primary pigment cells optically insulate each conecell group, and together with the conecell group, they form the “core” structure of the ommatidium, whereas secondary and tertiary pigment cells form the ommatidium border (“frame”). This pattern is repeated >800 times within each retina. The cone cells form an “aperture” for focused light to be transmitted from the lens to the photoreceptors, a structure highly preserved across invertebrate species (9).
All of the cells adhere with each other through apicolateral contacts known as adherens junctions. It is known that cadherin molecules, which are localized at the adherens junctions, play an important role in cell geometric packing. All retinal cells express Ecadherin, and cone cells additionally express Ncadherin, which is specifically localized at cone–cone adherens junctions (10, 11). Loss of either cadherin in individual cells severely alters the shape of the mutant cells and their wildtype (normal) neighbors (10, 11).
The cone cell geometries of Drosophila ommatidia resemble the conformations of unconstrained soapbubble aggregates (10). Because soap foams are surface energy minimizers (12, 13), this suggests that conecell geometry is a consequence of the cells' interfacial energies (14). However, a cell membrane is not a uniformtension film (like a soap bubble), but is composed of a lipid bilayer and embedded molecules, and it is physically associated with an actin cytoskeletal cortex (15). Because a cell membrane is elastic in nature due to the bilayer and cortex, the net effect is a variable membrane tension (15–18). We therefore have developed an energy functional that combines biological membrane elasticity and adhesion, and we have used this functional to model cell geometries in the Drosophila retina.
Results and Discussion
To simplify the model, we make several assumptions. First, because the morphogenesis of the ommatidium pattern does not involve cell migration or sorting (19), we assume that cells do not rearrange neighbors when forming the final geometric pattern. Second, because the adherens junction is only ≈50 nm tall, the structure is twodimensional to a good approximation, and we can consider mechanical forces and energies in this plane only. Any forces out of plane have to fulfill separate, independent force balances that do not enter into the model. In this twodimensional model, cell membranes are characterized by their crosssectional perimeter lengths L _{i}, and the strain on cell i is governed by the deviation of L _{i} from the equilibrium length L _{0i}. Third, all cell membranes are assumed to share a uniform stretch modulus K _{A} (18). Fourth, the model assumes uniform coverage of E and Ncadherin molecules (where present) that reduces the surface energy per area by a fixed amount Γ_{E} and Γ_{N}, respectively. Fifth, we do not expect the crosssectional area of the cells to change, because such changes would involve bulk rearrangement of cytoplasm outside the narrowly defined adherens junction plane and also because, at least for small deformations, such changes are higherorder contributions to the total energy compared with perimeterlength changes [see supporting information (SI) Text ]. Then, a minimal model for an energy functional is where each cell i contacts its neighboring cells j. Because we are only interested in the shape (not the absolute size) of the pattern, K _{A} was scaled out and we introduced a relative membrane stretch Δ_{i} = (L _{i} − L _{0 i})/L _{0 i}. The relative cadherin adhesion strengths are γ_{E} = Γ_{E}/K _{A} and γ_{N} = Γ_{N}/K _{A}, active along the edges of length L _{ij} shared by cells i and j, if both cells express the respective cadherin (a homophilic interaction, represented by the product of Kronecker delta functions). This energy functional is supplemented with Lagrange multipliers (20) to enforce the constant crosssectional areas of the cells.
We quantified the observed shapes of cells within ommatidia by using a number of geometric characteristics. We measured these characteristics in Drosophila pupal retinas (48 h old; see Materials and Methods). To obtain unbiased, quantitative data, we digitized the cell interfaces from microscopic images by using algorithms from the Matlab imaging toolbox (21), and we fit the skeletonized interfaces to circular arcs (see Materials and Methods), consistent with our model that predicts uniform tension along each edge. Several geometric characteristics were measured. To simplify the measurements, we treated ommatidia as structures with fourfold symmetry (Fig. 2). Anterior and posterior cone cells are treated equivalently and are labeled C1, equatorial and polar cone cells are type C2, and primary pigment cells are type P. The geometric characteristics that were measured were the following: (i) the ratio of the width of a secondary cell attached to one P cell (w _{1}) to that of a secondary cell attached to two P cells (w _{2}); (ii) the center edge length L _{cen} between C2 cells; (iii) the distinct angles θ_{1}–θ_{4} between core edges; (iv) the relative crosssectional areas S _{i}/S _{tot}, where S _{tot} is the area of the entire ommatidium. The measured data, averaged over ∼100 measurements, are shown in Table 1, demonstrating that the geometry is extremely well defined, with only a few percent of variation in angles and lengths between different ommatidia.
We noted that both the skeletonized images and fits showed angles between primary and secondary pigment cells very close to 90°. Thus, the tension in secondary cells is much higher than in primary cells, resulting in a mechanical decoupling of the “frame” from the “core.” Our modeling naturally proceeds in three steps, fixing parameters starting from the frame and moving inward to the core.
In the first step, we optimized the membrane tension of secondary and tertiary pigment cells τ_{f} in the frame, to match the measured w _{1}/w _{2} (Table 1). Our second step analyzed the θ_{1} angle (Fig. 2), which governs the attachment of the cone cells to the P cells, and found that it strongly depends on L _{0P}, the equilibrium perimeter length of the P cells. The observed θ_{1} is reproduced by setting L _{0P}/(2(πS _{P})^{1/2}) ≈ 1.40 ± 0.01, independent of other parameters, confirming experimental evidence that the P cell crosssections in equilibrium (i.e., in P cells detached from their neighbors) are elongated (10). A similar fit of L _{0C1} and L _{0C2} is unnecessary, because detached cone cells show very nearly circular crosssections (10), so that we can set L _{0Ck}/(2(πS _{Ck})^{1/2}) = 1.00, where k = 1,2. The third and final step takes the core angles θ_{2},θ_{3}, θ_{4}, and L _{cen} as inputs to optimize the adhesion strength parameters γ_{E} and γ_{N} (Fig. 2). The membrane edge lengths are in quasistatic equilibrium, so that force balances between membrane edges exist wherever three membranes join. Enforcing zero total force in each vertex (force balance) relates the membrane stretch parameters to the four angles θ_{i} by Here, we introduced A _{i} ≡ cos(θ_{i}). Eqs. 2 – 5 form a linear system whose solution space can be analyzed by linear algebra (see Materials and Methods). We found that the experimental A _{i} lead to coefficients of O(1). This indicates that cadherin adhesion contributes significantly to membrane tension, being comparable to membranecytoskeleton elasticity. We also found that, if the angles θ_{i} were arbitrary, solutions would only exist for γ_{N} = 0 (see Materials and Methods). To allow for the observed finite γ_{N}, the angle cosines have to fulfill the constraint This condition does not depend on A _{1}, which was already determined. We inserted the measured angles into Eq. 6 and found that the relative error between the left and righthand sides of the equation is only ≈7%, giving a strong indication that we found the correct solution. Thus, θ_{4} is not an independent parameter, but determined by the other angles.
For the simulation of the global ommatidium shape, we used Surface Evolver (22), a program specifically developed to find minima of interfacial energy functionals such as Eq. 1 . Given an initial configuration of cells of random shapes and sizes (Fig. 3 A), and by using the constraints on crosssectional areas and equilibrium membrane perimeters, Surface Evolver found a consistent mechanical equilibrium configuration, dependent on the cadherin adhesion strengths γ_{E} and γ_{N}, but independent of the initial configuration (see SI Text ). The quality of this solution was monitored by combining the relative errors of the remaining quantities θ_{2}, θ_{3}, and L _{cen} in a penalty function that does not depend on subjective impressions of “realistic” shapes, The superscript SE denotes quantities obtained through Surface Evolver simulation, whereas the unsuperscripted quantities are those measured experimentally. The contour plot of the penalty function, F _{P} as a function of (γ_{E}, γ_{N}), found γ_{E} ≈ 0.025, γ_{N} ≈ 0.032 to be the adhesion parameters that gave the best modeling solution (Fig. 3 B). Strikingly, this level of adhesion contribution to surface energy is consistent with the relative level of adhesion energies experimentally measured for cadherins in cells (23). The optimized parameters γ_{E} and γ_{N} resulted in a model ommatidium that does exceptionally well in quantitatively reproducing the observed geometries (Table 1 and Fig. 3 C). Deviations from the optimum resulted in markedly deformed ommatidia (Fig. 3 B), whose characteristics were inconsistent with experimental measurements. In particular, for smaller values of γ_{E} and γ_{N}, a reorganization of cone cells became energetically favorable with the C1 cells being neighbors. This is never experimentally observed in normal ommatidia at this stage of development, although it is seen in ommatidia at earlier stages of development (8). It suggests that an increase in adhesion by E and Ncadherin might be an important factor in driving the morphogenetic reorganization of cone cells during pupal development.
We further validated the model by applying it to situations in which one or more cone cells specifically lack Ncadherin. With the parameters γ_{E}, γ_{N} unchanged, the only change in Eq. 1 is the absence of Ncadherin in specific cone cells, for which the corresponding δ_{i,N} are switched from one to zero. These modeling solutions compared well in qualitative terms with the cell geometries observed in ommatidia composed of comparable cone cells experimentally missing the Ncadherin gene (Fig. 4 A). More quantitative comparisons with data from these mosaic mutants were not possible because the number and quality of available mutant images were too low. In another crucial test of the model we varied the number of cone cells in each ommatidium, keeping all other parameters equal. When this was done, unique solutions were obtained that showed striking differences from the normal geometric pattern (Fig. 4 B). These solutions were very similar to the patterns observed experimentally in a Roi mutant background, where the number of cone cells varies from one to six. In the case of five Roi cone cells, a single geometry predominates, which is also observed as the modeling solution (Fig. 4 B). In the case of six Roi cone cells, three topologically different configurations are observed 97% of the time (119 of 123 examined). Surface Evolver also finds three local equilibrium shapes closely corresponding to these patterns (Fig. 4 C). The differences in energy between these configurations, as determined by Surface Evolver, are very small compared with the energy required to deform one configuration into another, so that their coexistence is indeed expected.
The ability of our physical model to explain both normal eye geometry and the eye geometries within different genetic mutants supports the underlying simple physical model for cell and tissue mechanics. Going beyond cell sorting or qualitative organization, all observed quantitative features of ommatidium geometry have been faithfully reproduced with a physical model that requires very few free parameters. The adhesion strengths of E and Ncadherin predicted by the model are close to the measured contributions of cadherins to surface tension in vivo and play a critical role in determining cell geometry. Bulk elastic energy terms (those allowing for volume shifts in the cell body and changes of crosssectional area in the adherens junction) proved unnecessary for an accurate model—only the interfacial energies of membrane area changes are needed. In the future, cell body elasticity may have to be included to model the morphogenetic development of the retina, which involves more severe deformations and substantial crosssectional area changes. Overall, our findings show that—at least for this system—interfacial energy functionals as influenced by cell adhesion molecules, are a simple, quantitative, and efficient way to understand and predict the shapes of adhering cells in tissues.
Materials and Methods
The experimental images used for this work consisted of confocal fluorescence micrographs (24) and cobalt sulfidestained micrographs (10, 24–29). In all, 60 different wildtype ommatidia were analyzed. By using image analysis algorithms built in or modified from Matlab 7.1 (21), the images were (i) converted to grayscale, (ii) thresholded, and (iii) skeletonized to convert membrane images to singlepixel width. We developed separate Matlab routines to extract the positions of vertices (meeting points of three membrane edges) from the skeleta, and to perform a bestfit of arcs of circles to each of the edges by using a goldenmean bisection algorithm (30). With the vertex positions and circular arcs given, it is easy to extract membrane edge lengths and angles between edges by planar analytical geometry. The entire image analysis for all sample ommatidia takes a few hours on a PC.
Surface Evolver simulations were performed with Surface Evolver 2.26c (22) under both Linux and Windows. Results were checked for consistency between the systems and for insensitivity to initial conditions or the details of the grid refinements used by the Surface Evolver program. The interface energy functional (Eq. 1 ) was programmed into Surface Evolver by using method instances keeping track of the total length of each cell membrane, and the cadherin adhesion strengths were treated as constant (negative) surface tensions. The evaluation of the penalty function (Eq. 7 ) and the search for its minimum were performed with Mathematica (31).
The system of Eqs. 2 – 5 was analyzed as a linear equation with the matrix and the vector of variables As a system of m = 4 equations with n = 5 variables, the solution space is expected to be onedimensional for a generic matrix of rank r = m. Indeed, one easily finds that the only solutions of Eq. 8 for arbitrary A _{i} are the multiples of (1/2, 1/2, 1/2, 1, 0), that is, these solutions have γ_{N} = 0. Because it is known that Ncadherin plays a role in determining ommatidium geometry, these solutions are not physical. Consequently, the matrix must have lower rank r < m. A standard singular value decomposition of M (30) shows that rank r = m − 1 requires the condition (Eq. 6 ) above. The resulting twodimensional solution space can then be written as where the coefficients are functions of the angle cosines, namely, This establishes an analytical relation between the relative membrane extensions Δ_{i} and the angles between the membranes. The solution shapes obtained by Surface Evolver fulfill all Eqs. 8 – 14 within numerical accuracy, establishing that, at given γ_{E}, γ_{N} > 0, the energy minimization process indeed finds this solution of rank m − 1. For the experimentally measured values of A _{1} through A _{3} and the resultant A _{4}, the values are f _{P} ≈ 0.58, f _{C1} ≈ 1.39, and f _{C2} ≈ 1.21, showing that the contributions of both solution base vectors in Eq. 11 are of a comparable order of magnitude. We attempted to include curvature (bending) energy terms in the energy functional, but these terms did not contribute significantly to the overall energy or to the final equilibrium shape.
Acknowledgments
We thank Francois Graner and Jos Kaefer for their insight and advice on alternative simulation methods, Takashi Hayashi for discussions and provision of image material, Jon Widom for reading the manuscript, and Kenneth Brakke for crucial help with Surface Evolver. This work was supported by the National Institutes of Health (R.W.C.), a 3M Nontenured Faculty Grant (S.H.), and a National Science Foundation Research Training Grant in Applied Mathematics (S.E.).
Footnotes
 ^{†}To whom correspondence may be addressed. Email: sascha{at}northwestern.edu or rcarthew{at}northwestern.edu

Author contributions: S.H. and R.W.C. designed research; S.H. and S.E. contributed new reagents/analytic tools; S.H., S.E., and R.W.C. analyzed data; and S.H. and R.W.C. wrote the paper.

The authors declare no conflict of interest.

This article contains supporting information online at www.pnas.org/cgi/content/full/0711077105/DC1.
 © 2008 by The National Academy of Sciences of the USA
References
 ↵
 ↵

↵
 Tilney LG ,
 Saunders JC

↵
 Steinberg MS
 ↵

↵
 Brodland GW

↵
 Mizuno D ,
 Tardin C ,
 Schmidt CF ,
 MacKintosh FC

↵
 Wolff T ,
 Ready D
 Bate M ,
 MartinezArias A

↵
 Strausfeld NJ ,
 Nassel DR
 Autrum H
 ↵

↵
 Nern A ,
 et al.

↵
 Weaire DL ,
 Hutzler S

↵
 Plateau JAF

↵
 Käfer J ,
 et al.
 ↵

↵
 Safran SA
 ↵

↵
 Boal DH
 ↵

↵
 Landau LD ,
 Lifshitz EM

↵
 Biran AB ,
 Breiner MMG
 ↵
 ↵
 ↵

↵
 Crew JR ,
 Batterham P ,
 Pollock JA

↵
 Takatsu Y ,
 et al.

↵
 Freeman M
 ↵
 ↵

↵
 Press WH ,
 Flannery BP ,
 Teukolsky SA ,
 Vetterling WT

↵
 Wolfram S
References
 ↵
 ↵

↵
 Tilney LG ,
 Saunders JC

↵
 Steinberg MS
 ↵

↵
 Brodland GW

↵
 Mizuno D ,
 Tardin C ,
 Schmidt CF ,
 MacKintosh FC

↵
 Wolff T ,
 Ready D
 Bate M ,
 MartinezArias A

↵
 Strausfeld NJ ,
 Nassel DR
 Autrum H
 ↵

↵
 Nern A ,
 et al.

↵
 Weaire DL ,
 Hutzler S

↵
 Plateau JAF

↵
 Käfer J ,
 et al.
 ↵

↵
 Safran SA
 ↵

↵
 Boal DH
 ↵

↵
 Landau LD ,
 Lifshitz EM

↵
 Biran AB ,
 Breiner MMG
 ↵
 ↵
 ↵

↵
 Crew JR ,
 Batterham P ,
 Pollock JA

↵
 Takatsu Y ,
 et al.

↵
 Freeman M
 ↵
 ↵

↵
 Press WH ,
 Flannery BP ,
 Teukolsky SA ,
 Vetterling WT

↵
 Wolfram S