Abstract
As animals develop, tissue bending contributes to shape the organs into complex threedimensional structures. However, the architecture and packing of curved epithelia remains largely unknown. Here we show by means of mathematical modelling that cells in bent epithelia can undergo intercalations along the apicobasal axis. This phenomenon forces cells to have different neighbours in their basal and apical surfaces. As a consequence, epithelial cells adopt a novel shape that we term “scutoid”. The detailed analysis of diverse tissues confirms that generation of apicobasal intercalations between cells is a common feature during morphogenesis. Using biophysical arguments, we propose that scutoids make possible the minimization of the tissue energy and stabilize threedimensional packing. Hence, we conclude that scutoids are one of nature's solutions to achieve epithelial bending. Our findings pave the way to understand the threedimensional organization of epithelial organs.
Introduction
Epithelial cells are the building blocks of metazoa. These bricks display columnar, cubic, or squamous shapes and organize in simple or multilayer arrangements. Faithful execution of the body plan during morphogenesis requires a complex reshaping of epithelial tissues to achieve organ development. In this context, the transition from planar epithelial sheets to cylindrical, ellipsoidal, or spherical forms, involves fundamental reorganization of the cells along their apicobasal axes. The coordination of these individual cell shape changes has been shown to induce large tissue rearrangements^{1,2,3,4,5}.
As for tissue cellular organization, the apical surface of cells has been assumed to be a faithful proxy for their threedimensional (3D) shape. Consequently, epithelial cells have been depicted as prisms with polygonal apical and basal faces. For example, during tissue invagination processes, such as the Drosophila mesoderm furrowing or vertebrate neurulation, epithelial cells change their shape from columnar to the socalled “bottle” form^{1,2,3}. When schematized, the bottle shape is pictured as a variation of a prism, the frustum, i.e., the portion of a pyramid that remains between two parallel planes^{6}. Frusta display apical and basal polygonal faces with the same number of sides but with a different area^{1,2,3}. Thus, it is generally assumed that the cell organization in the apical surface drives the epithelial 3D architecture.
The arrangement of cells in the apical surface of the epithelium has been extensively analysed from biophysical, mechanical, and topological viewpoints^{1,7,8,9,10,11,12,13,14,15,16}. These studies have been essential to understand fundamental morphogenetic processes, such as convergent extension, tissue size and shape control, and organogenesis. Topologically, the apical surface of epithelial sheets is arranged similarly to Voronoi diagrams. The Voronoi formalism has been shown to be useful to understand the mechanisms underlying tissue organization in the plane of the epithelium^{7,17}. Moreover, any curved surface and 3D structure can be partitioned by means of Voronoi cells using computational geometry tools^{18,19,20}.
Several groups have tried to go beyond the twodimensional description of tissues combining computational models and experimental systems^{21,22}. This has been done by analysing the apical surface of 3D structures^{23,24} or by developing lateral vertex models to study epithelial invaginations^{25,26}. Recently, studies have focused on understanding 3D curved epithelia^{27,28}. Khan et al. quantified epithelial folding by tracking individual cells during Drosophila gastrulation and showed intercalations in the plane of the epithelium and shape changes^{29}. Other studies have addressed the emergence of curved 3D structures (e.g., tubes and spheroids) by means of numerical simulations^{3,21,22,30,31,32,33,34,35,36,37,38}. Notably, in all these works epithelial cells are, anew, described and modelled as either prisms or frusta. However, there is evidence that epithelial cells are able to contact different neighbouring cells at different depths along the apicobasal axis of the cell (contrary to the prism/frustum paradigm). The appearance of these intercalations along the apicalbasal axis has been observed in the columnar epithelium of Drosophila imaginal discs^{39} or during Drosophila germband extension^{40}^{,}^{41} and has been also modelled computationally in the context of a planar tissue^{42}.
Altogether, there is a gap of knowledge about the 3D packing of epithelial cells in curved tissues and, by extension, about the associated morphogenetic processes that create these structures. In addition to this fundamental aspect of morphogenesis, the ability to engineer tissues and organs in future critically relies on the ability to understand, and then control, the 3D organization of cells^{43,44}. Here we combine a theoretical/computational framework with experimental data to quantify and characterize the 3D structure of curved epithelia. Remarkably, our approach unveils new aspects of cells and their packing properties. Thus, in curved tissues such as tubes, vaults, or spheroids, epithelial cells adopt a previously undescribed geometrical shape that makes their packing energetically efficient. Here we propose that such geometrical conformation is one of nature's solutions to epithelial bending.
Results
A tubular model reveals apicobasal cell intercalations
To investigate whether 3D packing in curved epithelia could be simply explained by prismatic cells adopting a frustalike shape (Fig. 1a, b), we computed a cylindrical epithelium that mimics the architecture of epithelial tubes and glands (Methods section, Fig. 1c). In this model, the ratio between the outer cylinder radius (R_{b}, basal radius) and the inner cylinder radius (R_{a}, apical radius) determines how large the relative expansion of the basal surface is with respect to the apical surface. We called this magnitude the surface ratio of the tube: R_{b}/R_{a} (Fig. 1c). The model assumes that apical and basal surfaces behave as Voronoi diagrams. A Voronoi diagram is generated by a set of seeds in space. A Voronoi cell is then defined by the unique region of space that contains all points closer to a given seed than to any other. First, we generated a Voronoi diagram on the inner surface (apical). Then, the seed of each Voronoi cell was projected to the closest point on the surface of an outer cylinder (basal) generating a second Voronoi diagram in the basal surface of the tube. Interestingly, the Voronoi diagrams that emerged from apical and basal surfaces showed different topologies: their cells had different neighbours on each surface (Fig. 1c). This change was due to the asymmetric enlargement of the area between basal and apical surfaces that caused the redistribution of the distances between the corresponding seeds. The distance between two neighbouring seeds along the longitudinal direction of the tube was the same in apical and basal surfaces (Fig. 1c, green and red cells). On the contrary, the distance between neighbouring seeds along the transverse direction increased on the basal surface (Fig. 1c, blue and yellow cells).
This topological argument entails a fundamental consequence: a change in the identity of neighbouring cells from the apical to the basal surface is not compatible with the view of epithelial cells as either prisms or frusta (Fig. 1a, b). In fact, in a fourcell motif (Fig. 1c), the neighbour exchange means a cell intercalation process along the apicobasal axis. This spatial transition would be analogous to the temporal T1 transition leading to cell intercalations in a number of morphogenetic events, such as convergent extension^{14,45,46}. Therefore, in this manuscript, we will refer to this topological change in fourcell motifs as apicobasal intercalations or apicobasal transitions.
We then analysed in detail the theoretical geometrical shapes or solids participating in this topological transition. The exchange of neighbours involved the appearance of a vertex along the apicobasal axis of these solids, at the location where the transition takes place (Fig. 1d). To the best of our knowledge, this type of geometric figures had been previously undescribed. We coin the term scutoid to name this geometrical solid for its resemblance with the scutum/scutellum shapes in the thorax of some insects such as the Cetoniidae subfamily of beetles (Fig. 1e). Scutoids, in contrast to prisms/frusta, are not necessarily convex, and present nonplanar lateral faces (Fig. 1d, f, g, and Methods). This is due to the fact that the boundary between cells follow geodesic trajectories (the shortest path connecting two points)^{20}. Importantly, this ensures that scutoids display concave surfaces that allow the 3D packing of cells in a curved tissue (Fig. 1d, f, g).
Scutoids appear in tubular epithelia
Our model predicts that scutoids should exist in actual curved epithelia. To determine this, we examined the third larval instar salivary gland of Drosophila, a model for tubulogenesis that is extensively used, where cubic secretory cells arrange around a lumen to form a cylinder^{47,48} (Fig. 2a). We developed a method to identify and track in depth (i.e. along the apicobasal axis) the epithelial cells in this gland (Fig. 2a, Supplementary Movie 1, Methods). In agreement with our predictions, we identified numerous apicobasal transitions (716 cells examined, 75 ± 17% of scutoids, surface ratio: 7 ± 3, Methods and Fig. 2b). Using confocal images, the 3D reconstruction of cell volumes confirmed the appearance of “scutoidal” shapes (Figs. 1g, 2c).
To assess additional geometrical features underlying the emergence of apicobasal transitions, we exploited further the predictive capabilities of the Voronoi tubular model. (Fig. 3a, Supplementary Fig. 1 and Methods). First, we computed nine exterior expansions of the tube to span a wide range of surface ratio conditions. We obtained a series of increasing surface ratios with the formula \(1{/(1  x)}\) taking x the values {0.1, 0.2, …, 0.9} (Fig. 3a, b and Methods). The expansion of the exterior face of the tube caused the basal Voronoi cells to increase proportionally their area, to tile the surface (Fig. 3a). As a baseline, we obtained the model for the planar epithelium \(\left(\frac{R_{b}}{R_{a}} = 1 \right)\) that represented the case without a topological change between the two surfaces (Supplementary Fig. 1). At the same time, we examined the influence of cell density by testing a different number of cells (40, 80, 200, 400, and 800 cells) on the same surface (Fig. 3a, b and Methods). We quantified the percentage of cells presenting at least one apicobasal transition (Fig. 3b and Supplementary Data 1). Such percentage increased as a function of the basal surface enlargement in all the cases, reaching 100% when \(\frac{R_{b}}{R_{a}} = 10\) (Fig. 3b). The increment rate was similar for the models with 80, 200, 400, and 800 cells, and slightly smaller for 40 cells. In addition, we measured the number of transitions per cell. The expansion of the basal surface produced an increasing number of cells participating in several intercalations. Again, the numbers were lower in the case of 40 cells (Supplementary Fig. 2 and Supplementary Data 1). These results confirmed that the change of surface ratio, leading to variations of the cellular aspect ratio, is the main contributor to the emergence of apicobasal transitions and that cell density affects the appearance of scutoids in a minor way.
Since the apicobasal intercalations necessarily imply changes in the edge shared by cells, we also investigated the geometrical traits of fourcell motifs. Note that the aforementioned edge has to change its size and rotate during an apicobasal transition (Fig. 3c and Methods). To that end, we first collected information from fourcells motifs, with and without apicobasal transitions (n = 40 larval salivary glands, 146 motifs with transitions, 656 motifs without transitions, Methods). The measurements were taken using the confocal images of the outer (basal) surface of the glands. We observed that the angle of the edge between the contacting cells was smaller than 30º in most of the transitions analysed; in contrast, the no transitions cases showed a much wider range of angles (Fig. 3d). In addition, we found that the length of the edge was shorter in fourcell motifs presenting apicobasal intercalations (Fig. 3d and Supplementary Data 2). As expected, we obtained equivalent results in the corresponding Voronoi tubes that show a similar percentage of scutoids (1/0.6 surface ratio) (Fig. 3e). In this regard, the analysis of the different Voronoi tubes revealed two clear results: (i) for small surface ratios the transitions appeared mainly when the edge angle was between 0º and 15º with respect to the transverse axis of the cylinder (measured in the basal surface); (ii) an increase of the surface ratio led the motifs with transitions to have larger angles and edge lengths (Supplementary Fig. 2). In all these cases, there were significant statistical differences between transitions and no transitions regarding edge angle and edge length distributions (Kolmogorov–Smirnov test; Supplementary Data 1) indicating that in tubular geometries, the appearance of apicobasal transitions is not random but associated to certain geometrical traits of the fourcell motifs.
A linetension minimization framework to analyse 3D packing
The geometric characteristics of the solids that appear due to the apicobasal cell transitions suggest the existence of a minimization principle governing the cellular packing motifs. From that perspective, we extended the theoretical model to analyse the biophysical origin and consequences of the appearance of the apicobasal transitions (Fig. 4a, b). Given that the latter is ultimately spatial T1 transitions, our approach focuses on the minimal model able to generate neighbour exchanges. In that regard, it has been shown that linetension energy variations, i.e., the energetic changes due to the remodelling of the membrane shared by cells in a given apicobasal plane, are necessary and sufficient to induce T1 transitions^{49} (Discussion).
Thus, a cell packing configuration in any plane perpendicular to the basalapical direction is circumscribed to a quadrilateral (Fig. 4a). Our theoretical model relies on an idealized situation where packing is described in terms of two limiting stable configurations that are characterized by the length of their edges, l_{w} and l_{h}. Here, for simplicity, we also use l_{w} and l_{h} to refer to the packing configurations. Therefore, when estimating the energy from experimental measurements (Methods), we assumed that an actual configuration can be effectively represented as a linear combination of those fundamental modes (Fig. 4a). In addition, a third, unstable, configuration is possible if l_{w} = l_{h} = 0, (fourfold vertex configuration, Fig. 4a, b). By focusing on the tensile forces controlling cell interactions, an energy functional can be defined^{35,46,50,51,52}. Under these conditions, the energy landscape of a particular motif can be determined, and the stability of packing configurations can then be ascertained as a function of the length of the edge and the cellular aspect ratio, see Fig. 4b (Methods). We note that the changes in the aspect ratio, \({\it{\epsilon }} = \frac{h}{w}\), as a readout of the surface ratio anisotropy, relies on the assumption that cellular motifs are not affected by neighbouring configurations.
According to our Voronoi model, the larger the surface ratio of a tube, the more scutoids develop (Fig. 3b). Also, our data indicate that the shorter the edge length, the more likely it is that a transition appears (Fig. 3d, e). These two phenomena can be qualitatively explained by our linetension minimization framework since, on the one hand, the energy landscape is further explored in terms of the aspect ratio (surface ratio anisotropy) and, on the other hand, for a given tube, scutoids (\(l_w \leftrightarrow l_h\) transitions) are favoured when, close to the unstable fourfold vertex configuration, the surface ratio anisotropy enables jumping above the energy barrier (Fig. 4b and Methods). The experimental edge configurations and the values spanned by the aspect ratio of cellular motifs leading to no transitions or transitions in Voronoi tubes and salivary glands are in agreement with our energetic considerations about these data (Fig. 4c, Supplementary Fig. 3, and Methods). This suggests an analysis of transitions and no transitions from the viewpoint of the path followed in the energy landscape to explore the attractors and enable energy relaxation. That path can be effectively characterized by basalapical trajectories in the E_{h}, E_{w} energy space (Methods), thus accounting for the change of the weight of each fundamental packing configuration as one moves along the apicobasal axis of the cell. We computed the energy of cell configurations in the apical and basal surfaces in salivary glands and Voronoi tubes (Fig. 4d, e, Methods). In the case of no transitions (frustra packing), we obtained short, disordered, local rearrangements corresponding to fluctuations around the same energetic attractor. On the other hand, in the case of transitions, we found long trajectories due to jumps over the fourfold vertex energy barrier to reach a new attractor and achieve further stability (Fig. 4b, d, e). Interestingly, these trajectories showed a clear directionality as one moves from the basal to the apical surface of the tissue (from l_{h} to l_{w}) in agreement with what is expected from the theoretical energy landscape.
Scutoids are a general feature of curved epithelia
Our model predicts that apicobasal cell intercalations may appear whenever there is an asymmetric enlargement of the basal and apical areas. Thus, we predicted that the appearance of scutoids in curved epithelia should not to be an exclusive of the larval salivary gland but rather, a general cornerstone of epithelial architecture. In this regard, we analysed the epithelial folds in the Drosophila developing embryo (4 embryos analysed, 214 cells examined, 50 ± 15% of scutoids, surface ratio: 1.6 ± 0.2, Methods and Supplementary Fig. 4a–c). We also selected, for the same embryos, 5 regions of interest at locations without folds and a smaller surface ratio (1.2 ± 0.1). In this case, the number of scutoids decreased (734 cells examined, 15 ± 4% of scutoids, Methods). These results confirmed, in an independent context, a correlation between surface ratio anisotropy and the appearance of scutoids in curved epithelia.
A tubular geometry, or a fold, is a particular case where one of the main curvatures is null, whereas the other is constant^{53}. To extend our conclusions to a general framework, we explored other curved epithelia that have nonnull curvatures along their symmetry axes (closer to a sphere or a spheroid). The Drosophila egg chamber is a spheroidal developmental structure, where a monolayer of follicle cells surrounds a group of germ cells. During development, the egg chamber elongates from a quasispherical shape to an apparent spheroidal form. This dynamic process is orchestrated through a collective migration of the follicular cells^{54,55,56}. We obtained and processed confocal stacks from stage 4 (mild elongation) and stage 8 (complete elongation) egg chambers (Fig. 5a, b, Methods) and the data revealed the presence of apicobasal transitions in the follicular epithelium at both stages (Fig. 5c and Methods). In addition, we found that the percentage of scutoids decreased during the egg chamber maturation process (stage 4, 13 egg chambers analysed, 719 cells examined, 20 ± 8% scutoids; stage 8, 14 egg chambers analysed, 1283 cells examined, 10 ± 3% scutoids, Methods). This raises interesting questions about the importance of nonstationary processes during development (Discussion).
Contrary to the larval salivary glands, but likewise the locations without folds of the Drosophila embryo (Supplementary Fig. 4a–c, f and Supplementary Data 2), we did not detect a preferential orientation of the edge angle in fourcells motifs leading to apicobasal transitions in egg chambers (Fig. 5d, Supplementary Data 2 and Methods). This was also the case in the quasispherical multilayered Zebrafish epithelium at 50% epiboly stage, where while we registered apicobasal cell intercalations (7 embryos analysed, 6724 cells examined, 2.9 ± 1.5% of scutoids, Supplementary Fig. 4d, e, Supplementary Data 2 and Methods), the orientation of the edge of fourcells motifs did not display a clear bias (Supplementary Fig. 4g and Supplementary Data 2).
A Voronoi model for spheroidal epithelia
The aforementioned experimental evidence suggests that scutoids are a general feature of many curved epithelia. In order to further investigate the 3D packing in tissues that have nonnull curvatures along both symmetry axes, we implemented a Voronoi spheroidal model (Methods). We note that in the particular case of the sphere, the basal to apical surface reduction is isotropic. Consequently, there is no change of the aspect ratio and the cells are expected to have the frusta shape^{53}.
First, we tested the appearance of scutoids under different combinations of the length of the spheroid radii and cell apicobasal sizes (Supplementary Fig. 5a–c and Supplementary Table 2). We observed that the increase in the percentage of scutoids correlated positively with the cell apicobasal length and with the ratio of the length of the axes. However, in the case of the sphere (i.e., when the length of the symmetry axes matched), we did not observe any scutoid regardless of the apicobasal length of the cell (Supplementary Fig. 5a and Supplementary Table 2). These data support our energetic model, where apicobasal transitions are driven by the anisotropy of the main curvatures ratios between both surfaces, which ultimately leads to changes in the aspect ratio.
We computationally mimicked stage 4 and 8 egg chambers by using the measured average values of the symmetry axes lengths and the cell apicobasal size (Fig. 5e, f, Supplementary Data 3 and Methods). In the case of the Voronoi spheroidal model stage 4, the percentage of scutoids registered was smaller than in its actual tissue (130/15,285 simulations/cells analysed, 6 ± 4% scutoids). This percentage increased in the Voronoi spheroidal model stage 8 (30/7051 simulations/cells analysed, 11 ± 4% scutoids, Methods). It is worth to notice that in the actual egg chamber the number of scutoids decreased as this developmental structure became more prolate (i.e., a spheroid elongated towards the poles) from stage 4 to stage 8 (see Discussion). The second difference with respect to the actual egg chambers was the appearance of the bias, due to the differential curvature, in the orientation of the “transitions” as in tubular geometries (see Discussion, Supplementary Fig. 5d, e and Supplementary Table 1).
Finally, we computed the edge configurations and the values spanned by the aspect ratio of cellular motifs and also the energetic trajectories in the case of the transitions and no transitions for both the Voronoi spheroidal model and the egg chambers (Fig. 5g–h and Supplementary Fig. 3c–f). The results of our analyses were compatible with those in the case of tubular geometries, including the long trajectories in the case of transition packing and the local fluctuating character of no transition packing. Interestingly, the directionality of the trajectories in the case of transitions in the actual egg chambers suggested a role of active developmental processes in the 3D tissue packing (see Discussion). Altogether, our results support the hypothesis that the surface ratio anisotropy is able to drive energetic transitions to provide tissue stability, and therefore, scutoids facilitate 3D packing of curved epithelia.
Discussion
Throughout embryogenesis, epithelia undergo dramatic morphogenetic changes that involve tissue bending and invagination. In general, the formation of any curved epithelium requires a fundamental 3D tiling adjustment: the same cells have to pave two surfaces that could vary largely in area^{57}. Until now, this was envisioned to occur by a reduction on the surface area of the cells in just one side: i.e., an apical (or basal) constriction that confers a bottle shape to the cells^{1,2,3,58} (Fig. 1a, b). Our results rule out this possibility as a general explanation of tissue packing. Instead, we propose a model in which curved epithelia are not formed just by prisms or frusta when a certain surface ratio asymmetry is exceeded (Fig. 6a, b). This prediction is supported by the identification of cell intercalations along the cellular apicobasal axis in diverse curved epithelia (Figs. 2 and 5, Supplementary Fig. 4). Therefore, our results reveal that curved tissue packing relies on another solution that manages cell shape and packing: apicobasal intercalations that give rise to scutoids. These previously undescribed geometric figures are not necessarily convex solids and allow the accommodation between concave and convex surfaces in the apicobasal transitions. In our model, we typically observe that contacts between the cells can be concave (Fig. 1g). This feature has been observed in diverse epithelia^{59} and is very clear on the basal surface of the larval salivary gland (Fig. 2b, c).
Our Voronoi tubular model predicts and explains the appearance of the apicobasal transitions and key topological features. We have found that the increment of the surface ratio positively correlates with the percentage of scutoids (Fig. 3b). This behaviour was also observed in the Drosophila larval salivary glands and in the Drosophila embryo folds. The agreement between the results obtained in the Voronoi tubular models and in these tissues also includes the preferential orientation of the motifs that are involved in an apicobasal transition (Fig. 3d, e). Yet, we acknowledge limitations to explain all experimental data by our modelling approach. For a given surface ratio, a larger number of scutoids develops in the Voronoi tubular model compared to the case of the larval salivary gland. We argue that these differences could be explained by the biophysical characteristics of the cells, such as their stiffness or surface stress, that are neglected in the Voronoi model. In practice, cellular reorganization must actively overcome energy barriers whereas a purely geometrical description assumes an idealized (frictionfree) scenario.
The comparison between the egg chambers and the spheroid models brings very interesting conclusions. In the case of the stage 4 egg chamber, the percentage of scutoids is larger in the epithelium than in the Voronoi model. Moreover, in stage 4 follicular epithelium, we do not find a relationship between the orientation of the fourcells motifs and the appearance of apicobasal transitions (Fig. 5d). From stage 1 to 8, the follicular epithelium of the egg chamber is a very dynamic tissue that is rotating due to a coordinated migration that induces the elongation of the whole structure^{54}. These morphogenetic events are accompanied by an increase in size mediated by cell proliferation that ceases at stage 6. Consequently, in this context, our analyses were performed under highly dynamic conditions in which the tissues had not yet adopted their final morphology (contrary to the stationary condition in the larval salivary gland). Interestingly, it has been recently shown how epithelial cells can form migrating basolateral protrusions during cell intercalation in the Drosophila germband extension process^{40}. These protrusions would initiate the tissue rearrangement, and imply, like scutoids, a different topology in the apical and basal surfaces. Thus, we propose that the observed differences are a consequence of the active cell and developmental processes. This hypothesis is supported by the double directionality of the energetic trajectories in the actual egg chambers (compare insets in Fig. 5g, h). In contrast, in the case of Voronoi spheroids, the trajectories are clearly driven by energetic cues following the most stable state as a function of the surface ratio anisotropy. In general, we reckon that 3D epithelial packing is a multifactorial process and that the plasticity of the epithelial cells in some conditions (e.g., long columnar cells) will be coordinated with the geometrical constraints to accommodate cells in the tissue. All these findings open new promising directions, yet outside of the scope of this study.
Tissue bending implies the difference in size between the apical and basal surfaces of epithelial cells, which in turn induces the emergence of apicobasal intercalations. As we mentioned above, apicobasal cell intercalations are analogous to the wellknown T1 transitions^{14,45,46} where the role played by time in the latter is played by space (apicobasal radial coordinate) in the former (Fig. 2b). T1 transitions have been shown to enable cellular rearrangements that stabilize tissues^{60,61}. In this context, our tensile energetic model recapitulates previous results showing that changes in the cellular aspect ratio promote cellular rearrangement^{62}. Here we propose that scutoids help to stabilize 3D packing in curved epithelia. To explore this idea our model focuses on the minimal energetic considerations that can induce neighbour reorganization: tensile forces. As in the case of T1 transitions, turgor pressure, surface stress, and the contractility of the cortex (including the effect of the actomyosin ring at the apical junctional complex in epithelial cells) could certainly contribute to shaping the conditions under which a scutoid develop. Moreover, our energetic model considers isolated cellular environments (i.e., configurations) thus disregarding the effect of the surrounding cells. Notwithstanding the fact that experimentally characterizing these additional energetic contributions is challenging, our results support the premise that the coupling between geometrical changes due to bending and linetension energy variations is a fundamental driving force in scutoid formation.
Additionally, our model facilitates a general intuition to understand 3D packing in epithelia. In the simple case of a tube, this can be understood in terms of the surface ratio (the relation between the dimensionless radii \(\hat R_{a}\) and \(\hat R_{b}\)); (Fig. 6a). When \(\hat R_{b}{/}\hat R_{a} = 1\), apical and basal surfaces present the same area, the epithelium is flat and formed by prisms. If \(\hat R_{b}{\mathrm{/}}\hat R_{a} \approx 1\), i.e., similar but not equal, like in the case of a tiny tube formed by a squamous epithelium (topleft corner), the differences between apical and basal surfaces are small and the whole tissue is formed by frusta. Note that this sort of epithelial packing occurs as long as \(\hat R_{a}\) and \(\hat R_{b}\) increase simultaneously (Fig. 6a, topleft to bottomright diagonal). On the contrary, when \(\hat R_{a} \gg \hat R_{b}\) or viceversa (topright corner and bottomleft corner) our model predicts that the tissue will be exclusively formed by scutoids (Fig. 6a). In intermediate situations, our model predicts different proportions of frusta and scutoids (Fig. 6a).
In the case of geometries involving two nonnull principal curvatures (e.g. spheroids), the same arguments apply: the anisotropy of the surface ratio is able to drive energetic transitions (Fig. 6b). In spheres, where there is not any change in the cellular aspect ratio as in planar epithelia, the appearance of scutoids due to tensile forces is not expected. Also, in tissues with a spheroidal geometry (e.g., egg chamber), we would expect that the more prolate the shape of the tissue, the larger the percentage of scutoids. Finally, the model also considers tubes as a limiting case where one of the curvatures is null. In that situation, the effect of the surface ratio increase over the percentage of scutoids is more pronounced.
In this work, we provide a more accurate and detailed view of epithelial 3D packing. It now becomes important to investigate how this new framework affects distinct processes during development and disease, from the morphogenesis of organs to the environment at the onset of tumour formation^{1,11,50,63}. Altogether, our study paves the way to understand the biomechanics of morphogenesis in developing organisms and sheds light on the underlying logic of 3D cellular organization. This is fundamental not only for an understanding of tissue architecture during development and disease, but to the fields of tissue and organ engineering^{43}.
Methods
Geometry of Voronoi diagrams in curved surfaces
From the mathematical point of view, in order to model a tissue, we defined: the space, the metric, and, finally, the seeds. With these ingredients, we built a Voronoi diagram, whose cells will be a model for the epithelial cells in a curved tissue (Fig. 1c).
Regarding the space, we started with a given surface S, then for each point X (u.v) ∈ S, there exists a vector N(u,v) of length 1 that is orthogonal to the tangent plane in X(u,v) to S. Thus, for each λ∈ [0,1], it is possible to define a new surface S_{λ} parallel to S in such a way that any point of S_{λ} is X_{λ}(u,v) = X(u,v) + λN(u,v) (a point in one of the surfaces has an equivalent point in each one of the parallel surfaces).
The metric in each one of the surfaces previously defined is just the length of the shortest geodesic joining two points. In the case of the cylinder, the geodesics are helices in the cylinder. In the case of the sphere, the geodesics are arcs of greatcircles.
We defined every seed starting in a point on the apical surface. That point defined a segment between the basal and the apical surfaces by means of its normal (given the point X(u,v), the segment is X(u,v) + λN(u,v), λ∈ [0,1]). The intersection of these line segments with a given surface determined a seed. Thus, in order to generate all the seeds, in a first step we have chosen n points on the apical surface, then the n segments that were generated by them, and, finally, the intersection of those segments with every surface S_{λ} defined the seeds for that surface.
The next step was to compute the Voronoi diagrams of the seeds obtained in each one of the parallel surfaces. We linked the Voronoi regions corresponding to seeds on the same segment obtaining a 3D figure. In some cases, we obtained figures with vertices not located neither in apical nor in basal surfaces. These solids are not prisms or frusta (Fig. 1f, g). Moreover, it must be pointed out that this geometrical figure does not need to be convex.
Immunohistochemistry and confocal imaging of the epithelia
Flies were grown at 25 °C by employing standard culture techniques. Zebrafish were maintained and bred according to standard procedures^{64}. For each type of sample different immunohistochemistry conditions were used.
For Drosophila larval salivary glands we used the sqhGFP line^{65} and Cy3labeled phalloidin (Sigma) to label the cell contours of the epithelial cells. Salivary glands from third instar larvae were dissected in PBS and fixed with 4% paraformaldehyde in PBS for 20 min. The samples were washed three times for 10 min with PBT (PBS, 0.3% Triton) before phalloidin incubation for 1 h. Larval salivary glands were mounted using FluoromountG (Southern Biotech) using double sided tape to avoid the squashing of sample^{66}. Images were taken with a Nikon Eclipse TiE laser scanning confocal microscope. The images were captured using ×20 dry (for the complete glands that were segmented) and ×40 immersion objectives and exported as a 1024 × 1024 pixels TIFF file.
For Drosophila embryos, embryos collected for 6–8 h after egg laying were dechorionated and fixed in 1:1 formaldehyde 4% in PBS:nheptane for 20 min at room temperature and stained according to standard protocols. For basolateral membrane staining embryos were incubated at room temperature 4 h with antiDisc large antibody (1:200 in PBS—0.1% Tween—0.1% BSA; mouse, Developmental studies hybridome bank, 4F3). Secondary antibody was coupled to Alexa488 (Molecular Probes, 1:500). Images were taken on a SPE Leica confocal microscope equipped with a ×20 or ×40 ACSAPO oil objective and processed using FIJI and Adobe Photoshop CS6.
For Drosophila egg chambers, to label the follicle cells membrane the ubiquitously expressed membrane marker ResilleGFP (II) was used^{67}. For ovaries staining 1–2dayold females were fattened on yeast for 48 h before dissection. Ovaries were dissected in PBT (phosphatebuffered saline + 0.1% Tween 20). Fixation was performed incubating egg chambers with 4% paraformaldehyde in PBS on ice for 20 min. For nuclei labelling fixed ovaries were incubated with Hoechst (Molecular Probes ^{TM}) 1:1000 in PBT. The antibody goat antiGFPFITC (abcam, ab6662) was used 1:500. Individual ovarioles without the muscle sheath were mounted in Vectashield (Vector Laboratories). Egg chamber images were captured with a Leica TCSSPE confocal equipped with a ×40 (1.15 NA) ACSAPO oil objective.
For Zebrafish embryos membranes and nuclei of zebrafish embryos were labelled in vivo by microinjecting a bicistronic mRNA (GAP43GFP:h2bRPF) at onecell stage. mRNA for microinjection was synthesised using the mMessage Machine kit (Ambion), following manufacturer’s instructions. Embryos were grown up to pregastrula stage, fixed in 4% paraformaldehyde and embedded in 1.5% lowmelting point agarose (as in ref. ^{68} for imaging under an inverted confocal microscope “NIKON A1R+ in vivo”, using a 20X/0.75 dry objective. Images were processed using FIJI/ImageJ. All procedures comply with the guidelines from the European Community and Spanish legislation for the experimental use of animals.
Drosophila larval salivary gland measurements
In order to understand the cellular organization in a 3D salivary gland, we measured several geometrical properties from selected motifs. Importantly, to estimate these measurements we assumed that the glands have a cylindrical shape. These motifs were chosen considering the fourcells that were contacting alongside all the sections. Furthermore, at least two of these cells had to be central cells: a central cell was defined as a cell situated at the middle of the gland and its whole area was completely visible from basal to apical. As a rule of thumb, we counted as a valid motif when there were no fourfold vertex configurations in any of the two surfaces. We applied this consensus to all the images from actual tissues. Once the valid motifs were selected, we measured four properties related to epithelial architecture, directly in the confocal stack. We took as a reference the central edge that two cells of the motif were sharing in the basal (outer) surface (Fig. 3c). They were the surface ratio of the gland, the central edge length, the central edge angle of the motif, and the percentage of central cells with “scutoidal” shape
Surface ratio estimation: We measured the average width of the lumen (2R_{a}) and the average width of the salivary gland (2R_{b}) at the region of the selected motif by using the FIJI’s Straight Line tool. The surface ratio was then calculated as R_{b}/R_{a}.
Edge angle and length: We defined the angle of the motif’s intercellular edge, which is the edge shared by cells that have three neighbours within the motif, and the orientation it had with respect to the transverse axis. We quantified the edge angle and the edge length using FIJI’s angle and FIJI’s straightline tools respectively. (Fig. 3c, d).
Percentage of scutoids: We quantified the percentage of central cells participating in at least one transition along the apicobasal axis using all the sections of the confocal stack.
In addition, three geometric properties in the salivary glands images were measured to feed the linetension minimization model (Supplementary Data 3). They were the sum of all the edges length into the motif packing configuration (L_{T}), width (w) and height (h).
Total length of edges (L_{T}): We gauged and summed the length of every intercellular edge in the packing configuration from each fourcells motif (Fig. 4a, c). We used FIJI’s straightline tool.
Width (w) and height (h): We measured the “width” and “height” from each fourcells motif in basal and apical surfaces. The “width”/“height” is defined as the distance between two vertices in the quadrilateral that defines a fourcells configuration along the longitudinal/transverse axis (Fig. 4a, c). The average “width” \(\left\langle w \right\rangle\) and the average “height” \(\left\langle h \right\rangle\) of a motif were calculated as the average of each measurement: w_{1}and w_{2}, and h_{1} and h_{2} (Fig. 4a, c), which were estimated using FIJI straightline tool.
Drosophila embryo measurements
We analysed confocal stacks imaging the whole length of the Drosophila embryo epithelial during gastrulation stage. First, we selected two types of regions of interest (ROI), one corresponding to the flatter embryo surface, and a second one corresponding to the embryonic folds (Supplementary Fig. 4a–c). For each one of these ROIs, we counted the total number of cells and the percentage of scutoids using FIJI.
In the folds, we measured edge length and edge angle in fourcells motifs taking as a reference the basal surface. We chose ten fourcells motifs (at most) where "transition" occurred and another ten with "no transition”. The edge angle and length were measured similarly to the larval salivary gland (Supplementary Fig. 4f). The surface ratio was calculated by measuring the area of the fourcells motif in basal (using FIJI’s polygon selection tool) and dividing it by the area of the fourcells motif in apical. The same procedure was followed with the ROIs corresponding to the flatter surface of the embryo (Supplementary Fig. 4f). In this case, edge length and edge angle in fourcells motifs were measured on the apical surface. Then, we calculated the surface ratio by dividing the apical area by the basal area (since in this case, the epithelium presents a basal constriction).
Drosophila egg chamber measurements
To further investigate the presence of scutoids in Drosophila egg chamber, we have measured several geometric aspects of motifs with a transition and with no transition in both apical and basal layers. In particular, we calculated the value for the parameters w_{1}, w_{2},h_{1} and h_{2} and L_{T} (see linetension minimization model) along with edge angle and length (Fig. 4a, c), following the same procedure from the larval salivary gland. These motifs were captured from a visible central region of interest (ROI), discarding the cells from the borders of the egg chamber. We took all the follicular cells that appeared just before any nurse cell was completely visible. Moreover, we manually counted the number of cells that were found on the defined ROI and the total number of cells that appeared previous to the polar cells. We also estimated the number of scutoids (cells with at least a scutoidal side) and the cells which were not scutoids.
In order to estimate the geometrical shape of the egg chamber and its influence in the scutoids appearance, we measured the radii of the Xaxis, in anterior and posterior, and the Yaxis. Furthermore, to appraise the mean cell height of each sample, we measured the height of several cells that were located, laterally, on the region of the transverse axis, where the curvature was maximum.
Zebrafish embryo measurements
We inspected the first two layers of cells from zebrafish embryo samples at 50% epiboly. We first captured at most 10 cellular motifs with apicobasal transitions and 10 with no transition. To further analysed them, we segmented the fourcell motifs in their apical and basal frames, capturing the following measurements: edge angle and length regarding the Xaxis and surface ratio estimation (getting the ratio between the area of the motif in apical and in basal) (Supplementary Fig. d, e, g). We also counted the number of cells in the first two layers, along with the percentage of scutoids.
Voronoi tubular model
We have designed a geometric model to simulate curved epithelia. The model was based on generating a Voronoi diagram in a cylindrical shape. This process required that Voronoi cells located in an extreme of abscissa axis (the transverse axis of the cylinder image) were continuous in another abscissa extreme^{20} (Supplementary Fig. 1). We have used Matlab R2014b (Mathworks) as our computational tool. Taking as a reference the previous 2D Voronoi model established in^{7}, we implemented a Centroidal Voronoi Tessellation (CVT) on the lateral surface of the cylinder. To this end, we developed the cylinder lateral surface in the Euclidean plane and applied Lloyd's algorithm to a random Voronoi diagram^{69,70}.
These are the detailed steps for the method: (a) We randomly located a finite number of seeds (40, 80, 200, 400, and 800) in a defined Euclidean space (512 × 4096 pixels) for 20 different images. (b) We tripled the images with all the seeds in the direction of cylindrical transverse axis (abscissas axis), getting an image with size 1536 × 4096 pixels, and the number of seeds tripled (c) Voronoi regions were delimited by proximity over all the seeds, producing a Voronoi diagram. (d) The position of all the Voronoi cell centroids was calculated, and the centroids located between in the central 513 and 1024 abscissas coordinates were stored for next steps. (e) Images were cropped, capturing the [5131024] x [14096] pixels interval, giving as result 512 × 4096 pixels images. (f) The procedure from (a) to (e) was repeated using the calculated centroids as new seeds in (a) to perform Lloyd's algorithm. This algorithm was repeated four times for each image, achieving 5 diagrams for each initial random image.
To study the tubular model, we chose the Voronoi 5 diagram of the CVT, since it resembled a homogeneous distribution of the cells^{7}. The 20 Voronoi 5 diagrams that represented the cylinder were used to create a new lateral surface of a new cylinder with a larger radius and the same height (Fig. 3a and Supplementary Fig. 1). In this way, we were able to generate a model that represented the inner (initial) and the outer (expansion) lateral surfaces of a tube. The surface expansion was established in 9 steps using the following formula: \(f(x) = 1{\mathrm{/}}(1  x{\mathrm{/}}10),{\mathrm{I}} \subset {\Bbb N},{\mathrm{I}} = \left[ {1,9} \right],\quad \forall x \in {\mathrm{I}}:1 \le x \le 9\). The calculations were made by maintaining fixed the R_{a} length (inner surface) and increasing proportionally the R_{b} length (outer surface, Supplementary Fig. 1).
To partition the outer surface with Voronoi cells we used the following strategy: we multiplied the size of abscissas axis of the original image (apical) by the surface ratio (R_{b}/R_{a}) (Supplementary Fig. 1). In a similar way, abscissas coordinates (X) of the Voronoi seeds were multiplied by the same factor, getting a new position of the seeds in the basal surface that corresponds to the given surface ratio expansion. Using these seeds, we applied again (b), (c), and (e) steps, using the size of the enlarged image. The Voronoi cells corresponding to the original seed and their projections were tracked and labelled with a specific colour allowing their identification along different surfaces. (Figs. 1f–g and 3a and Supplementary Fig. 1).
Voronoi tubular model measurements
The Voronoi tubular model was used to measure the edge length and the edge angle of all the intercellular edges (Fig. 3c). These calculations were made on the basal surface and repeated for each one of the nine surface ratios. The angles were calculated taking as reference the transverse axis of the cylinder image (abscissas axis of a planar image). For each one of these intercellular edges, we established their involvement in apicobasal transitions (Fig. 3b, Supplementary Fig. 2a and Supplementary Data 1) and obtained the percentage of scutoids. We discarded the edges shared by nonvalid cells (cells located at the tube tip), and we only quantified scutoids along the set of valid cells (not located at the tube tip in any surface, apical or basal). Finally, 200 random measurements were taken from transition and no transition edges to represent the scatter polar graphics (angle vs edge length) (Fig. 3e and Supplementary Fig. 2b–e).
Additionally, a set of geometric properties were quantified to feed the linetension minimization model (Supplementary Data 3): the sum of all the edges length within the motif packing configuration (L_{T}), width (w) and height (h) (Fig. 4a, c).
To calculate the mentioned parameters, we captured the vertices involved in the packing configuration of fourcells motifs and we measured the distances between them. The fourcells motifs used in our analysis kept their fourcells together in apical and basal surfaces and we discarded the motifs that did not satisfy this requirement. We also disregarded the motifs with nonvalid cells on any of these surfaces and we classified our measurements: fourcells motifs with and without apicobasal transition. Similar to the process followed in the actual tissue images, and in order to avoid artefacts, we discarded the motifs presenting fourfold vertex configuration in one of the two surfaces. To make this process automatic, we only analysed the motifs with an edge length longer than 4 pixels. This number was selected after our analyses in the computed Voronoi spheres, where the basal to apical surface reduction is isotropic and all cells must have frusta shape^{53}. However, in the cases of fourcell motifs with very short edges, we found false positives (apicobasal transitions) in the sphere models. According to our analyses, 4 pixels set the minimum threshold for which false positives were removed. For consistency, we applied this condition to both tubular and spheroid models.
First, we got all the possible energetic measurements that satisfied out restrictions, discriminating between transition and no transition motifs. These datasets were used to represent the cellular motifs’ aspect ratio (Supplementary Fig. 3b). Second, we chose 100 samples of energetic measurements for both types of fourcells motifs: we randomly selected, at most, 10 measurements (5 from transition edges and 5 from no transition edges) from each cylindrical Voronoi Diagram (20 different diagrams). In summary, we took a maximum of 100 measurements in transition edges and 100 in no transition edges, for each surface ratio. This data set was used to map the energetic trajectories from basal to apical (Fig. 4e).
Voronoi spheroidal model
A 3D Voronoi model was developed simulating the behaviour of a curved epithelium with two nonnull principal curvatures. First, we implemented spheroidal models mimicking the Drosophila egg chamber in the stages 4 and 8 of development. Second, we explored several spheroidal structures (in terms of its radii and apicobasal length of the cells) in order to study the cell packing with different degrees of curvature. The inputs of the Voronoi spheroidal model were its three axial radii, number of cells (N) and apicobasal length of the cells.
We developed a pipeline to create different spheroidal shapes in terms of its axial radii, keeping constant Y and Z radii, while modifying X. By defining the axial radii Y and Z of the inner spheroid layer as the unitary value, we got the following spheroids: 10 random realizations of a sphere (X radius inner layer = 1), a balloonlike spheroid (X radius inner layer = 1.5) and a more prolate spheroid, ‘Zeppelin’shaped, (X radius inner layer = 2), with three different apicobasal length of the cells 0.5, 1 and 2 (Supplementary Fig. 5a–c); 130 random realizations of spheroid stage 4 (X radius inner layer = 1.26, cell height = 0.22, N = 200) and 30 random realizations of spheroid stage 8 (X radius inner layer = 2.16, cell height = 0.15, N = 450) (Fig. 5e, f), whose dimensions are proportional to the measurements taken in the Drosophila egg chamber at stage 4 and stage 8, respectively.
The pipeline followed these steps: (1) We put N seeds on the apical surface with a minimum distance between them (depending on the surface area), this distance avoided the overlapping and homogenized the seeds distribution along the ellipsoid. Note that N represents the maximum number of seeds because, in some cases, the surface could not be filled by all the seeds and satisfy the condition about the minimum distance between them simultaneously. (2) We extrapolated the seeds to the basal surface using the ellipsoid equation:
where parameters are: a, b, and c represent the radii along the X, Y, and Z axes of the ellipsoid respectively, and x, y, and z correspond to the X, Y, and Z coordinates in the ellipsoid surface, respectively. the case of the spheroid, b = c. (3) We determined the segments connecting the seeds in apical and their corresponding ones in basal. (4) We implemented the Voronoi algorithm, in the 3D space, taking as seeds the mentioned apicobasal segments to computed the Voronoi cells. (5) Finally, we applied again the ellipsoid equation, capturing only the surfaces of interest (apical and basal).
Voronoi spheroidal model measurements
The Voronoi spheroidal model was used to extract similar measurements that those taken in the Voronoi tubular model. Due to the impossibility to measure accurately in the 3D space, we took some measurements from the spheroid in an analogous way that was carried out for the confocal images of Drosophila egg chambers. This procedure is explained in detail below.
To transform each 3D spheroid into four 2D images, we computed two maximum intensity projections along the Zaxis (first from Z = zRadius to Z = 0, and second from Z = zRadius to Z = 0) and two others along the Yaxis (first from \(Y = {\mathrm{yRadius}}\) to Y = 0, and second from \(Y =  {\mathrm{yRadius}}\) to Y = 0) of each layer (outer and inner). After the projections were created, we homogenized the cell's border, eroding the cells and transforming their border to a width of one pixel. This procedure led to images similar to the ones obtained with the confocal microscope. Afterwards, we chose a region of interest (ROI) like the ones used in the Drosophila egg chambers measurements. This ROI excluded the cells from the border of the spheroid projections and only selected as valid cells, those in the central region (approximately from \(X =  \frac{2}{3}\ast {\mathrm{xRadius}}\) to \(X = \frac{2}{3}\ast {\mathrm{xRadius}}\)). We then selected all the fourcells motifs that composed the ROI and measured the same parameters analysed in the Voronoi tubular model: the intercellular edge angle (with respect to the transverse axis of the ellipsoid projection), the intercellular edge length, the percentage of scutoids (Supplementary Table 2), the sum of the edges length (L_{T}), the width (w) and the height (h) of the fourcells motifs (Fig. 4a, c). The energetic measurements followed exactly the same constraints that we applied to the Voronoi tubular model: we chose the cellular motifs present in both apical and basal layers, excluding the nonvalid motifs (containing nonvalid cells). We also distinguished between transition and no transition motifs. As explained above, we only selected fourcell motifs with the edge length longer than 4 pixels (see Voronoi tubular model measurements section). Furthermore, we extracted two sets of data: one with all the possible motifs used for plotting the cellular motifs’ aspect ratio; and another composed by 100 samples (selected randomly) to illustrate the energetic trajectories from basal to apical (Fig. 5h). To represent the geometric traits into polar scatter graphs, we randomly selected 200 edge length and angle measurements from transition edges and 200,200 from no transition edges (Supplementary Fig. 5d, e).
Additionally, we measured the maximum and minimum curvatures of the 3D ellipsoid coordinates in the inner and outer surfaces^{53}. We analysed the curvatures from \(X =  \frac{2}{3}\ast {\mathrm{xRadius}}\) to \(X = \frac{2}{3}\ast {\mathrm{xRadius}}\), covering the remaining axis entirely and matching the ROI of the projected images.
Kolmogorov–Smirnov statistical test
We performed a twosample Kolmogorov–Smirnov decision test to check if two samples came from the same continuous distribution (Supplementary Table 1). In particular, we analysed whether the edge angle and the edge length had different distributions in the motifs leading to transitions and to no transition. For this purpose, we have used the Matlab function ‘kstest2’, using as source data the same measurements represented in the polar scatters plots (Figs. 3d, e, 5d, Supplementary Fig. 2b–e, Supplementary Fig. 4f–g and Supplementary Fig. 5d–e).
3D segmentation of Drosophila larval salivary glands
We designed a protocol to segment, identify and track every epithelial cell in the confocal images taken from the larval salivary gland (Fig. 2a). In all the cases, the images were confocal stacks containing the cells contours along the height of the epithelial cells. The protocol had several steps: (1) Segmentation: We used Trainable Weka Segmentation^{71} (FIJI) a machine learningbased plugin. We labelled two categories: the cells’ outlines, and the inner part of the cells. The outlines included the perimeter of the cells and the lumen. The inner part included the body of the cells and the rest of the image. After the training step, we applied the plugin to all the planes of the confocal stack, getting a first set of segmented images. (2) Manual correction: this step was necessary to reflect rigorously the real cell outlines in the segmented images. This was especially important in the segmentation of the lumen of the larval salivary gland. The process was performed using Adobe Photoshop CS6. (3) Image homogenization: we developed a Matlab R2014b (Mathworks) script to execute and accomplish a final postprocessing to equalize the size of the cells’ border in order to improve the segmented images. (4) 3D segmentation through the confocal stack: we built a semiautomatic method that tracked the cells throughout plane \(Z\) of the confocal stack. This script, which was developed in Matlab, integrated the information from the positions and the areas of the cells at different stacks. Despite these calculations, it needed an additional step of human curation. As final output, we obtained segmented stacks where every cell was labelled (Fig. 2a and Supplementary Movie 1).
3D reconstructions
Using the data from the mathematical model, we performed the 3D reconstruction of a fourcell motif to visualize the surfaces of the cells along the apicobasal axis. We used the Voronoi cylinder images with different surface ratios as the source to generate the Voronoi tube. We applied the cylinder Matlab function to get a 3D representation of the whole model. We chose the original image that represented the inner part of the tube, a basal expansion image with surface ratio = 2.5, and we additionally built the intermediate layers between these surfaces. Cells were identified in each layer and we used Matlab’s alphaShape function to calculate the tightest fitting shape for each cell, producing a volume that embodied a 3D Voronoi cell (Fig. 1f, g). To represent the Voronoi tubular simulations with different surfaces ratios as a real 3D cylinder, we made use of Matlab’s function ‘wrap’ (Fig. 3a).
Using a similar approach, we obtained a 3D cell representation from motifs located at the Drosophila salivary gland. We selected a cellular motif from a segmented stack of images (see above). In order to model the 3D shape of each cell, first, we captured the boundary pixels of a cell in all the layers. Since the source images were obtained from plane confocal sections, we included a correction step where \(Z\) coordinates were extrapolated to a cylinder using an inferred cylinder radius for each layer. Once the cell boundaries were obtained, we defined the convex hull of the points which shaped the cell using boundary Matlab function. Finally, we reconstructed the 3D shape using the surface Matlab function (Fig. 2c). To illustrate the shape the spheroids that we analysed, and the contained cells, we performed a 3D reconstruction (Fig. 5e, f and Supplementary Fig. 5a–c). For this purpose, we have used the function ‘alphaShape’ to get the shape of each cell. The source data being the coordinates we got from applying the Voronoi algorithm to the initial seed segments thus representing cells’ height.
Linetension minimization model
The existence of tensile forces controlling cell interactions, either driven by cadherin and/or actomyosin activity, can be described through the following linetension energy functional:
where σ stands for the line tension, that we assume to be constant, and H for the Heaviside step function. Using w as a characteristic length scale and the energy of the unstable (fourfold vertex) configuration, \(E_0 = 2\sigma \sqrt {h^2 + w^2}\), as an energetic reference, the dimensionless energy of a fourcell configuration reads,
where \({\it{\epsilon }} \in (0,\infty )\) stands for the aspect ratio h/w and \(l \in (  1,{\it{\epsilon }})\), such that the negative and positive values represent l_{w} and l_{h}, respectively.
If dissipative forces are neglected, the relaxational dynamics of packing configurations due to tensile forces reads,
and stable configurations (force balance) are prescribed by the energy minima (Fig. 4b). As a function of the aspect ratio, different regimes can be distinguished: if \({\it{\epsilon }} \in \left( {0,\frac{1}{{\sqrt 3 }}} \right)\) there is a single global minimum located at l < 0, if \({\it{\epsilon }} \in \left( {\frac{1}{{\sqrt 3 }},\sqrt 3 } \right)\) there are two global minima located (one at l < 0 and another one at l > 0); finally, if \({\it{\epsilon }} \in (\sqrt 3 ,\infty )\) there is a single global minimum located at l > 0. Namely, if \({\it{\epsilon }} \in (0,\frac{1}{{\sqrt 3 }})\) then the only stable configuration for packing is l_{w} and if \({\it{\epsilon }} \in (\sqrt 3 ,\infty )\) the only stable packing configuration is l_{h}. In case \({\it{\epsilon }} \in \left( {\frac{1}{{\sqrt 3 }},\sqrt 3 } \right)\) both packing configurations are energetically stable.
On the one hand, no transition packing (i.e., frustra) is the only energetic stable configuration if along the radial coordinate between the apical and basal surfaces the range of variation of the aspect ratio is kept either smaller than \(\frac{1}{{\sqrt 3 }}\) or larger than \(\sqrt 3\). On the other hand, transition packing (i.e., scutoids) is the only possible stable configuration if along the radial coordinate there is a change in the aspect ratio from a value smaller than \(\frac{1}{{\sqrt 3 }}\) to a value larger than \(\sqrt 3\) (or viceversa). Under other conditions in terms of the variation of the aspect ratio (as a readout of the surface ratio anisotropy), we expect a mix of packing configurations since, as a function of the radial coordinate, a transition between configurations, \(l_w \leftrightarrow l_h\), can occur.
The experimental (in vivo or in silico) dimensionless values of the energy of a packing configuration was estimated by,
where \(\hat L_T\) is the total length of the packing configuration (all edges defining the motif inside a parallelogram) in units of \(\left\langle w \right\rangle\). The fundamental components of the energy, \(\hat E_w\) and \(\hat E_h\), are determined by a decomposition in terms of the director cosines,
such as \(\hat E_{{\rm exp}}^2 = \hat E_w^2 + \hat E_h^2\) and where θ is the measured angle formed by the transition edge. As for the value of the aspect ratio, we measured \(\left\langle w \right\rangle = (w_1 + w_2){\mathrm{/}}2\) and \(\left\langle h \right\rangle = (h_1 + h_2){\mathrm{/}}2\) (the average width and height of the configuration motifs, respectively) in the apical and basal surfaces (Fig. 4a).
In Supplementary Fig. 3, to characterize the experimental motifs in terms of l (either l_{h} or l_{w}) and \({\it{\epsilon }}\), we classify the edge configuration as a function of the value of the angle. If the latter is larger than 50 degrees we assumed an l_{w} configuration, if smaller than 40 we assumed an l_{h} configuration. The intermediate cases between 40 and 50 degrees were discarded in our analysis to avoid artefacts. The density plot was obtained estimating the probability density of couples, \(\left( {l,{\it{\epsilon }}} \right)\), with l normalized by the value of \(\left\langle w \right\rangle\), and applying a Gaussian kernel to smooth the results.
Code availability
The code is available at: https://github.com/ComplexOrganizationOfLivingMatter/Epithelia3D/tree/master/InSilicoModels/paperCode.
Data availability
The authors declare that all relevant data supporting the findings of this study are available within the paper and its Supplementary information files. Additional data are upon request.
Change history
08 October 2018
The original version of this Article contained an error in ref. 39, which incorrectly cited ‘Fristrom, D. & Fristrom, J. W. in The Development of Drosophila melanogaster (eds. Bate, M. & MartinezArias, A.) II, (Cold spring harbor laboratory press, 1993)’. The correct reference is ‘Condic, M.L, Fristrom, D. & Fristrom, J.W. Apical cell shape changes during Drosophila imaginal leg disc elongation: a novel morphogenetic mechanism. Development 111: 2333 (1991)’. Furthermore, the last sentence of the fourth paragraph of the introduction incorrectly omitted citation of work by Rupprecht et al. The correct citation is given below. These errors have now been corrected in both the PDF and HTML versions of the Article. Rupprecht, J.F., Ong, K.H., Yin, J., Huang, A., Dinh, H.H., Singh, A.P., Zhang, S., Yu, W. & Saunders, T.E. Geometric constraints alter cell arrangements within curved epithelial tissues. Mol. Biol. Cell 28, 35823594 (2017).
References
 1.
Lecuit, T. & Lenne, P. F. Cell surface mechanics and the control of cell shape, tissue patterns and morphogenesis. Nat. Rev. Mol. Cell. Biol. 8, 633–644 (2007).
 2.
Pearl, E. J., Li, J. & Green, J. B. A. Cellular systems for epithelial invagination. Philos. Trans. R. Soc. Lond. B. Biol. Sci. 372, pii: 20150526 (2017).
 3.
Davidson, L. A. Epithelial machines that shape the embryo. Trends Cell. Biol. 22, 82–87 (2012).
 4.
Escudero, L. M., Bischoff, M. & Freeman, M. Myosin II regulates complex cellular arrangement and epithelial architecture in Drosophila. Dev. Cell 13, 717–729 (2007).
 5.
Pilot, F. & Lecuit, T. Compartmentalized morphogenesis in epithelia: from cell to tissue shape. Dev. Dyn. 232, 685–694 (2005).
 6.
Schneider, P. J. & Eberly, D. H. Geometric tools for computer graphics. 10 Distance in 3D. Morgan Kaufmann publishers 397–398 (Boston, 2003).
 7.
SanchezGutierrez, D. et al. Fundamental physical cellular constraints drive selforganization of tissues. EMBO J. 35, 77–88 (2016).
 8.
Heller, D. et al. EpiTools: an opensource image analysis toolkit for quantifying epithelial growth dynamics. Dev. Cell 36, 103–116 (2016).
 9.
Mao, Y. et al. Differential proliferation rates generate patterns of mechanical tension that orient tissue growth. EMBO J. 32, 2790–2803 (2013).
 10.
Gibson, M. C., Patel, A. B., Nagpal, R. & Perrimon, N. The emergence of geometric order in proliferating metazoan epithelia. Nature 442, 1038–1041 (2006).
 11.
Farhadifar, R., Roper, J. C., Aigouy, B., Eaton, S. & Julicher, F. The influence of cell mechanics, cellcell interactions, and proliferation on epithelial packing. Curr. Biol. 17, 2095–2104 (2007).
 12.
Voronoi, G. F. Nouvelles applications des paramètres continus à la théorie de formes quadratiques. J. für die reine und Angew. Math. 134, 198–287 (1908).
 13.
Lewis, F. T. The correlation between cell division and the shapes and sizes of prismatic cells in the epidermis of cucumis. Anatom. Rec. 38, 341–376 (1928).
 14.
Zallen, J. A. & Zallen, R. Cellpattern disordering during convergent extension in Drosophila. J. Phys. Condens. Matter. 16, S5073–S5080 (2004).
 15.
Heisenberg, C. P. & Bellaiche, Y. Forces in tissue morphogenesis and patterning. Cell 153, 948–962 (2013).
 16.
CanelaXandri, O., Sagués, F., Casademunt, J. & Buceta, J. Dynamics and mechanical stability of the developing dorsoventral organizer of the wing imaginal disc. PLoS Comput. Biol. 7, e1002153 (2011).
 17.
Honda, H. Description of cellular patterns by Dirichlet domains: the twodimensional case. J. Theor. Biol. 72, 523–543 (1978).
 18.
Okabe, A., Boots, B., Sugihara, K. & Chiu, S. N. Spatial Tesselations. Concepts and Applications of Voronoi Diagrams (2009). https://doi.org/10.1002/0471721182.scard
 19.
Preparata, F. P. & Shamos, M. I. Computational geometry: an introduction. 6 Proximity: Variants and Generalizations SpringerVerlag New York. 47, 241–255 (1985).
 20.
Grima, C. I. & Márquez, A. Computational Geometry on Surfaces. (Springer Netherlands, 2001). https://doi.org/10.1007/9789401598095
 21.
Alt, S., Ganguly, P. & Salbreux, G. Vertex models: from cell mechanics to tissue morphogenesis. Philos. Trans. R. Soc. Lond. B. Biol. Sci. 372, 20150520 (2017).
 22.
Fletcher, A. G., Osterfield, M., Baker, R. E. & Shvartsman, S. Y. Vertex models of epithelial morphogenesis. Biophys. J. 106, 2291–2304 (2014).
 23.
Monier, B. et al. Apicobasal forces exerted by apoptotic cells drive epithelium folding. Nature 518, 245–248 (2015).
 24.
Osterfield, M., Du, X., Schupbach, T., Wieschaus, E. & Shvartsman, S. Y. Threedimensional epithelial morphogenesis in the developing Drosophila egg. Dev. Cell 24, 400–410 (2013).
 25.
Štorgel, N., Krajnc, M., Mrak, P., Štrus, J. & Ziherl, P. Quantitative morphology of epithelial folds. Biophys. J. 110, 269–277 (2016).
 26.
Hočevar Brezavšček, A., Rauzi, M., Leptin, M. & Ziherl, P. A model of epithelial invagination driven by collective mechanics of identical cells. Biophys. J. 103, 1069–1077 (2012).
 27.
Goodwin, K. & Nelson, C. M. Generating tissue topology through remodeling of cellcell adhesions. Exp. Cell Res. https://doi.org/10.1016/j.yexcr.2017.03.016 (2017).
 28.
Kim, H. Y., Varner, V. D. & Nelson, C. M. Apical constriction initiates new bud formation during monopodial branching of the embryonic chicken lung. Development 140, 3146–3155 (2013).
 29.
Khan, Z., Wang, Y.C., Wieschaus, E. F. & Kaschube, M. Quantitative 4D analyses of epithelial folding during Drosophila gastrulation. Development 141, 2895–2900 (2014).
 30.
Bielmeier, C. et al. Interface contractility between differently fated cells drives cell elimination and cyst formation. Curr. Biol. 26, 563–574 (2016).
 31.
Misra, M., Audoly, B., Kevrekidis, I. G. & Shvartsman, S. Y. Shape transformations of epithelial shells. Biophys. J. 110, 1670–1678 (2016).
 32.
Okuda, S., Inoue, Y., Eiraku, M., Adachi, T. & Sasai, Y. Vertex dynamics simulations of viscositydependent deformation during tissue morphogenesis. Biomech. Model Mechanobiol. 14, 413–425 (2015).
 33.
Honda, H. & Nagai, T. Cell models lead to understanding of multicellular morphogenesis consisting of successive selfconstruction of cells. J. Biochem. 157, 129–136 (2015).
 34.
Misra, M., Audoly, B. & Shvartsman, S. Y. Complex structures from patterned cell sheets. Philos. Trans. R. Soc. B. Biol. Sci. 372, 20150515 (2017).
 35.
Hannezo, E., Prost, J. & Joanny, J.F. Theory of epithelial sheet morphology in three dimensions. Proc. Natl Acad. Sci. USA 111, 27–32 (2014).
 36.
Murisic, N., Hakim, V., Kevrekidis, I. G., Shvartsman, S. Y. & Audoly, B. From discrete to continuum models of threedimensional deformations in epithelial sheets. Biophys. J. 109, 154–163 (2015).
 37.
Tan, R. Z., Lai, T. & Chiam, K.H. The role of apical contractility in determining cell morphology in multilayered epithelial sheets and tubes. Phys. Biol. 14, 46003 (2017).
 38.
Bokka, K. K., Jesudason, E. C., Warburton, D. & Lubkin, S. R. Quantifying cellular and subcellular stretches in embryonic lung epithelia under peristalsis: where to look for mechanosensing. Interface Focus 6, 20160031 (2016).
 39.
Condic, M. L., Fristrom, D. & Fristrom, J. W. Apical cell shape changes during Drosophila imaginal leg disc elongation: a novel morphogenetic mechanism. Development 111, 23–33 (1991).
 40.
Sun, Z. et al. Basolateral protrusion and apical contraction cooperatively drive Drosophila germband extension. Nat. Cell. Biol. 19, 375–383 (2017).
 41.
Rupprecht, J. F. et al. Geometric constraints alter cell arrangements within curved epithelial tissues. Mol. Biol. Cell 28, 3582–3594 (2017).
 42.
Honda, H., Nagai, T. & Tanemura, M. Two different mechanisms of planar cell intercalation leading to tissue elongation. Dev. Dyn. 237, 1826–1836 (2008).
 43.
Yin, X. et al. Engineering stem cell organoids. Cell. Stem. Cell. 18, 25–38 (2016).
 44.
Hendow, E. K. et al. Biomaterials for hollow organ tissue engineering. Fibrogenes Tissue Repair 9, 3 (2016).
 45.
Lecuit, T. & Wieschaus, E. Junctions as organizing centers in epithelial cells? A fly perspective. Traffic 3, 92–97 (2002).
 46.
Bertet, C., Sulak, L. & Lecuit, T. Myosindependent junction remodelling controls planar cell intercalation and axis elongation. Nature 429, 667–671 (2004).
 47.
Maybeck, V. & Roper, K. A targeted gainoffunction screen identifies genes affecting salivary gland morphogenesis/tubulogenesis in Drosophila. Genetics 181, 543–565 (2009).
 48.
Girdler, G. C. & Roper, K. Controlling cell shape changes during salivary gland tube formation in Drosophila. Semin. Cell. Dev. Biol. 31, 74–81 (2014).
 49.
Spencer, M. A., Jabeen, Z. & Lubensky, D. K. Vertex stability and topological transitions in vertex models of foams and epithelia. Eur. Phys. J. E. 40, 2 (2017).
 50.
Bi, D., Lopez, J. H., Schwarz, J. M. & Manning, M. L. A densityindependent rigidity transition in biological tissues. Nat. Phys. 11, 1074–1079 (2015).
 51.
Bosveld, F. et al. Modulation of junction tension by tumor suppressors and protooncogenes regulates cellcell contacts. Development 143, 623–634 (2016).
 52.
Lecuit, T., Lenne, P.F. & Munro, E. Force generation, transmission, and integration during cell and tissue morphogenesis. Annu. Rev. Cell. Dev. Biol. 27, 157–184 (2011).
 53.
Do Carmo, M. P. Differential geometry of curves and surfaces. Dover Phoenix Ed. 2, xiv, 474 (1976).
 54.
Cetera, M. & HorneBadovinac, S. Round and round gets you somewhere: collective cell migration and planar polarity in elongating Drosophila egg chambers. Curr. Opin. Genet. Dev. 32, 10–15 (2015).
 55.
HorneBadovinac, S. & Bilder, D. Mass transit: epithelial morphogenesis in the Drosophila egg chamber. Dev. Dyn. 232, 559–574 (2005).
 56.
Haigo, S. L. & Bilder, D. Global tissue revolutions in a morphogenetic movement controlling elongation. Science 331, 1071–1074 (2011).
 57.
Nelson, C. M. On buckling morphogenesis. J. Biomech. Eng. 138, 21005 (2016).
 58.
Wen, F. L., Wang, Y. C. & Shibata, T. Epithelial folding driven by apical or basallateral modulation: geometric features, mechanical inference, and boundary effects. Biophys. J. 112, 2683–2695 (2017).
 59.
Brodland, G. W. et al. CellFIT: a cellular forceinference toolkit using curvilinear cell boundaries. PLoS ONE 9, e99116 (2014).
 60.
Collinet, C., Rauzi, M., Lenne, P. & Lecuit, T. Local and tissuescale forces drive oriented junction growth during tissue extension. Nat. Cell. Biol. 17, 1247–1258 (2015).
 61.
Osterfield, M., Du, X., Schüpbach, T., Wieschaus, E. & Shvartsman, S. Y. Threedimensional epithelial morphogenesis in the developing Drosophila egg. Dev. Cell 24, 400–410 (2013).
 62.
Chen, H. H. & Brodland, G. W. Celllevel finite element studies of viscous cells in planar aggregates. J. Biomech. Eng. 122, 394–401 (2000).
 63.
Levayer, R., Hauert, B. & Moreno, E. Cell mixing induced by myc is required for competitive tissue invasion and destruction. Nature 524, 476–480 (2015).
 64.
Westerfield, M. TheZebrafish Book. A Guide for the Laboratory Use of Zebrafish (Danio rerio). 5th Edn, (Univ. Oregon Press, Eugene, 2007).
 65.
Royou, A., Sullivan, W. & Karess, R. Cortical recruitment of nonmuscle myosin II in early syncytial Drosophila embryos: its role in nuclear axial expansion and its regulation by Cdc2 activity. J. Cell. Biol. 158, 127–137 (2002).
 66.
Aldaz, S., Escudero, L. M. & Freeman, M. Dual role of myosin II during Drosophila imaginal disc metamorphosis. Nat. Commun. 4, 1761 (2013).
 67.
Morin, X., Daneman, R., Zavortink, M. & Chia, W. A protein trap strategy to detect GFPtagged proteins expressed from their endogenous loci in Drosophila. Proc. Natl Acad. Sci. USA 98, 15050–15055 (2001).
 68.
HernandezBejarano, M. et al. Opposing Shh and Fgf signals initiate nasotemporal patterning of the retina. Development 3933–3942 (2015). https://doi.org/10.1242/dev.125120
 69.
Lloyd, S. Least square quantization in PCM. IEEE Trans. Inform. Theory 28, 129–137 (1982).
 70.
Du, Q., Faber, V. & Gunzburger, M. Centroidal voronoi tessellations: applications and algorithms. SIAM Rev. 41, 637–676 (1999).
 71.
ArgandaCarreras, I. et al. Trainable Weka Segmentation: a machine learning tool for microscopy pixel classification. Bioinformatics. https://doi.org/10.1093/bioinformatics/btx180 (2017).
Acknowledgements
M.D.M.B. and A.V.E. work was funded by the Ministerio Español de Ciencia y Tecnología (grants numbers BFU201348988C21P and BFU201680797 to M.D.M.B.). A.V.E. is supported by a FPI studentship (BES2014068850) from the Ministerio de Economía y Competitividad. S.S. is supported by grant of the MICINN/FEDER to J.C.G. Hombría (BFU201676528P). F.C. work is funded by Spanish Government Grant BFU201455918 and BBVA Foundation Personal Grant IN[16]_BBM_BAS_0078. J.B. acknowledges core funding by Lehigh University. L.M.E. and P.G.G. are supported by the Ramón y Cajal program (PI13/01347); L.M.E., A.M., and C.G. work is funded by the Ministry of Economy, Industry, and Competitiveness grant BFU201674975P cofunded by FEDER funds. P.V.M. is supported by a contract cofunded by the Asociación Fundación Española contra el Cáncer and the Seville University. A.M.C. is supported by Fundación Pública Andaluza Progreso y Salud (Consejería de Salud), ref. PI00332014. C.F. and A.T. is supported by a contract from Sistema Nacional de Garantía Juvenil and Programa Operativo de Empleo Juvenil 2014–2020. We are thankful to Dr. Nicolas Gompel for the picture of the Cetoniidae: Protaetia (Potosia) speciose to illustrate the term “scutoid”. We thank Dr. Paco Martín, Dr. José LópezBarneo, Dr. Alberto Pascual, Dr. Colin Adrain and Dr. Matthew Freeman for helpful comments and corrections to the manuscript.
Author information
Affiliations
Contributions
L.M.E. designed the study with help from A.M., C.G., and J.B. P.G.G. and P.V.M. wrote the software for the mathematical model and performed analyses. A.T., C.F., A.M.C., M.L., A.V.E., M.B.G., O.S.P.H., F.C., S.S., and M.D.M.B., collaborated in the obtaining, measurement of parameters and analysis of the confocal images. All authors participated in the interpretation of results, discussions, and the development of the project. J.B. and L.M.E. wrote the manuscript with input from all authors.
Corresponding authors
Ethics declarations
Competing interests
The authors declare no competing interests.
Additional information
Publisher's note: Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Rights and permissions
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/.
About this article
Cite this article
GómezGálvez, P., VicenteMunuera, P., Tagua, A. et al. Scutoids are a geometrical solution to threedimensional packing of epithelia. Nat Commun 9, 2960 (2018). https://doi.org/10.1038/s41467018053761
Received:
Accepted:
Published:
Further reading

The Voronoi theory of the normal liver lobular architecture and its applicability in hepatic zonation
Scientific Reports (2021)

Establishment of a morphological atlas of the Caenorhabditis elegans embryo using deeplearningbased 4D segmentation
Nature Communications (2020)

Selfsustained planar intercalations due to mechanosignaling feedbacks lead to robust axis extension during morphogenesis
Scientific Reports (2020)

Universal hidden order in amorphous cellular geometries
Nature Communications (2019)

A 3D cell shape that enables tube formation
Nature (2018)
Comments
By submitting a comment you agree to abide by our Terms and Community Guidelines. If you find something abusive or that does not comply with our terms or guidelines please flag it as inappropriate.