Skip to main content
Download PDF

Extracting multiple surfaces from 3D microscopy images in complex biological tissues with the Zellige software tool

BMC Biology volume 20, Article number: 183 (2022) Cite this article

Abstract

Background

Efficient tools allowing the extraction of 2D surfaces from 3D-microscopy data are essential for studies aiming to decipher the complex cellular choreography through which epithelium morphogenesis takes place during development. Most existing methods allow for the extraction of a single and smooth manifold of sufficiently high signal intensity and contrast, and usually fail when the surface of interest has a rough topography or when its localization is hampered by other surrounding structures of higher contrast. Multiple surface segmentation entails laborious manual annotations of the various surfaces separately.

Results

As automating this task is critical in studies involving tissue-tissue or tissue-matrix interaction, we developed the Zellige software, which allows the extraction of a non-prescribed number of surfaces of varying inclination, contrast, and texture from a 3D image. The tool requires the adjustment of a small set of control parameters, for which we provide an intuitive interface implemented as a Fiji plugin.

Conclusions

As a proof of principle of the versatility of Zellige, we demonstrate its performance and robustness on synthetic images and on four different types of biological samples, covering a wide range of biological contexts.

Background

The interplay between gene regulatory networks and physical forces in driving collective cell behaviors is key to tissue morphogenesis during development and to tissue homeostasis throughout life. Recent quantitative studies of epithelial morphogenesis have begun to unravel the basic cellular and physical principles of tissue development, by providing the tools to integrate multiple scales of tissue dynamics [1,2,3,4]. These tools are instrumental to quantify how cell shape changes, cell divisions, cell rearrangements, and cell extrusions contribute to tissue remodeling, and to establish data-driven computational models of tissue morphogenesis.

Quantitative analysis of an epithelium starts with the extraction of its apical surface from 3D-microscopy images (z-stacks of xy-optical sections) encompassing the volume immediately surrounding the epithelium. To this end, the abovementioned studies relied on specific preparations of the specimen, allowing to expose the entire epithelial surface labeled with junctional fluorescent markers to reveal the network formed by epithelial cell-cell contacts. Once the epithelial surface has been extracted, automated cell segmentation and cell contour tracking tools can be used to follow the dynamics of every cell within the epithelium. However, extracting the epithelial surface remains a difficult task because this surface is usually not flat (it is best modelled as a curved surface, or 2D submanifold embedded in 3D space), and is often surrounded by other biological structures such as cell layers, acellular membranes, extracellular matrix, and vesicles that hamper its visualization and reconstruction. Even more difficult situations occur when the epithelium surface cannot be readily exposed such as in epithelial organoids in which the epithelial surface faces the lumen inside the cyst. More generally, there is a need for efficient tools allowing one to discriminate a surface of interest from surrounding structures of high intensity and contrast. This is required to study how epithelial cell dynamics account for tissue shape changes in complex specimens, in a wide range of contexts including developmental biology, host-pathogen interaction, and tissue regeneration studies.

Another experimental limitation of previous studies is that some of the structures surrounding the epithelium exert external physical constraints that are known to critically affect epithelial morphogenesis by directing cellular dynamics and signalling pathways [1, 5,6,7]. To understand the physical forces controlling tissue morphogenesis, it is thus essential to also characterize how the dynamics of these extra-epithelial surfaces relate to that of the epithelium (see [8] for review). This calls for the development of dedicated tools allowing the automated extraction of information from several surfaces of interest in a given sample, since the sheer volume of the data challenges any attempt at a manual analysis [9].

Several surface extraction tools have been developed, some of which are available as open access software, such as PreMosa [10], FastSME [11], and LocalZProjector [12]. These tools focus on the extraction of a single, near-horizontal epithelial layer, which is assumed to (i) be sufficiently smooth, (ii) show enough contrast against surrounding background signals, and (iii) cover the entire image field-of-view. Specifically, it is assumed that the fluorescent marker used to label the epithelial cell network should provide the highest contrast in the image and allow to select it out from autofluorescent extracellular structures such as the cuticle in flies, or other acellular membranes in mammalian epithelia. The surface is then localized using heuristic algorithms based on the detection of the voxels of maximum contrast and/or brightness. However, applying these tools on more complex biological images with several epithelia of weaker contrast often leads to incorrect localization of the surface of interest, and its blending with the nearby unwanted biological structures.

MinCostZ on the other hand is the only available open-source tool that allows the extraction of up to two surfaces from a 3D stack, and imposes explicit continuity constraints on the reconstructed surfaces. MinCostZ surface extraction relies on a previously developed formulation of the problem as a graph-cut optimization [13]. It is implemented as an ImageJ plugin [14], taking as control parameters, the number of surfaces to extract, the maximum slope, and the range of distances allowed between the surfaces, as well as some user-defined cost function that should reflect the characteristics of the surfaces in term of signal intensity, contrast, and texture. Despite its interest, this approach remains computationally costly and difficult to apply in practice due to the non-trivial choice of the cost function and the need to know beforehand the relative positions of the surfaces to be extracted.

Alternatively, one can segment the surfaces of interest by using supervised machine learning tools such as the software solutions Weka [15] or Ilastik [16], as proposed in the ImSAnE surface reconstruction framework [17]. A deep learning approach, using a network of the U-net type to segment the voxels belonging to a single surface of interest, has also recently been reported [18]. While promising as they can provide state of the art segmentations of epithelial surfaces in difficult imaging conditions, machine learning approaches require the prior manual annotation of a sufficiently large set of surfaces to generate suitable training sets. This process can be very time consuming, often necessitating several rounds of trials and errors to obtain satisfactory results, without guarantees to be generic, i.e., to generalize to a wide range of datasets. So far, no solution to the multiple surface extraction problem has been proposed, which is satisfactory both in terms of genericity and ease of use.

However, such a tool is highly desirable for modern biology studies. Indeed, tissue organization in the context of developmental biology emerges from the interaction of several neighboring structures [8] through the interplay of molecular signals [19, 20], as well as electrical [21, 22], hydraulic [23, 24], and mechanical contact interactions [25, 26]. The importance of such interaction is exemplified in embryonic explanted tissue cultures that develop abnormally when separated from their neighboring structures [27]. Similarly, in the context of tissue engineering, stem cell-derived aggregates harbor various types of tissues surrounding the genuine organoid, and these tissues presumably influence organoid shape, fate, and differentiation (see [28] for review). The ability to simultaneously study the dynamics of neighboring structures together with the structure of interest is therefore essential for an integrated understanding of tissue development, and for any attempt to harness tissue self-organization in vitro.

Here, we introduce Zellige, a tool based on a novel constructive approach that allows the automatic extraction of a non-prescribed number of surfaces from a 3D image. Unlike most other tools, Zellige does not rely on the localization of voxels of highest contrast or brightness. Rather, it performs the surface reconstruction in two steps: it first classifies voxels most likely belonging to surface structures by a combination of maxima detection along the z-axis and the application of two brightness and contrast classifiers; it then extracts each surface present in the image volume with an assembly algorithm implementing explicit continuity constraints. This two-step approach allows the reconstruction of multiple surfaces of different contrast characteristics to be selectively reconstructed.

To do this, the user is only required to adjust a small set of intuitive control parameters, a task largely facilitated by a user-friendly interface implemented as a plugin for the open-source Fiji platform [29]. We tested the performance and robustness of Zellige for multiple surface extraction by applying it to synthetic and 3D microscopy images containing multiple surfaces of interest of widely varying texture and contrast. Among the four biological samples we considered, three previously required laborious manual extraction of the epithelial surface (fly pupa histoblast, mouse cochlea, and inner ear organoid) therefore greatly challenging the study of the cellular mechanisms involved in the development of these large organs [9]. We also analyzed a bronchial epithelium specimen that highlights the use of Zellige in a recent study to address how SARS-CoV-2 infection interferes with the epithelial function [30]. These experiments demonstrate the ability of the approach to extract several surfaces of potentially very low contrast, selectively from other highly contrasted and complex structures, with a single set of reconstruction parameters. A sensitivity analysis also reveals a high robustness of Zellige against small variations of these parameters. This will make it a tool of choice in terms of versatility and ease of use for the investigation of biological surfaces.

Implementation

Zellige was devised with the goal of achieving accurate segmentation of multiple biological surfaces from 3D confocal images. Unlike other existing surface extraction tools, it makes no assumption on the number of surfaces to be extracted and does not require the surfaces of interest to be the structures of highest contrast in the image. Zellige is written in Java, relying on the ImgLib2 library [31] and is distributed as a Fiji plugin, with a graphical user interface (GUI) designed to allow users to quickly find a good set of extraction parameters for a given image.

Zellige extracts each surface present in the image in the form of a height map (or z-map), that is, a mapping z=h(x,y), which associates to each point (x,y) over which the surface projects, the z-coordinate of the unique voxel (x,y,z) belonging to the surface. Each extracted height map is then used to produce a 2D projection of the 3D stack restricted to a small subvolume (of user-selected width) centered around the corresponding surface. To achieve a robust extraction, Zellige proceeds in two algorithmic steps, which are only outlined below (see Additional file 1: Supplementary Note 1 for implementation details) [32].

In the first step, or surface voxel selection step, a segmentation is applied to the 3D image to select voxels that likely belong to a surface of interest (Fig. 1 and Additional file 2: Figure S1A). These putative surface voxels are detected as local maxima of image intensity along the z-axis, after applying a 3D Gaussian filter to the image to reduce noise artifacts, followed by two independent binary classifiers, one based on voxel contrast and the other one on voxel intensity (see Additional file 1: Supplementary Note 1 for details). Five adjustable parameters control the selection step: two threshold parameters (TA and Totsu) control the strength of the binary classifiers applied on contrast and intensity, respectively, and three parameters (Smin, σxy, and σz) control clean-up operations applied at the end of the classification (removal of small isolated spots, and local averaging along the xy plane and the z-axis, respectively).

Fig. 1
figure 1

Flowchart of Zellige’s algorithmic steps. Surface voxel selection (step 1), surface assembly in the form of a height map (step 2), and subsequent projection localized to the height map, are schematically depicted in the case of a 3D image containing 4 surfaces of interest

Full size image

In the second step, or surface assembly step, an iterative algorithm is used to extract the height maps of each of the surfaces present in the image (Fig. 1 and Additional file 2: Figure S1A). The assembly starts by grouping together contiguous chains of putative surface voxels within each orthogonal xz section (alternatively, within each yz section) of the 3D image, in order to form a set of building blocks referred to as orthogonal surface elements (OSEs). These building blocks are then used to assemble the surfaces, in a process analogous to jigsaw puzzles, where OSEs adjacent to the surface boundary are added if they match this boundary, and rejected otherwise, according to matching criteria that impose explicit continuity constraints to the surface under construction. (See Additional file 1: Supplementary Note 1 for details on these matching criteria and on the construction of the OSEs.) The surface assembly proceeds until no matching OSE can be found (Additional file 2: Figure S1). In order to increase the robustness of the surface assembly step, Zellige applies it in two rounds, constructing the OSEs within yz sections only and performing the assembly along the x-axis during the first round, and constructing OSEs within the xz sections and proceeding along the y-axis during the second round (see Additional file 1: Supplementary Note 1). Each round is controlled by 3 adjustable parameters: a threshold parameter (0 ≤ TOSE ≤1) sets a minimum size for the building blocks that can be used as seeds to initiate the assembly of a surface; and two other parameters (0 ≤ R ≤ 50 and 0 ≤ C ≤ 1) set the matching constraints used to accept or reject the addition of OSEs to a surface. The assembly step is thus controlled overall by 6 parameters, i.e., two groups of 3 parameters (TOSE1, R1, C1) and (TOSE2, R2, C2) controlling the first and second assembly rounds, respectively.

Finally, the height maps of each of the reconstructed surfaces are used to obtain a corresponding 2D projection (Fig. 1). In practice, a maximum projection restricted to a subvolume of width δz (δz being a user-defined parameter) centered around the surface of interest is performed.

Results and discussion

Proof of concept of multiple surface extraction on a synthetic image

Figure 2 shows the results produced by Zellige on a phantom 3D image [33] containing three distinct synthetic surfaces generated as described in Additional file 1: Supplementary Note 2 [34]. The three surfaces are extracted with little errors. We assessed the quality of the reconstruction by comparing each of the height maps produced by Zellige to the corresponding ground truth (GT) height map, which is exactly known in this case (Fig. 2A–C, and Additional file 1: Supplementary Note 3) [12, 35]. For the three surfaces, the absolute value |h(x,y) – hGT(x,y)| of the difference between the reconstructed and the GT height maps is ≤ 1 (in pixel units) over >99% of the GT pixels (Fig. 2D–E), with a root mean square error value (RMSE; defined over the pixels common to the reconstructed and the GT height maps, see Eq. 1 in Additional file 1: Supplementary Note 3) of ≤ 0.6 in pixel units, showing that the surface localization is highly accurate. In addition, the coverage, which measures the proportion of the reconstructed surface relative to the GT, is near 100% for the three surfaces. To achieve these results, the control parameters of the two steps of the surface extraction were adjusted manually to some adequate reference values (see Additional file 1: Supplementary Table S1) using the Zellige Fiji interface. Only the parameters controlling the voxel classification step (amplitude and Otsu threshold parameters TA and Totsu, minimal island size Smin, and smoothing parameters σxy and σz) did actually require a modest adjustment. The parameters of the surface assembly step (parameters TOSE1, R1, C1 and TOSE2, R2, C2 of the 1st and 2nd construction rounds, respectively) were set to their default reference values (see Additional file 3: Figure S2 and Additional file 1: Supplementary Note 4) and did not need to be adjusted.

Fig. 2
figure 2

Multiple surface extraction on a synthetic 3D image. A,B The image contains 3 phantom surfaces (S1, S2, S3) of different shapes (sinusoidal, flat, and paraboloidal, respectively), and different textures (surface S1 has constant intensity, while surfaces S2 and S3 are supported by Voronoi meshes of different cell sizes). C 3D representations of the height maps extracted by Zellige (in green) and of the ground truth (GT, in blue) height maps of surfaces S1, S2, and S3. D Error maps displaying the distance along the z-axis between the reconstructed and GT height maps for surfaces S1, S2, and S3. E Projections of the 3D image localized to the different surfaces S1–S3 (maximum intensity projections over a subvolume of a width δz=1 pixel above or below the corresponding height maps). Upper and lower panels show the projections based on the GT and the reconstructed height maps, respectively

Full size image

Thus, using a single set of control parameters, Zellige can extract multiple surfaces of various shapes and textures with very little error, without requiring the user to provide information about their number or relative position, nor about their shape or texture characteristics.

Performance of Zellige on biological samples

Example 1: extracting multiple surfaces from an image of a pupal fly specimen

Over the past few decades, the Drosophila model has been invaluable to decipher the molecular and cellular mechanisms underlying organ embryogenesis [36, 37]. Epithelium morphogenesis studies not only revealed the importance of mechanical stresses (including stress boundary conditions) and planar polarity signalling on cell dynamics to generate tissues of reproducible sizes and shapes, it also highlighted the importance of extracellular matrix attachments in constraining the tissue stresses that guide patterning [1, 38]. At the pupal stage, the fly undergoes dramatic remodeling of its larval organs into adult organs. Large-scale tissue flows initiate at a timing that coincides with molting, when the epithelium contracts away from the overlaying cuticular sac, a protective acellular membrane that imposes mechanical constraints to the tissue.

Figure 3 shows the results of applying Zellige on a 3D image of a Drosophila pupa acquired with a spinning disk confocal microscope [39]. The sample expresses Ecadherin-GFP, a fluorescent marker of cell-cell junctions, and encompasses a portion of the pupa’s abdomen and a small portion of its wing. Four surfaces of interest can be identified, with varying signal intensities, noise levels, and features (Fig. 3A,B). The abdomen is formed of an epithelium (surface S2) overlaid by a cuticle (surface S1). Lying just beneath these two surfaces, one can observe globular structures showing in some places a higher intensity than the signal coming from the surfaces. The wing also consists of an epithelium of low-intensity signal (surface S4), and an overlying cuticle (surface S3). These two surfaces are relatively flat, except for surface S3 which is very steep near one of its edges.

Fig. 3
figure 3

Fly specimen. A,B Volume rendering (A) and orthogonal sections (B) of a 3D image of fly embryo taken around 24 h after puparium formation, covering a portion of the abdomen (showing histoblast cells and larval cells), and a portion of the developing wing. Scale bar 50 μm. Four surfaces of interest may be identified in the dataset (of dimensions 1200 × 1200 × 51 voxels): surfaces S1 and S2 are relatively close to one another and located within overlapping z-ranges (8 ≤ z ≤ 50 and 20 ≤ z ≤ 50, respectively). Surfaces S3 and S4 (located in the z-ranges 42 ≤ z ≤ 50 and 9 ≤ z ≤ 50, respectively) are relatively far from each other and can nearly be separated by a plane. C 3D representations of the height maps extracted by Zellige (in green) and of the ground truth height maps (GT, in blue) of surfaces S1–S4. The reconstructed height maps of all surfaces S1–S4 cover >93% of the area of the corresponding GT (cf. Additional file 3: Figure S2 and Additional file 1: Supplementary Table S1). To reduce the staircase artifacts (more or less visible depending on the surface) due to the digitization of the GT and reconstructed height maps, all height maps were smoothed with a 2D Gaussian filter with a standard radius of 5 pixels (cf. Additional file 1: Supplementary Note 1). D Error maps (color-coded distance along the z-axis between the reconstructed and the GT height maps) plotted for each of the reconstructed surfaces. The large majority of pixels on the reconstructed height maps (98%, 96%, 91%, and 99% for surfaces S1 to S4, respectively) display errors of <2 pixels. The height maps of surfaces S1, S2, S4 show subpixel accuracy on average (RMSE < 1), while that of surface S3 is slightly less accurate (RMSE = 1.25). E Projections of the 3D image localized to the different surfaces S1–S4 (in this and all subsequent figures, these are maximum intensity projections over a subvolume of width δz=±1 pixel above or below the corresponding height maps). Upper and lower panels show the projections based on the GT and the reconstructed height maps, respectively. Scale bar 50 μm

Full size image

Figure 3C shows a 3D graphical representation of the height maps reconstructed by Zellige (green) and those reconstructed by an expert biologist (blue), taken as ground truth (GT). While these height maps clearly show greater roughness than those of the synthetic surfaces presented earlier, they could again be obtained with a single set of control parameters that were adjusted manually with the Zellige interface (see Additional file 1: Supplementary Table S1). We observe an excellent match between the four reconstructed and corresponding GT height maps, despite the rather complex topography of surfaces S1 and S2 (with slopes reaching up to 45°), and the near-vertical inclination of surface S3 at its boundary. Yet, small deviations may be seen in the regions of highest slope of the surfaces. Some of these deviations are likely attributable to uncertainties in the definition of the GT height maps, whose accuracy depends on the expert.

Figure 3D shows the differences between the reconstructed and GT height maps, plotted as color-coded error maps. These differences are <2 (in pixel units) for most pixels, while some regions of higher error values can be seen locally in surfaces S1 and S2, and at the boundary of surface S3. Note that for surfaces S2 and S4, which contain junctional epithelial meshes composed of larger and smaller cells, respectively, the GT height map encompass not only the mesh but also the interior of the cells, where no junctional signal is detected. The distance calculated inside the cells is thus more subjected to intensity fluctuations, especially for surface S2. Nonetheless, the RMSEs of surfaces S1, S2, and S4 are less than 1, showing that on average the reconstructed height maps match the corresponding ground truth with subpixel accuracy (i.e., RMSE < 1). The higher RMSE (1.25) of surface S3 is largely due to the region of steep slope at the edge of this surface (yellow region on the error map for this surface, Fig. 3D). The coverage of the reconstructed height maps is excellent (≥ 96%) for surfaces S1 and S2, and slightly lower, but still very good (≥ 93%) for the smaller surfaces S3 and S4. Figure 3E shows the 2D projections of the 3D image obtained for each of the reconstructed surfaces and for the corresponding ground truths. The inaccuracies visible on the error maps (see Fig. 3D) do not significantly impact these projections, which appear very similar to the projections obtained with the corresponding ground truths. Thus, while the biological sample contains significant noise and shows a much more variable contrast (especially with the presence of high-intensity globular structures underneath surfaces S1 and S2), Zellige makes it possible to segment these surfaces selectively, with a quality of segmentation comparable to that obtained by manual expert segmentation.

This possibility brings several perspectives that are not offered by single surface extraction algorithms. First, it opens the possibility to systematically study the tissue axial movements (along z) relative to the cuticle during molting, allowing for example to gain insights into the early tissue contraction of the wing hinge that acts as a mechanical inducer over the wing blade [8]. Second, Zellige makes it possible to automatically extract structures such as the abdomen epithelium, which is usually segmented manually [9], due to the difficulty to separate the large larval cells from the cuticle mesh and from other globular structures (such as fat bodies or macrophages) present underneath the epithelium. All these structures become intertwined when using other extracting tools. In this context, Zellige opens new opportunities to study collective cell behavior during epithelial morphogenesis in vivo, and to integrate in the analysis the surrounding surface-like structures involved in the mechanics of the system.

Example 2: extracting a thin cochlear epithelium surface from a multilayer dataset

As the first model in which planar cell polarity signalling was shown to be conserved in vertebrates [40], the mammalian auditory organ, the cochlea, is arguably our most valuable model to study epithelial patterning and morphogenesis beyond the fly and zebrafish [41, 42]. Cochlear morphogenesis involves complex and tightly controlled patterning processes during which the cochlear sensory epithelium extends and develops its characteristic coiled snail shape, while adopting a striking cellular mosaic organization, with graded changes of morphogenetic parameters along the cochlea [43, 44]. These morphogenetic processes are well recapitulated in organotypic cultures, on the condition that the mesenchyme that underlies the epithelium be preserved. The cultures are then amenable to live imaging [42, 45], pharmacological [46], and genetic manipulations.

Figure 4A shows a confocal swept-field microscope acquisition of an embryonic mouse cochlea [47]. The sample contains only one surface of interest, the cochlear epithelium, but this surface lays on top of a thick tangled mesh of non-epithelial cells originating from the mesenchyme. The whole biological tissue is stained for filamentous actin (F-actin) using phalloidin. The epithelium surface presents a non-uniform signal included in a small z-range (6 ≤ z ≤ 10), and a mesh of very heterogeneous size. Between sections z = 10 and z = 14, one can observe the basolateral region of the epithelial cells, also stained for F-actin. The particularity of this sample is that the mesenchyme presents an intense and contrasted signal over a wide range of z-values (14 ≤ z ≤ 43). This makes it challenging to extract the surface of the epithelium, which is characterized by low intensity and low contrast.

Fig. 4
figure 4

Cochlea specimen. A Volume rendering of a 3D confocal swept-field image of the mouse cochlear embryo on embryonic day E14.5. The dataset (of dimensions 1024 × 1024 × 45 voxels) shows a portion of the sensory epithelium (at the topmost sections of the stack) and the underlying non-cellular layer of mesenchyme on which the organ develops. Both structures are stained with phalloidin to reveal F-actin. The surface of interest is the epithelium surface, harboring the sensory and supporting cells under differentiation. The mesenchyme layer is not (strictly speaking) assimilable to a surface, but it produces a strong background signal nearby the surface of interest, hampering its extraction. B 3D representations of the height map extracted by Zellige (in green) and the GT height map (in blue), of the epithelium surface. C Color-coded error map of the reconstructed height map, which shows subpixel accuracy (errors <1) over a large majority (83%) of pixels, while being close to 1 pixel on average (RMSE ~ 1.1). D Projections localized to the GT height map of the epithelium surface (left most panel), and to the height maps extracted with the four different algorithms: FastSME, LocalZProjector, PreMosa, and Zellige. Only Zellige correctly extracts the surface of the epithelium in this example. Scale bar 40 μm

Full size image

Figure 4B shows a 3D representation of the height map reconstructed by Zellige and the corresponding GT height map (again reconstructed manually by an expert). On the corresponding error map (Fig. 4C), most (~83%) pixels of the reconstructed height map show subpixel accuracy (with distances < 1 to the corresponding pixels of the GT height map). The errors are greater in regions where the cell size is larger, as well as in the area where the signal intensity is very low. However, they remain smaller than 2 for > 95% of the pixels. This result is consistent with the low value of the RMSE (1.1). The surface is also reconstructed with an excellent coverage (> 99%, Additional file 1: Supplementary Table S1).

As the sample contains a single epithelial surface of interest, we compared the performance of Zellige with three other software that can extract only a single surface (Fig. 4D). The projections of PreMosa, FastSME, and LocalZProjector completely miss the epithelium, despite our attempts to find an optimal set of parameters using a multi-dimensional parameter screening (as done previously in Herbert et al. [12]). Only regions of high contrast corresponding to the mesenchyme are projected. In contrast, Zellige generates a projection very close to the ground truth. This demonstrates the efficiency of Zellige to selectively extract a low contrast surface, despite the presence of several structures of higher contrast. Indeed, Zellige detects every structure as a possible surface seed without any assumption on its contrast, and only extends this seed into a surface if enough spatial continuity is found in the surrounding signal. This feature allows to separate individual surfaces from other structures spatially, which should greatly facilitate the analysis of live imaging experiments.

Example 3: extracting a single bronchial epithelial surface rendered abnormally rough by SARS-CoV-2

Recently, we used Zellige to extract the surface of a primary culture of bronchial epithelial cells following infection by the SARS-CoV-2 virus [30]. The infection causes the surface of the epithelium to become abnormally rough due to cell damages as seen from discontinuities within the cell layer. The sample we chose from this study is a 3D confocal image of the epithelium responding to SARS-CoV-2 infection [48] (Fig. 5A). The surface of interest in this image corresponds to the layer of epithelial cells stained for the tight junction protein Zona Occludens-1 (ZO-1). The surface roughness causes the network of junctions to extend over the height of the z-stack, with a signal of varying intensity (Fig. 5A). In addition, the junctional network remains non-planar even at the level of a single cell, hence violating the smoothness condition commonly assumed to hold in the context of epithelial surface extraction. We also observe the presence of nearby punctiform structures of high contrast that are mainly located outside of the epithelium surface. This sample therefore provides an example of a surface with a complex landscape, interspersed with a constellation of signals which may interfere with the surface segmentation.

Fig. 5
figure 5

Human bronchial epithelial cells infected by SARS-CoV-2. A Volume rendering and individual sections of a confocal 3D image of a primary culture of bronchial epithelial cells 4 days after it was infected by the SARS-CoV-2 virus. The dataset (of dimensions 1024 × 1024 × 15 voxels) covers a portion of the epithelium immunostained for the tight junction protein ZO-1. Notice the roughness of the epithelium surface and the presence of anomalous bulges resulting from the SARS-CoV-2 infection. Scale bar 10 μm. B 3D representations of the height map extracted by Zellige (in green) and the GT height map (in blue), of the epithelium surface. C Color-coded error map of the reconstructed height map. Despite its roughness, the surface of interest is reconstructed with subpixel accuracy over the majority (71%) of pixels, as well as on average (RMSE ~ 0.81). D Projections localized to the GT height map of the epithelium surface (leftmost panel), and the height maps extracted with the four different algorithms: FastSME, LocalZProjector, PreMosa, and Zellige. Scale bar 30 μm

Full size image

The 3D representation of the reconstructed and corresponding GT height maps (Fig. 5B) makes it possible to appreciate the roughness of the surface of interest. The two height maps overlap quite satisfactorily. As shown on the error map (Fig. 5C), a large majority (71.1%) of pixels of the reconstructed height map show errors smaller than 1 pixel (~ 96% of them showing errors smaller than 2 pixels). The error is however larger in regions where the cell size is larger, as well as in areas where the ZO-1 signal intensity is very low, preventing a complete reconstruction of the junctions. Nevertheless, the overall RMSE remains small, with a value of 0.81 (Fig. 5). The coverage of the reconstructed surface is also excellent (98% of the GT height map), despite the abovementioned discontinuities. Figure 5D shows the comparison of the projection generated by Zellige by manually adjusting parameters to those produced by PreMosa, FastSME, and LocalZProjector after a search of an optimal set of parameters (see Herbert et al. [12]). The projection generated by PreMosa misses many junctions of the epithelial network, but it removes quite well the punctiform signal originating from other optical sections. FastSME performs better than PreMosa in reconstructing the junctions, but it produces a projection where the punctiform signal remains strong. In contrast, Zellige and LocalZProjector manage to both reconstruct the surface well and to filter out the punctiform signal quite effectively. This result demonstrates the efficiency of Zellige to extract a surface with complex topography by excluding intense and contrasted spurious signals away from the epithelium surface.

Example 4: extracting the apical and basal layers of a dome-shaped epithelium (developing inner ear organoid)

Organoids are stem cell-derived and self-organizing 3D tissue structures that can mimic certain organ structures. They have emerged as promising in vitro models for developmental biology research, as well as biomedical translational research applications. Here we take the example of mouse stem cell-derived inner ear organoids that form vesicular structures composed of an epithelium harboring sensory cells. These organoids are part of a cellular aggregate that also contains other tissues such as the mesenchyme [49] adjacent to the organoids. The epithelial cells of the forming inner ear organoids acquire a basal-apical polarity, with their apical side facing the lumen of the organoid, and their basal side facing outwards. The apical junctional network of the epithelium is difficult to visualize in microscopy images as it is seen from below, through the basal layer. Another difficulty is the spherical geometry of the vesicle system, which makes the epithelial surface of interest difficult to extract in regions of high inclination relative to the focal plane.

Figure 6 shows the result of applying Zellige on a 3D confocal microscopy image of half of a developing inner ear organoid at 14 days of culture, a stage at which markers characteristic of the mouse otic vesicle can be detected [50]. The sample was fixed and stained for F-actin to visualize all cellular structures including the epithelium. Two surfaces of interest can be identified, namely the basal side of the epithelium and the apical junctional network (Fig. 6A,B). Both surfaces are mesh-like structures of high inclination, high signal intensity, and high contrast. The vesicle lumen also contains cell debris of high intensity and contrast that are not part of any surface of interest.

Fig. 6
figure 6

Organoid specimen. A,B Volume rendering (A) and orthogonal sections (B) of a confocal 3D image of a (half of) inner ear organoid, which has been fixed and stained with phalloidin to reveal F-actin. The dataset (of dimensions 520 × 465 × 35 voxels) includes two dome-shaped epithelial surfaces of interest, forming the apical (inward) and basal (outward) sides of the organoid. C 3D representations of the height map extracted by Zellige (in green) and the GT height map (in blue), of the epithelium surface. D Color-coded error maps of the reconstructed height maps for the apical (left) and basal (right) epithelial surfaces of the organoid. The surfaces of interest are reconstructed with an error of < 2 pixels over a large majority (96% and 93% for the apical and basal surfaces, respectively) of pixels, as well as on average (RMSE ~ 0.8 and 1.1 for the apical and the basal surfaces, respectively). E Projections localized to the GT height maps of the epithelium surface (panels on the left), and the height maps extracted by Zellige (panels on the right). F Using the DeProj tool (Herbert et al. [12]) to generate cell morphology measurements corrected for the projection factor associated with the local surface slope. Cells of the surface S2, which represents the epithelial apical surface, were segmented using the TissueAnalyzer tool (Etournay et al. [3]). The height map obtained from Zellige and the segmentation mask were used to calculate geometrical parameters of the surface and its constituent cells. From left to right: the local slope, the local mean curvature, the projection error on cell area due to the local slope, and the corrected cell area in a 3D representation (axis ratio 1:1:3), for the reconstructed surface S2. Scale bar 100 μm

Full size image

Figure 6C shows a 3D graphical representation of the height maps reconstructed by Zellige and those reconstructed by an expert biologist, taken as ground truth (GT). Due to their dome-shaped topography, the manual segmentation of these surfaces was rather laborious and is more likely prone to errors in the regions of high inclination. Figure 6D shows that the distance between the two height maps is <2 for the large majority (> 92%) of pixels, while larger error values occur locally in regions of near-vertical inclination of the surfaces. The RMSEs of both apical and basal surfaces are close to 1, showing that on average the reconstructed height maps match the corresponding GT with about 1 pixel of accuracy. The coverage of the reconstructed height map is nearly 100% for the basal surface, and ~88% for the apical surface. Figure 6E shows the 2D projections of the 3D image obtained for each of the reconstructed surfaces, as compared to the projected GT height maps. Thus, for the extraction for these highly inclined surfaces, Zellige produces height maps of a quality comparable to those obtained from manual expert segmentation.

In this example, Zellige could be combined with a 2D cell tracking framework such as TissueMiner [3] to perform cell shape and topology analysis. Note that, using complementary tools such as DeProj [12], geometric distortions introduced in the projected surfaces by the epithelium inclination could be corrected for (Fig. 6F). This approach could provide a means to quantitatively address how an inner ear organoid epithelium patterns at the cellular and organoid scales, while quantifying the epithelial thickness changes due to cellular intercalation or cell shape changes in the depth of the epithelium. This would also permit to better characterize the variability of inner ear organoids within a given aggregate, and it could allow one to explore how the organoid interacts with surrounding tissues and how these interactions influence the differentiation of their constituent sensory cells.

A sensitivity analysis reveals the robustness of Zellige in extracting surfaces from biological images

To evaluate further the quality and robustness of the segmentation obtained by Zellige, we carried out a sensitivity analysis of the reconstruction on each of the samples tested. This analysis consisted in varying one control parameter at a time (Additional files 3, 4, 5, 6 and 7: Figure S2-S6, Fig. 7, Additional file 1: Supplementary Note 4), while keeping the other parameters fixed at a nominal value (Additional file 1: Supplementary Table S1). The RMSE and coverage of each of the reconstructed surfaces were evaluated and plotted as a function of the value of the parameter that was varied.

Fig. 7
figure 7

Summary of the sensitivity analysis. The intervals indicated in gray for each parameter and each of the images tested correspond to the parameter values for which the reconstruction satisfies high-quality criteria defined by RMSE ≤ 1.5 and coverage ≥ 85%. Black marks indicate the reference value obtained by manual adjustment for each image (cf. Additional file 1: Supplementary Note 2). A Parameters of the surface selection step. B Parameters of the surface assembly step

Full size image

Additional file 4: Figure S3 shows the results of the sensitivity analysis carried out on the image of Fig. 3 (Example 1, pupal fly specimen), when varying the parameters controlling step 1, i.e., the selection of putative surface voxels (parameters TA, Totsu, Smin, σxy, and σz). As can be seen, the variations of the two classification threshold parameters TA, Totsu, and of the minimum island size Smin in their respective intervals do not substantially modify the RMSE and the coverage of the reconstructed surfaces, whose values remain roughly constant for TA ≤ 8, Totsu ≤ 12, and over the entire Smin interval (Additional file 4: Figure S3A). Within these intervals, the RMSEs of all the reconstructed surfaces remain ≤1.5, while the coverage values are >95% for surfaces S1 and S2, and >85% for surfaces S3 and S4. The surfaces S1 and S3 show lower signal intensity and lower contrast than surfaces S2 and S3, making them more difficult to extract. Surface S4 has the lowest contrast and fails to be reconstructed if the classification threshold values are too stringent (namely for Totsu > 12, or TA > 8). Nevertheless, the intervals of stability of TA and Totsu (that is, the intervals of values over which a high-quality extraction of all surfaces is obtained) remain relatively wide (cf. Fig. 7). The smoothing parameters σxy and σz also have some effect on the reconstruction of the surfaces. When σxy is less than 3, the RMSE is higher for the mesh-like epithelial surfaces S2 and S4 (formed by the junctional network of the epithelium). A minimal smoothing along the axial direction is also important to ensure that the reconstructed surfaces are not too fragmented, preventing their complete reconstruction. Yet, σz should not be chosen too large either, to avoid merging nearby surfaces along the z-axis. In this case, the closely positioned surfaces S1 and S2 are well-separated if setting σz close to 1, but they become merged when σz > 2. In general, the surface construction parameters have little effect on RMSE and coverage (Additional file 4: Figure S3B). The sensitivity analysis on this challenging specimen shows a good robustness of Zellige to extract the four surfaces with a single set of parameters, each of which can be chosen in a reasonably wide interval considering the other fixed.

The results of the sensitivity analyses performed with the other biological image stacks (examples 2, 3, and 4 described above) are shown in Additional files 5-7: Figures S4-S6 (Additional file 1: Supplementary Note 4). These results are summarized in Fig. 7, which show the stability intervals over which the extracted height maps satisfy the criteria RMSE ≤ 1.5 and coverage ≥ 85%, for which the reconstruction can be considered of high quality. Overall, the stability intervals for the two classification threshold parameters (TA and Totsu) are narrower for specimens containing a surface of low signal intensity and low contrast relatively to other structure of higher signal, as in examples 2 and 3. The graphical user interface of Zellige allows the user to adjust TA and Totsu interactively, making it intuitive to search for reasonable values. We found that 2 ≤ σxy ≤ 3 and σz = 1 generally give high-quality results for all tested specimens (Fig. 7A). We therefore expect only little adjustment to be required by the user on the smoothing parameters from their default permissive values (set to σxy = 1 and σz = 1). Values of σz that are too large may lead to the merging of surfaces with a nearby structure of high contrast (surface or else), as it happens for the epithelium surface in the cochlea specimen, which merges with the underlying mesenchyme signal when σz is greater than 2 (Fig. 4 and Additional file 5: Figure S4). The effect is even more pronounced for closely positioned surfaces (Figs. 2 and 6 and Additional files 3, 7: Figures S2, S6).

Regarding the control parameters of step 2 (the assembly step), for the four examples presented except for example 1, the quality of the reconstructed surfaces is stable and high over the main part of their respective intervals of variation, deteriorating occasionally only when extreme values are used for these parameters (see Additional files 3, 4, 5, 6 and 7: Figures S2-S6, Fig. 7B and Additional file 1: Supplementary Note 4). Example 1 poses particularly stringent constraints on the control parameters of the reconstruction due to the requirement of reconstructing the 4 surfaces of different contrast and texture present in this sample using the same set of parameters.

Finally, the computation times for running Zellige on a given dataset ranged between a few seconds and a minute on a standard PC computer (see Additional file 8: Figure S7 and Additional file 1: Supplementary Note 5), except in a few exceptional cases corresponding to extreme values of the control parameters (not shown). As a safeguard, a stopping criterion could be implemented so as to exit the run (declaring the current parameter values invalid) if the surface assembly computation exceeds a user-prescribed duration.

Overall, the sensitivity analysis indicates that the surface extraction performed by Zellige is robust to variations of the control parameters of step 1 (surface voxel selection step). In general, the reconstruction is more sensitive to the amplitude threshold parameter (TA), and the Otsu threshold parameter (Totsu) should be kept sufficiently low for samples containing surfaces of intensity close to the background. However, in some cases such as in our example 3, the opposite is true. Thus, the two threshold parameters play somewhat complementary roles, and the possibility to adjust them independently is useful in practice to be able to cover as many cases as possible. A smoothing along xy appears necessary to correctly reconstruct the surfaces supported by a junctional mesh. Not surprisingly, best results are obtained when the radius of the Gaussian filter used for this (parameter σxy) is adapted to the mesh (or cell) size. Likewise, a smoothing along z is beneficial, but the extent of this smoothing (parameter σz) should not be too large to avoid causing the fusion of nearby surfaces. With a few exceptions, the values of the RMSE and coverage show little sensitivity to the values of the parameters controlling step 2 (surface assembly step), at least once putative surface voxels have been properly selected. In the presence of several surfaces of potentially very different sizes, the parameter controlling the fraction of OSE sizes allowed for OSE seeds (parameter TOSE1) should be relatively large (≥0.5 or greater, i.e., allowing more than 50% of the largest OSE sizes for seeds) to allow the extraction of a surface of small size (for example, to extract the surface S4 of example 1, which covers less than 20% of the xy-field of view, TOSE1 must be larger than 0.6). Extreme values (close to 1) for the connectivity rates (C1 and C2) are too stringent and lead to a drop in the coverage of the reconstructed surfaces. To sum up, we see that the most critical parameters for a satisfactory extraction of the different surfaces are those controlling step 1. In most cases, the parameters controlling step 2 do not need to be adjusted and can be fixed to their default reference values.

Conclusion

We have developed Zellige, a new tool to extract multiple surfaces from 3D fluorescence microscopy images. Zellige automatically finds surfaces by first identifying putative voxels that are likely to belong to a biological surface, and second by assembling a surface through connection of adjacent voxels satisfying natural proximity constraints. By using Zellige on synthetic epithelium images, we have shown that it accurately reconstructs a surface with excellent performances in terms of both the distance to the ground truth height map and the surface coverage (Fig. 2). Zellige can deal with complex images containing multiple surfaces, with computation times not exceeding a few tens of seconds on a standard computer. Importantly, the user is not required to specify the number of surfaces to be extracted. In the Drosophila specimen (Fig. 3), the software readily extracts the 4 surfaces of interest that could be identified. Since Zellige detects putative surface voxels in the first step by combining local and global thresholds, it can discriminate between a surface of low intensity and some nearby surface or other structures of higher contrast, as is the case in the mouse cochlear embryo (Fig. 4). With this difficult dataset, we could also confirm Zellige’s robustness against very low signal-to-noise levels. The constructive approach of surface region growing used by Zellige in its second step enables it to circumvent the surface smoothness requirement that is classically assumed by other surface extraction tools. For instance, it could reconstruct the highly irregular surface of a bronchial tissue infected by SARS-CoV-2 (Fig. 5).

The robustness and flexibility of Zellige come at a price, namely, the requirement to specify 12 parameters when running the surface extraction. However, the sensitivity analysis we performed shows that adjusting only 4 of these parameters is enough in practice to handle a wide range of image types. These parameters correspond to intuitive notions (e.g., thresholding and smoothing), which makes Zellige particularly easy to use. The Fiji interface that we implemented to perform this adjustment should make Zellige even more user-friendly and effective for biological applications.

To our knowledge, Zellige is the only open-source tool that can extract an unspecified number of epithelial surfaces from a 3D volume, possibly larger than two. This unique feature is especially useful in complex images that could be processed only by specialized tools before. For instance, Zellige can extract surfaces with projections on the xy plane that completely overlap, such as the basal and apical epithelia in the organoid image of Fig. 6. Previously, such an image could be processed only by tools that relied on segmenting a mesh around the object surface, such as MorphoGraphX or ImSAnE [17, 51].

The flexibility and robustness of Zellige should allow to considerably relax the constraints that were previously imposed on the sample preparation and the image acquisition steps by the subsequent analysis. Indeed, Zellige can accommodate any number of surfaces in the acquired volume, overlapping or not, and of different contrast features. Zellige also showed excellent robustness against image noise. This should make it particularly useful in imaging contexts that are not easily amenable to automated analysis, such as intravital imaging. Finally, it is worth noting that Zellige is a generalist and modular method. With some adaptation of the surface voxel selection step, it could be used with imaging modalities beyond the scope of this article, for instance, in extracting the irregular and noisy surfaces of biological objects imaged with 3D electron microscopy images.

Methods

Human bronchial epithelium imaging

The data used in Fig. 5 were taken from the recent study [30] to which we refer for the preparation and imaging of human bronchial epithelium cultures. In brief, MucilAirTM were purchased from Epithelix (Saint-Julien-en-Genevois, France) and cultured for at least 4 weeks to reconstruct a differentiated human bronchial epithelium in vitro and stained as previously described. Images of the cultures were acquired using an inverted Zeiss LSM 710 confocal microscope controlled by the ZEN pro 2.3 software. Z-stack images of whole-mount samples were acquired with a Zeiss Plan Apochromat ×63 oil immersion lens (NA=1.4). The image used here was published in Robinot et al. [30] under the CC-BY-4 license.

Drosophila imaging

Flies were raised at 25°C under standard conditions. Pupae were collected for imaging as described previously [52]. Ecad::GFP flies [53] were used for live imaging as previously described [1]. In brief, images were acquired with a spinning disk microscope from Gataca Systems driven by the MetaMorph software. The system is equipped with an inverted Nikon TI2E stand, a motorized XYZ stage, and a Nikon Plan Apo ×60 oil immersion (NA=1.4) lens and with a Prime95B camera.

Cochlea imaging

The inner ears from wild-type (C57BL/6) mice were rapidly dissected from temporal bones at embryonic stages E14.5 in HEPES-buffered (10 mM, pH 7.4) Hanks’ balanced salt solution and fixed in 4% paraformaldehyde, 1 h at room temperature. Specimens were permeabilized and stained for phalloidin-Atto 565 (Sigma) as previously described [44]. Fluorescence images were obtained with a swept-field confocal microscope (Opterra II) from Brucker. This system is equipped with a Nikon Plan Fluor ×60 oil immersion lens (NA=1.4).

Inner ear organoid imaging

ESCs derived from blastocyst-stage embryos of R1 mice (mESCs) (ATCC, SCRC-1036) were maintained in feeder-free culture on 0.1% w/v gelatin (Sigma) coated substrates using LIF-2i medium as established previously [54]. The organoids were generated following the previously published protocol [49, 54]. Aggregates were harvested at day 14 and fixed in 4% v/v PFA (Electron Microscopy Sciences) overnight at 4 °C. After blocking (PBS; 10% v/v normal goat serum; 0.1% v/v Triton X-100), the aggregates were stained for phalloidin Atto 565 (1:1000) (Sigma) overnight at RT on a shaker and washed three times with PBS containing 0.1% v/v Triton X-100 for 1 h each at RT. Prior to imaging, the aggregates were incubated in a modified version of ScaleS solution containing 4 M Urea (Sigma), 40% w/v D-Sorbitol (Sigma), and 0.1% v/v Triton X-100, for 3–5 days to clarify the tissue. Finally, the aggregates were whole-mounted using the ScaleS solution and imaged using a confocal laser scanning microscope (A1R HD25, Nikon) equipped with a Nikon ×25 silicon oil immersion lens (NA=1.05).

Availability of data and materials

The project homepage below contains the source code, installation instructions, and documentations.

Zellige

• Project name: Zellige.

• Project homepage: https://gitlab.pasteur.fr/ida-public/zellige-core

• URL for the Fiji plugin (Update > Manage update sites): https://sites.imagej.net/Zellige/

• Gitlab Branch: master

• Operating systems: Platform independent.

• Programming language: Java.

• Compiled in Java8

• Other requirements : Runs from Fiji [29].

• License: BSD 2

• Any restrictions to use by non-academics: License needed.

Scripts to create Phantoms

• Project name: Phantoms.

• Project homepage / Archived version: https://doi.org/10.5281/zenodo.6414596

• Operating systems: Platform independent.

• Programming language: MATLAB.

• License: BSD 2

• Any restrictions to use by non-academics: License needed.

Data sets

Image 3D stacks, ground truth height map, and height maps produced by Zellige are available on Zenodo under the CC-BY license: https://zenodo.org/communities/zellige/ [33, 39, 47, 48, 50].

References

  1. Etournay R, Popović M, Merkel M, Nandi A, Blasse C, Aigouy B, et al. Interplay of cell dynamics and epithelial tension during morphogenesis of the drosophila pupal wing. eLife. 2015;4:e07090. https://doi.org/10.7554/eLife.07090.

    Article  PubMed  PubMed Central  Google Scholar 

  2. Guirao B, Rigaud SU, Bosveld F, Bailles A, López-Gay J, Ishihara S, et al. Unified quantitative characterization of epithelial tissue development. eLife. 2015;4:e08519. https://doi.org/10.7554/eLife.08519.

    Article  PubMed  PubMed Central  Google Scholar 

  3. Etournay R, Merkel M, Popović M, Brandl H, Dye NA, Aigouy B, et al. TissueMiner: a multiscale analysis toolkit to quantify how cellular processes create tissue dynamics. eLife. 2016;5. https://doi.org/10.7554/eLife.14334.

  4. Saadaoui M, Rocancourt D, Roussel J, Corson F, Gros J. A tensile ring drives tissue flows to shape the gastrulating amniote embryo. Science. 2020;367:453–8. https://doi.org/10.1126/science.aaw1965.

    Article  CAS  PubMed  Google Scholar 

  5. Butler LC, Blanchard GB, Kabla AJ, Lawrence NJ, Welchman DP, Mahadevan L, et al. Cell shape changes indicate a role for extrinsic tensile forces in drosophila germ-band extension. Nat Cell Biol. 2009;11:859–64. https://doi.org/10.1038/ncb1894.

    Article  CAS  PubMed  Google Scholar 

  6. Collinet C, Rauzi M, Lenne P-F, Lecuit T. Local and tissue-scale forces drive oriented junction growth during tissue extension. Nat Cell Biol. 2015;17:1247–58. https://doi.org/10.1038/ncb3226.

    Article  CAS  PubMed  Google Scholar 

  7. Kong D, Wolf F, Großhans J. Forces directing germ-band extension in drosophila embryos. Mech Dev. 2017;144:11–22. https://doi.org/10.1016/j.mod.2016.12.001.

    Article  CAS  PubMed  Google Scholar 

  8. Villedieu A, Bosveld F, Bellaïche Y. Mechanical induction and competence in epithelial morphogenesis. Curr Opin Genet Dev. 2020;63:36–44. https://doi.org/10.1016/j.gde.2020.03.003.

    Article  CAS  PubMed  Google Scholar 

  9. Prat-Rojo C, Pouille P-A, Buceta J, Martin-Blanco E. Mechanical coordination is sufficient to promote tissue replacement during metamorphosis in drosophila. EMBO J. 2020;39:e103594. https://doi.org/10.15252/embj.2019103594.

    Article  CAS  PubMed  Google Scholar 

  10. Blasse C, Saalfeld S, Etournay R, Sagner A, Eaton S, Myers EW. PreMosa: extracting 2D surfaces from 3D microscopy mosaics. Bioinforma Oxf Engl. 2017. https://doi.org/10.1093/bioinformatics/btx195.

  11. Basu S, Rexhepaj E, Spassky N, Genovesio A, Paulsen RR, Shihavuddin A. FastSME: faster and smoother manifold extraction from 3D stack. In: 2018 IEEE/CVF conference on computer vision and pattern recognition workshops (CVPRW); 2018. p. 2362–23628. https://doi.org/10.1109/CVPRW.2018.00305.

    Chapter  Google Scholar 

  12. Herbert S, Valon L, Mancini L, Dray N, Caldarelli P, Gros J, et al. LocalZProjector and DeProj: a toolbox for local 2D projection and accurate morphometrics of large 3D microscopy images. BMC Biol. 2021;19:136. https://doi.org/10.1186/s12915-021-01037-w.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  13. Li K, Wu X, Chen DZ, Sonka M. Optimal surface segmentation in volumetric images-a graph-theoretic approach. IEEE Trans Pattern Anal Mach Intell. 2006;28:119–34. https://doi.org/10.1109/TPAMI.2006.19.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  14. Lombardot B. Min cost Z surface projection Fiji. Java. 2017; https://imagej.net/Minimum_Cost_Z_surface_Projection.

  15. Hall M, Frank E, Holmes G, Pfahringer B, Reutemann P, Witten IH. The WEKA data mining software: an update. ACM SIGKDD Explor Newsl. 2009;11:10–8. https://doi.org/10.1145/1656274.1656278.

    Article  Google Scholar 

  16. Sommer C, Straehle C, Kothe U, Hamprecht FA. Ilastik: interactive learning and segmentation toolkit. In: 2011 IEEE international symposium on biomedical imaging: from Nano to macro; 2011. p. 230–3. https://doi.org/10.1109/ISBI.2011.5872394.

    Chapter  Google Scholar 

  17. Heemskerk I, Streichan SJ. Tissue cartography: compressing bio-image data by dimensional reduction. Nat Methods. 2015;12:1139–42. https://doi.org/10.1038/nmeth.3648.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  18. Haertter D, Wang X, Fogerson S, Ramkumar N, Crawford J, Poss K, et al. DeepProjection: rapid and structure-specific projections of tissue sheets embedded in 3D microscopy stacks using deep learning; 2021. https://doi.org/10.1101/2021.11.17.468809.

    Book  Google Scholar 

  19. Green JBA, Sharpe J. Positional information and reaction-diffusion: two big ideas in developmental biology combine. Development. 2015;142:1203–11. https://doi.org/10.1242/dev.114991.

    Article  CAS  PubMed  Google Scholar 

  20. Miller CJ, Davidson LA. The interplay between cell signalling and mechanics in developmental processes. Nat Rev Genet. 2013;14:733–44. https://doi.org/10.1038/nrg3513.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  21. Levin M, Martyniuk CJ. The bioelectric code: an ancient computational medium for dynamic control of growth and form. Biosystems. 2018;164:76–93. https://doi.org/10.1016/j.biosystems.2017.08.009.

    Article  PubMed  Google Scholar 

  22. Duclut C, Prost J, Jülicher F. Hydraulic and electric control of cell spheroids. Proc Natl Acad Sci. 2021;118. https://doi.org/10.1073/pnas.2021972118.

  23. Ruiz-Herrero T, Alessandri K, Gurchenkov BV, Nassoy P, Mahadevan L. Organ size control via hydraulically gated oscillations. Development. 2017;144:4422–7. https://doi.org/10.1242/dev.153056.

    Article  CAS  PubMed  Google Scholar 

  24. Dumortier JG, Le Verge-Serandour M, Tortorelli AF, Mielke A, de Plater L, Turlier H, et al. Hydraulic fracturing and active coarsening position the lumen of the mouse blastocyst. Science. 2019;365:465–8. https://doi.org/10.1126/science.aaw7709.

    Article  CAS  PubMed  Google Scholar 

  25. Mammoto T, Mammoto A, Ingber DE. Mechanobiology and developmental control. Annu Rev Cell Dev Biol. 2013;29:27–61. https://doi.org/10.1146/annurev-cellbio-101512-122340.

    Article  CAS  PubMed  Google Scholar 

  26. Ladoux B, Mège R-M. Mechanobiology of collective cell behaviours. Nat Rev Mol Cell Biol. 2017;18:743–57. https://doi.org/10.1038/nrm.2017.98.

    Article  CAS  PubMed  Google Scholar 

  27. Närhi K. Embryonic explant culture: studying effects of regulatory molecules on gene expression in craniofacial tissues. In: Seymour GJ, Cullinan MP, Heng NCK, editors. Oral biology: molecular techniques and applications. New York: Springer; 2017. p. 367–80. https://doi.org/10.1007/978-1-4939-6685-1_21.

    Chapter  Google Scholar 

  28. Zeevaert K, Elsafi Mabrouk MH, Wagner W, Goetzke R. Cell mechanics in embryoid bodies. Cells. 2020;9:2270. https://doi.org/10.3390/cells9102270.

    Article  CAS  PubMed Central  Google Scholar 

  29. Schindelin J, Arganda-Carreras I, Frise E, Kaynig V, Longair M, Pietzsch T, et al. Fiji: an open-source platform for biological-image analysis. Nat Methods. 2012;9:676–82. https://doi.org/10.1038/nmeth.2019.

    Article  CAS  PubMed  Google Scholar 

  30. Robinot R, Hubert M, de Melo GD, Lazarini F, Bruel T, Smith N, et al. SARS-CoV-2 infection induces the dedifferentiation of multiciliated cells and impairs mucociliary clearance. Nat Commun. 2021;12:4354. https://doi.org/10.1038/s41467-021-24521-x.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  31. Pietzsch T, Preibisch S, Tomančák P, Saalfeld S. ImgLib2—generic image processing in Java. Bioinformatics. 2012;28:3009–11. https://doi.org/10.1093/bioinformatics/bts543.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  32. Otsu N. A threshold selection method from gray-level histograms. IEEE Trans Syst Man Cybern. 1979;9:62–6. https://doi.org/10.1109/TSMC.1979.4310076.

    Article  Google Scholar 

  33. Boutet de Monvel J, Trébeau C, Altay G, Tinevez J-Y, Etournay R. Zellige example dataset: synthetic image dataset. 2022. https://doi.org/10.5281/zenodo.6376584.

  34. Boutet de Monvel J. Scripts to create phantom images containg surfaces. Zenodo. 2022. https://doi.org/10.5281/zenodo.6414596.

  35. de Chaumont F, Dallongeville S, Chenouard N, Hervé N, Pop S, Provoost T, et al. Icy: an open bioimage informatics platform for extended reproducible research. Nat Methods. 2012;9:690–6. https://doi.org/10.1038/nmeth.2075.

    Article  CAS  PubMed  Google Scholar 

  36. Jülicher F, Eaton S. Emergence of tissue shape changes from collective cell behaviours. Semin Cell Dev Biol. 2017;67:103–12. https://doi.org/10.1016/j.semcdb.2017.04.004.

    Article  PubMed  Google Scholar 

  37. Morata G, Lawrence P. An exciting period of drosophila developmental biology: of imaginal discs, clones, compartments, parasegments and homeotic genes. Dev Biol. 2022;484:12–21. https://doi.org/10.1016/j.ydbio.2022.01.008.

    Article  CAS  PubMed  Google Scholar 

  38. Ray RP, Matamoro-Vidal A, Ribeiro PS, Tapon N, Houle D, Salazar-Ciudad I, et al. Patterned anchorage to the apical extracellular matrix defines tissue shape in the developing appendages of drosophila. Dev Cell. 2015;34:310–22. https://doi.org/10.1016/j.devcel.2015.06.019.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  39. Trébeau C, Boutet de Monvel J, Altay G, Tinevez J-Y, Etournay R. Zellige example dataset: drosophila pupal wing and abdomen; 2022. https://doi.org/10.5281/zenodo.6376566.

    Book  Google Scholar 

  40. Zallen JA. Planar polarity and tissue morphogenesis. Cell. 2007;129:1051–63. https://doi.org/10.1016/j.cell.2007.05.050.

    Article  CAS  PubMed  Google Scholar 

  41. Montcouquiol M, Rachel RA, Lanford PJ, Copeland NG, Jenkins NA, Kelley MW. Identification of Vangl2 and Scrb1 as planar polarity genes in mammals. Nature. 2003;423:173–7. https://doi.org/10.1038/nature01618.

    Article  CAS  PubMed  Google Scholar 

  42. Cohen R, Amir-Zilberstein L, Hersch M, Woland S, Loza O, Taiber S, et al. Mechanical forces drive ordered patterning of hair cells in the mammalian inner ear. Nat Commun. 2020;11:5137. https://doi.org/10.1038/s41467-020-18894-8.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  43. McKenzie E, Krupin A, Kelley MW. Cellular growth and rearrangement during the development of the mammalian organ of Corti. Dev Dyn. 2004;229:802–12. https://doi.org/10.1002/dvdy.10500.

    Article  CAS  PubMed  Google Scholar 

  44. Etournay R, Lepelletier L, Boutet de Monvel J, Michel V, Cayet N, Leibovici M, et al. Cochlear outer hair cells undergo an apical circumference remodeling constrained by the hair bundle shape. Dev Camb Engl. 2010;137:1373–83. https://doi.org/10.1242/dev.045138.

    Article  CAS  Google Scholar 

  45. Driver EC, Northrop A, Kelley MW. Cell migration, intercalation and growth regulate mammalian cochlear extension. Dev Camb Engl. 2017;144:3766–76. https://doi.org/10.1242/dev.151761.

    Article  CAS  Google Scholar 

  46. Qian D, Jones C, Rzadzinska A, Mark S, Zhang X, Steel KP, et al. Wnt5a functions in planar cell polarity regulation in mice. Dev Biol. 2007;306:121–33. https://doi.org/10.1016/j.ydbio.2007.03.011.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  47. Trébeau C, Boutet de Monvel J, Altay G, Tinevez J-Y, Etournay R. Zellige example dataset: mouse embryonic cochlea stained for F-actin; 2022. https://doi.org/10.5281/zenodo.6376542.

    Book  Google Scholar 

  48. Robinot R, Chakrabarti L, Michel V, Trébeau C, Boutet de Monvel J, Altay G, et al. Zellige example dataset: primary culture of human bronchial cells infected by SARS-CoV-2; 2022. https://doi.org/10.5281/zenodo.6376594.

    Book  Google Scholar 

  49. Koehler KR, Mikosz AM, Molosh AI, Patel D, Hashino E. Generation of inner ear sensory epithelia from pluripotent stem cells in 3D culture. Nature. 2013;500:217–21. https://doi.org/10.1038/nature12298.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  50. Altay G, Trébeau C, Boutet de Monvel J, Tinevez J-Y, Etournay R. Zellige example dataset: inner ear organoid epithelium; 2022. https://doi.org/10.5281/zenodo.6376582.

    Book  Google Scholar 

  51. Barbier de Reuille P, Routier-Kierzkowska A-L, Kierzkowski D, Bassel GW, Schüpbach T, Tauriello G, et al. MorphoGraphX: a platform for quantifying morphogenesis in 4D. eLife. 2015;4:e05864. https://doi.org/10.7554/eLife.05864.

    Article  PubMed Central  Google Scholar 

  52. Classen A-K, Aigouy B, Giangrande A, Eaton S. Imaging drosophila Pupal wing morphogenesis. In: Dahmann C, editor. Drosophila: methods and protocols. Totowa: Humana Press; 2008. p. 265–75. https://doi.org/10.1007/978-1-59745-583-1_16.

    Chapter  Google Scholar 

  53. Huang J, Zhou W, Dong W, Watson AM, Hong Y. Directed, efficient, and versatile modifications of the drosophila genome by genomic engineering. Proc Natl Acad Sci. 2009;106:8284–9. https://doi.org/10.1073/pnas.0900641106.

    Article  PubMed  PubMed Central  Google Scholar 

  54. Nie J, Koehler KR, Hashino E. Directed differentiation of mouse embryonic stem cells into inner ear sensory epithelia in 3D culture. In: Organ regeneration. New York: Humana Press; 2017. p. 67–83. https://doi.org/10.1007/978-1-4939-6949-4_6.

    Chapter  Google Scholar 

Download references

Acknowledgements

We thank Maia Brünstein of the Hearing Institute Bioimaging Core Facility - C2RT/C2RA and Florian Rückerl of the Photonic BioImaging - C2RT/C2RA, for sharing their expertise on light microscopy. We thank Romain Levayer for sharing flies and lab space with us. We thank Gilles Trébeau (https://fr.linkedin.com/in/gilles-trebeau-427601212) for helpful discussions regarding the Java frontend deployment of Zellige. We thank Lisa Chakrabarti, Rémy Robinot, and Vincent Michel for sharing data.

Funding

This work was supported by Fondation pour l’Audition (FPA-IDA-STARTING-GRANT); Institut Pasteur (PTR#272); DIM ELICIT’s grant from Région Ile-de-France (QP-IDF-DIM-ELICIT-2019); The French Agence Nationale de la Recherche (ANR-21-CE13-0038-013-01 and LabEx LIFESENSES ANR-10-LABX-65).

Author information

Authors and Affiliations

  1. Institut Pasteur, Université Paris Cité, Inserm, Institut de l’Audition, F-75012, Paris, France

    Céline Trébeau, Jacques Boutet de Monvel, Gizem Altay & Raphaël Etournay

  2. Institut Pasteur, Université Paris Cité, Image Analysis Hub, F-75015, Paris, France

    Jean-Yves Tinevez

Authors
  1. Céline Trébeau
    View author publications

    You can also search for this author in PubMed Google Scholar

  2. Jacques Boutet de Monvel
    View author publications

    You can also search for this author in PubMed Google Scholar

  3. Gizem Altay
    View author publications

    You can also search for this author in PubMed Google Scholar

  4. Jean-Yves Tinevez
    View author publications

    You can also search for this author in PubMed Google Scholar

  5. Raphaël Etournay
    View author publications

    You can also search for this author in PubMed Google Scholar

Contributions

CT designed the program and wrote the Java code. JYT supervised the Fiji plugin development and reviewed the Java code. JBM wrote the MATLAB scripts and generated all synthetic images. RE and GA performed some of the experiments presented here. CT, JMB, JYT, and RE jointly wrote the article. All authors participated in regular group discussions to develop the ideas presented here, which are the result of multiple cycles of tests on synthetic and biological data. The authors read and approved the final manuscript.

Authors’ information

Twitter handles: @jytinevez (Jean-Yves Tinevez)Twitter handles are transferred to the "Authors' information" section., @institutpasteur (the Institut Pasteur).

Corresponding authors

Correspondence to Jean-Yves Tinevez or Raphaël Etournay.

Ethics declarations

Ethics approval and consent to participate

RE holds a designer certificate of animal experimentation (level 2), allowing to perform experimental work on animals in strict accordance with the European directive 2010/60/EU, and French regulations. The Ethics Committee of the Institut Pasteur (Comité d’Ethique en Experimentation Animale - CETEA) has approved this study with the project identifier dha170006. This approval is based on careful compliance to the 3Rs principle in the care and use of animals (Annex IV - 2010/60/EU).

Consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

Additional information

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary information

Additional file 1: Supplementary Notes 1-5.

Note 1. Overview of Zellige implementation. Note 2. Generation of phantom 3D images. Note 3. Comparing the reconstructed and ground truth height maps. Note 4. Sensitivity Analysis. Note 5. Computational/processing time of Zellige. Table S1. Zellige parameters.

Additional file 2: Figure S1.

Zellige implementation. (A) The two main algorithmic steps of Zellige. Upper part: Surface voxel selection step (step 1). Lower part: Surface assembly step (step 2). (B) Determination of the amplitude A(p) of some local maximum p along the z-axis, relative to its two closest local minima p1 and p2. (C,D) Connectivity rules used to connect putative surface voxels together for the construction of orthogonal surface elements (OSEs). These rules are based on a 6-connectivity relationship defined within each orthogonal (xz or yz) section as shown in (C). This relationship is extended to allow the presence of single straight gaps, while being constrained to forbid the occurrence of forking points. The allowed local neighborhood configurations along an OSE are illustrated in (D) in the case of a construction within xz sections. The OSEs are formally defined as the connected components of the graph GOSE defined by these rules within each orthogonal section. (E) Compatibility rules used in the surface assembly step, illustrated here in the case where OSEs have been constructed within xz sections, and assembly proceeds along the y-axis. To validate the addition of a new OSE σ constructed within section y+1 (in blue) to the surface S under construction, whose intersection with section y is shown (surface line ly, in green), two quantities are computed: the overlap R(S,σ) is defined as the number of pixels of σ that share the x coordinate of some pixel of ly. The connectivity C(S,σ) is defined as the fraction of the overlapping pixels of σ that are 6-connected, in their respective yz section, to the corresponding point of ly (according to the 6-connectivity relationship shown in C, middle panel). In the depicted example, there are 2 OSEs (in green) with, for both of them, R(S,σ)=3 and C(S,σ)=3/3=1. The OSE will be added to S if R(S,σ) ≥ R0 and C(S,σ) ≥ C0, where R0 and C0 are tunable thresholds that set the stringency of the matching condition.

Additional file 3: Figure S2.

Sensitivity analysis of Zellige on the phantom image. (A) Surface voxel selection parameters. (B) Surface assembly parameters. Reference values are indicated by the dashed line.

Additional file 4: Figure S3.

Sensitivity analysis of Zellige on the pupal fly image. (A) Surface voxel selection parameters. (B) Surface assembly parameters. Reference values are indicated by the dashed line. (C) RMSE and coverage results obtained for different values of the σxy blur taken within and outside the interval defined by the parameter sweep in (A, “σxy” panel). (D) Visualization of the projections on the extracted surfaces obtained for the different σxy blur values shown in (C).

Additional file 5: Figure S4.

Sensitivity analysis of Zellige on the cochlear epithelium image. (A) Surface voxel selection parameters. (B) Surface assembly parameters. Reference values are indicated by the dashed line.

Additional file 6: Figure S5.

Sensitivity analysis of Zellige on the primary epithelium culture image. (A) Surface voxel selection parameters. (B) Surface assembly parameters. Reference values are indicated by the dashed line.

Additional file 7: Figure S6.

Sensitivity analysis of Zellige on the inner ear organoid image. (A) Surface voxel selection parameters. (B) Surface assembly parameters. Reference values are indicated by the dashed line.

Additional file 8: Figure S7.

Computational time analysis of Zellige. (A) Computation times (step 1 in blue, step 2 in brown) for processing the various images tested on a PC notebook computer with (processor Intel Core i9 2,4 GHz with 32 Gb of Ram). Apart for the case of a highly rough surface (culture specimen) the computation time is largely dominated by the surface voxel selection step (step 1). (B) The same computation times are re-plotted as a function of image size N (number of voxels). Note the linear growth of the computational time of step 1 as a function of N, while that of step 2 shows much slower growth.

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 licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence 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 licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and Permissions

About this article

Verify currency and authenticity via CrossMark

Cite this article

Trébeau, C., de Monvel, J.B., Altay, G. et al. Extracting multiple surfaces from 3D microscopy images in complex biological tissues with the Zellige software tool. BMC Biol 20, 183 (2022). https://doi.org/10.1186/s12915-022-01378-0

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s12915-022-01378-0

Keywords

  • Image analysis
  • Morphology
  • 3D imaging
  • Tissue imaging
  • Image segmentation
  • Surface extraction
  • Fiji plugin

BMC Biology

ISSN: 1741-7007

Contact us