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