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 great-circles.

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 procedures64. For each type of sample different immunohistochemistry conditions were used.

For Drosophila larval salivary glands we used the sqh-GFP line65 and Cy3-labeled 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 Fluoromount-G (Southern Biotech) using double sided tape to avoid the squashing of sample66. Images were taken with a Nikon Eclipse Ti-E 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:n-heptane 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 anti-Disc 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 ACS-APO 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 Resille-GFP (II) was used67. For ovaries staining 1–2-day-old females were fattened on yeast for 48 h before dissection. Ovaries were dissected in PBT (phosphate-buffered 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 anti-GFP-FITC (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 TCS-SPE confocal equipped with a ×40 (1.15 NA) ACS-APO oil objective.

For Zebrafish embryos membranes and nuclei of zebrafish embryos were labelled in vivo by microinjecting a bicistronic mRNA (GAP43-GFP:h2bRPF) at one-cell stage. mRNA for microinjection was synthesised using the mMessage Machine kit (Ambion), following manufacturer’s instructions. Embryos were grown up to pre-gastrula stage, fixed in 4% paraformaldehyde and embedded in 1.5% low-melting 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 four-cells 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 straight-line tools respectively. (Fig. 3c, d).

Percentage of scutoids: We quantified the percentage of central cells participating in at least one transition along the apico-basal axis using all the sections of the confocal stack.

In addition, three geometric properties in the salivary glands images were measured to feed the line-tension 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 four-cells motif (Fig. 4a, c). We used FIJI’s straight-line tool.

Width (w) and height (h): We measured the “width” and “height” from each four-cells motif in basal and apical surfaces. The “width”/“height” is defined as the distance between two vertices in the quadrilateral that defines a four-cells 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 straight-line 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 four-cells motifs taking as a reference the basal surface. We chose ten four-cells 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 four-cells motif in basal (using FIJI’s polygon selection tool) and dividing it by the area of the four-cells 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 four-cells 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 line-tension 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 X-axis, in anterior and posterior, and the Y-axis. 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 apico-basal transitions and 10 with no transition. To further analysed them, we segmented the four-cell motifs in their apical and basal frames, capturing the following measurements: edge angle and length regarding the X-axis 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 extreme20 (Supplementary Fig. 1). We have used Matlab R2014b (Mathworks) as our computational tool. Taking as a reference the previous 2D Voronoi model established in7, 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 diagram69,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 [513-1024] x [1-4096] 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 cells7. 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 apico-basal transitions (Fig. 3b, Supplementary Fig. 2a and Supplementary Data 1) and obtained the percentage of scutoids. We discarded the edges shared by non-valid 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 line-tension 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 four-cells motifs and we measured the distances between them. The four-cells motifs used in our analysis kept their four-cells together in apical and basal surfaces and we discarded the motifs that did not satisfy this requirement. We also disregarded the motifs with non-valid cells on any of these surfaces and we classified our measurements: four-cells motifs with and without apico-basal 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 shape53. However, in the cases of four-cell motifs with very short edges, we found false positives (apico-basal 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 four-cells 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 non-null 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 apico-basal 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 apico-basal 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 balloon-like spheroid (X radius inner layer = 1.5) and a more prolate spheroid, ‘Zeppelin’-shaped, (X radius inner layer = 2), with three different apico-basal 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:

$$\frac{x}{a} + \frac{y}{b} + \frac{z}{c} = 1,$$ (1)

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 apico-basal 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 Z-axis (first from Z = zRadius to Z = 0, and second from Z = zRadius to Z = 0) and two others along the Y-axis (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 four-cells 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 four-cells 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 non-valid motifs (containing non-valid cells). We also distinguished between transition and no transition motifs. As explained above, we only selected four-cell 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 surfaces53. 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 two-sample 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 Segmentation71 (FIJI) a machine learning-based 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 post-processing 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 semi-automatic 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 four-cell motif to visualize the surfaces of the cells along the apico-basal 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.

Line-tension minimization model

The existence of tensile forces controlling cell interactions, either driven by cadherin and/or acto-myosin activity, can be described through the following line-tension energy functional:

$$\begin{array}{*{20}{c}} E & = & {\sigma \left[ {l_w + 2\sqrt {\left( {w - l_w} \right)^2 + h^2} } \right]H\left( {l_w} \right)H\left( {w - l_w} \right) } \\ {} & {} & +{\sigma \left[ {l_h + 2\sqrt {(h - l_h)^2 + w^2} } \right]H(l_h)H(h - l_h)} \end{array},$$ (2)

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 four-cell configuration reads,

$$\hat E = \frac{{\left[ { - l + 2\sqrt {(1 + l)^2 + {\it{\epsilon }}^2} } \right]H( - l)H(1 + l) + \left[ {l + 2\sqrt {(h - l)^2 + 1} } \right]H(l)H(h - l)}}{{2\sqrt {{\it{\epsilon }}^2 + 1} }},$$ (3)

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,

$$\frac{{\partial l}}{{\partial t}} = - \frac{{\partial \hat E}}{{\partial l}},$$ (4)

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 vice-versa). 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,

$$\hat E_{{\rm exp}} = \frac{{\hat L_T}}{{2\sqrt {1 + \left\langle {{\it{\epsilon }}} \right\rangle ^2} }}$$ (5)

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,

$$\begin{array}{*{20}{c}} {\hat E_w}={\hat E_{{\rm exp}}{\rm sin}\theta } \\ {\hat E_h}={\hat E_{{\rm exp}}{\rm cos}\theta } \end{array}$$ (6)

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.