1.

INTRODUCTION

Biological soft tissues are functional substrates at the intersection between the cell and the organ-organism levels. From a functional perspective, they play critical roles at both scales: On one hand, they constrain how neighbouring individual cells can interact with each other. On the other hand, their woven structure provides emerging mechanical characteristics that become a critical part of their roles, e.g. defining essential properties in smooth, skeletal or cardiac muscles. When looking at the brain, histology becomes an integral part in the study of neuronal circuits – ensembles of interconnected neurons constrained in space that perform a certain computation. However, while tissues have long been studied at the macroscopic [1] and at the microscopic level [2], a multiscale study of biological tissues has not been technologically feasible until recently: resolving subcellular features across 1 mm³ requires leveraging 3D image datasets in the terabyte range.

Automation of volume electron microscopy [3-6] has unlocked the potential to bridge those two scales and already provided insights on how different tissues operate [7-9]. More recently, other nondestructive approaches have shown the capacity to operate in this regime, either employing light [10, 11] or X-rays [12, 13] as illumination source. In particular, X-ray holographic nanotomography (XNH) [14] has proven capacity to resolve subcellular features across landscapes that extend to the mm scale in samples with preserved ultrastructure and can be compatible with follow-up volume electron microscopy [12, 15, 16].

In turn, automated image segmentation algorithms have seen an unprecedented improvement in recent years and high-performance compute clusters have become widely accessible to the materials and life sciences domains [17, 18]. Pioneering efforts in volume electron microscopy have provided essential datasets that facilitated training the first generation of connectome-delivering classifiers [19-21].

However, while obtaining mm-scale volume datasets at subcellular resolution is becoming a routine approach, instance classifiers that generalise well across independently sourced datasets of the same tissue are not common, resulting in a bottleneck in the data analysis pipeline. Here we aim to address this bottleneck by developing a nuclei classifier for mouse olfactory bulb tissue.

Cells are the essential components of tissues, and therefore the distribution and shape of their nuclei defines essential properties of the tissue itself. But their manual annotation is unfeasible already at the mm-scale: 1 mm³ of mouse brain tissue contains on the order of 2*10⁵ nuclei and a 0.5 cm³ mouse brain contains 100*10⁶ [20, 22]. Moreover, cell nuclei display a characteristic set of morphological traits that make them good targets for automated image classifiers: they have a compact shape, with volumes within a narrow range, and the ultrastructure of their interior - the nucleoplasm - often displays a distinct pattern to the adjacent cytosol of the cell that contains it. Therefore, classifiers aimed at segmenting nuclei can be trained to recognise not only the thin nuclear membrane but also the nucleoplasm.

While several tools have been developed to detect subcellular features in volume EM images [23-25], specific solutions to automatically retrieve nuclei across mm-scale datasets acquired by XNH are not yet developed. Those classifiers will provide unique insights on the typical counts, sizes and distributions of cells within the context of tissues.

Here, we developed and trained an image classifier that reliably segments cell nuclei in mouse brain tissue volumes acquired with XNH. It operates integrated on a high performance computing cluster, providing straightforward scalability to larger volumes within the terabyte scale. It detects the shapes, volumes and positions of cell nuclei in tissue volumes from the mouse olfactory bulb with high precision and recall across multiple cell layers and generalises well to analogous datasets. This makes it an essential tool integrated in a broader image analysis pipeline devoted to mammalian systems neuroscience addressing how neuronal circuits compute information given their connectivity [16, 26]. Furthermore, it is also applicable to the study of the structure-function of other biological tissues.

2.

RESULTS

We first aimed to automatically detect nuclei in a 3D volume of mouse brain olfactory bulb acquired with X-ray holographic nanotomography. The dataset spanned across 0.456 mm³ with voxels of (100 nm)³, containing a continuous tissue volume of 0.212 mm³ and was obtained after a non-rigid stitching of 28 individual tiles (Fig. 1a, Table 1) [27]. We developed an image segmentation pipeline that provided an automated volumetric labelling for each nucleus inside a 0.197 mm³ mask that covered the vast majority of the tissue and avoided non-tissue areas (e.g. those regions occupied by bare resin or empty space around the sample) (Fig. 1b,c), computed individual meshes for each nucleus (Fig. 1d) and returned an indexed table with essential morphometric parameters of each segmented item: its volume, its surface and its sphericity (Fig. 1d) [28]. This segmentation pipeline consists of two major steps in sequence (Fig. 1e): an automated instance segmentation relying on a trained U-Net architecture [29] followed by a heuristic filtering of segmented instances based on the expected volumes and sphericities of nuclei.

Figure 1.

Workflow to segment cell nuclei in tissue volumes obtained with X-ray holographic nanotomography.

(a) XNH can provide continuous volumes mapping mm-scale landscapes with subcellular resolution. (b-c) Nuclei can be detected and segmented automatically, providing a quantitative insight on their morphologies and distribution in space (d). (e) The classifier used consists of a trained U-Net followed by a filter based on prior knowledge on the morphology of the targeted features. The arrowheads indicate correctly segmented nuclei (cyan) and segmentation artefacts discarded by the filtering step (yellow). GL, glomerular layer; EPL, external plexiform layer; MCL, mitral cell layer; GCL, granule cell layer.

The volume mapped in this dataset contains all histological layers in the olfactory bulb - each displaying a characteristic pattern of nuclei in terms of their individual size as well as in their spatial arrangement (density and distribution). We defined subvolumes where all nuclei contained were manually annotated at a 2-fold downscaled version of the datasets: (200 nm)³ voxels. A total of 30 subvolumes sampled 4 histological regions of interest: glomerular layer (GL, 10 subvolumes), external plexiform layer (EPL, 5 subvolumes), mitral cell layer (MCL, 10 subvolumes) and granule cell layer (GCL, 5 subvolumes) (Fig. 2a, b). Each of those subvolumes was (50 μm)³ in size to ensure large 20 μm-wide nuclei would be fully contained in the training data even in the regions where their density is the lowest, and therefore minimise split errors of the classifier. We make all the training and test data openly available for further reuse (see Data Availability section). The segmentation workflow consists of a 3D U-Net [29] to obtain a distance prediction, followed by a watershed algorithm [30] to generate an instance segmentation of cell nuclei (Fig. 2c). Due to the large size of the datasets, the prediction and segmentation was performed in smaller chunks which were merged afterwards.

Nuclei were positively identified in most cases (Fig. 2d, f) and rarely missed. However, certain features not constituting nuclei were also often erroneously classified as such (Fig. 2e). Manual inspection of the false positive segments revealed that these were features of diverse nature but with morphologies very distinct from nuclei. An unbiased analysis of detection accuracy confirmed this insight: truly detected nuclei clustered together in a parameter space defined by the volume and sphericity of the segmented features (Fig. 3e). An optimal double-threshold boundary therefore allowed filtering segmented nuclei from the pool of all segmented features: we defined both high-pass thresholds (volume and sphericity) by applying a grid-search in 2-fold cross-validation over 12 subvolumes of size (50 μm)³ (Methods). Together, the filtered result of the classifier provided a final segmentation for this dataset, consisting of 64166 nuclei (Fig. 3g) with

a precision and recall of above 90% across most regions (Fig. 3f). Each detected nucleus was stored with an associated identifier, a computed mesh and a summary of its key morphological characteristics (volume, surface and sphericity) (Fig. 3g, h, i).

Figure 2.

Architecture of the classifier

(a) Training data was obtained by manually annotating all nuclei present in 30 regions of interest each covering (50 μm)³, distributed across all relevant histological layers of the olfactory bulb tissue in the dataset. (b) Examples of manually annotated nuclei in the four layers mapped. (c) Diagram of the U-Net [29] and following segmentation steps. (d-f) This classifier provided a good result detecting nuclei (d), while some voxels not belonging to nuclei were sometimes also classified as such (e). Similarly, while most non-nuclear features were correctly ignored (f) in very rare occasions true nuclei were missed. (g) The 3D segmentations were later converted to standard meshes and stored as such, associated with the original dataset. (h) The measured sphericity of closed meshes extracted from segmented objects matched expected theoretical values for diverse regular polygon shapes across a broad range of sizes. GL, glomerular layer; EPL, external plexiform layer; MCL, mitral cell layer; GCL, granule cell layer.

Systems neuroscience - and tissue biology more generally - challenges hypotheses through experiments whose readouts are often vulnerable to sampling biases and sample preparation procedures. This vulnerability can be sometimes overcome by incorporating sufficient biological replicates in the experimental design. Two more independent datasets reporting continuous tissue volumes of 0.269 mm³ and 0.297 mm³ of the same brain region were acquired in a similar way with XNH, providing a cohort of 3 comparable datasets (Fig. 3d, j, p, Table 1). We applied the segmentation pipeline without modifications to the second dataset (Y391, 0.261 mm³ mask covering most tissue) (Fig. 3j, Table 1) and computed new dataset-specific thresholds for volume and sphericity using the same method as for C432 (see Methods) (Fig. 3k, Table 1). This provided an automated segmentation of 76820 nuclei with a precision and recall above 75% in most regions (Fig. 3l, m).

In order to provide a generalisable approach for further datasets, nuclei in the third dataset (Y489) (Fig. 3p, Table 1) were segmented using a battery of different classifiers incorporating new additional local ground truth (Fig. 3c). A total of 29 subvolumes each covering (50 μm)³ and distributed across layers of interest (10 GL, 5 EPL, 10 MCL and 4 GCL) were manually annotated for nuclei (Fig. 3a-b). Multiple classifiers were trained using different sets of training data: from the unmodified classifier to versions adding groups of 5 additional subvolumes sampling different histological layers from the Y489 dataset, to a classifier only trained with the new 29 subvolumes. All classifiers were used to independently detect all nuclei in the same Y489 dataset (Fig. 3c). All classifiers showed high precision (>95%) when evaluated on a test set of 14 subvolumes of (50 μm)³. While the unmodified classifier on the new dataset showed only a 90% recall, addition of 5 subvolumes to the subset of training data raised the overall recall rate to >95%. The GL and GCL are the two most demanding tissue layers for a nuclear classifier: there, cell nuclei are very densely packed, often with little cytoplasmic volume in cells whose cell bodies sit directly adjacent to each other (Fig. 3b). The amount of features available to avoid misclassification are therefore minimal - often only the nuclear and plasma membranes (as opposed to the histological layers EPL and MCL, where cytoplasmic content and neurites of the same and other neurons are often present between different nuclei). When observing the performance of those classifiers at the most stringent regions, the addition of all 29 new subvolumes with locally sourced new training data did not harm performance and in some cases it would seem to improve the segmentation of densely packed nuclear clusters in the GCL. Therefore, the final classifier for Y489 used was the one trained on all previous and new ground truth data (59 subvolumes (7.35*106 μm³) of training data). It provided an automated segmentation of 81779 nuclei (Fig. 3s), with a precision above 95% and recall above 90% across most regions (Fig. 3r).

Table 1.

Compute time used to run the instance segmentation pipeline and post-processing on each dataset at (200 nm)3 voxels (ie, 2-fold downscaling applied to original data). The classifier pipeline was run on a Slurm cluster on several CPU nodes equipped with an AMD EPYC 7763 CPU (hyperthreading enabled) of which we used up to 12 cores and 50 GB RAM each. The mentioned times include queuing times for sub-tasks.

Figure 3.

Generalisation to other datasets of the same imaging modality and similar characteristics

(a) Training data was also obtained from a third independently acquired dataset from a new sample that mapped the same brain region. (b) Examples of manually annotated nuclei in the four layers targeted. (c) Precision, recall and F1 score of the classifier operating on the new dataset when trained with increasing amounts of data sourced from the new dataset. (d) Volume rendering of the first dataset. (e) Objects segmented, distributed according to each item’s volume and sphericity. Prior knowledge allowed describing thresholds for both metrics, providing a filtered result that detects nuclei with high accuracy. (f) Precision and recall on the filtered results at multiple regions of the dataset, including high-resolution tile centers (ctr) as well as tile border regions (border) from all histological regions. (g) Volume rendering of the meshes of all nuclei segmented from this dataset. (h-i) Volumes (h) and sphericities (i) of the segmented nuclei in this dataset. (j-o) Same as d-i for the second dataset, segmented with the same classifier without any additional training data. (p-u) Same as d-i for the third dataset, using a classifier that incorporates all the newly sourced training data (performance shown in (c), 30+29 subvolumes). GL, glomerular layer; EPL, external plexiform layer; MCL, mitral cell layer.

3.

DISCUSSION

We have developed an image classifier that accurately segments cell nuclei in 3D maps of mouse brain tissue obtained with X-ray holographic nanotomography. This classifier not only provides satisfactory results within the dataset it has been trained on, but also generalises its performance to analogous datasets obtained independently with minimal addition of locally sourced training data. In this line, we observed a drop in recall when segmenting the Y391 dataset, which could be attributed to two factors: poorer signal-to-noise ratio (SNR) of the dataset itself (attributable to being reconstructed using tomograms taken at three instead of four distances) and no addition of locally sourced training data. In dataset Y489, where both factors were amended, we obtained F1 scores >95%. When comparing its performance within the same dataset, across regions with high SNR (tile center) and low SNR (tile border), as expected better performance is achieved where SNR is best while providing in all cases a sufficient yield for the intended analyses. Therefore, this classifier will be a useful component of an integrated analytical pipeline devoted to analysing mm-scale 3D maps of tissue ultrastructure in the mouse brain.

It operates on datasets that map the dense histology of samples with preserved ultrastructure: the contrast in the images is driven by the original distribution of proteins and lipids in space. Its performance on analogous independent datasets suggests that the classifier is robust enough to cope with variations of signal and noise patterns in the data. One would therefore expect this classifier to provide acceptable results on volumetric datasets reporting tissue ultrastructure obtained with any imaging technique that shares the essential features with XNH, such as the data being aligned within < 100nm precision in all directions, having isotropic voxels, and the contrast being driven by the distribution in space of all proteins and lipids, with an additional set of locally sourced training data comparable to the one reported in our third dataset. While focussed ion beam scanning electron microscopy and ptychographic X-ray tomography would most likely fit these requirements, other volumeEM techniques such as serial block-face electron microscopy could also be eligible provided that the quality of image registration meets the required precision. On the other hand, its performance on datasets obtained with imaging techniques that target serial sections (such as ssTEM, gridTape or ATUM-mSEM) will have to be assessed independently since these will incorporate slicing-specific artefacts such as cracks and folds that have not been included in the training data.

This classifier structure proved extremely efficient and robust to detect nuclei in XNH datasets. We expect other features mapped in these datasets - such as volume occupied by tissue, resin or blood vessels - to be detectable by specifically-trained classifiers with a similar architecture, possibly requiring adjustments related to the biological nature of the features. Moreover, promising classifiers that segment mitochondria in vEM datasets [24] could be leveraged and incorporated in the analytical pipeline to further broaden the insights one can quickly gather automatically from newly acquired XNH datasets.

However, segmentation of features in tissue is fundamentally limited by resolution: biological tissues are composed by cells and organelles, and both features are delimited by a 20 nm-thick bi-layered lipidic membrane. While resolving this essential 20 nm feature enables segmenting any cell or subcellular component [7, 31-33], relevant merging errors should be expected when aiming to segment organelles with coarser resolutions.

The analysed areas of the three datasets reported cover mouse brain olfactory bulb tissue volumes of 0.197 μm³, 0.261 μm³ and 0.282 μm³, and contain 64166, 76820 and 81779 nuclei with an average size of 212.9±0.5 μm³, 207.9±18.9 μm³ and 217.7±0.4 μm³ and and an average sphericity of 0.8114±2*10^-4, 0.7871±3*10^-4 and 0.7849±2*10^-4, respectively (mean ± SEM). Addressing these annotations manually would have been unfeasible: allowing 10 seconds to find, annotate and proof-read each individual nucleus would involve ~200 human-hours per dataset only to find their locations. Fully segmenting them in 3D (to extract their volumes and sphericities), allowing 24 min/nucleus, would take ~30*10³ human-hours per dataset. Recent developments in X-ray imaging and in volumeEM are boosting the capacity to generate these datasets at scale – enabling experimental timeframes compatible not only with the acquisition of singular proof-of-concept datasets but of controlled dataset cohorts like the one reported in this study. This technological capacity is already delivering novel insights on the function of biological soft tissues, led largely by the field of connectomics, focussed on the study of neuronal circuits. In mammals, neuronal circuits are embedded in volumes not smaller than 1 mm³ [34], and a whole mouse brain at connectomic resolution is regarded as a feasible aim within a decade [35, 36]. Processing the datasets that will be acquired in the upcoming years will necessarily require a robust data architecture that uses computational resources efficiently. This classifier ran at an average rate of 387 min per processed 0.1 mm³ (Table 1). We anticipate it being operative at this regime for up to 1 mm³ datasets. For larger volumes, there are opportunities for optimization, such as reducing scheduling overheads on the Slurm cluster and utilising GPUs for predictions. With these optimizations, we expect processing volumes of 10 mm³ (equivalent to the entire olfactory bulb of an adult mouse) to take a few days on a similarly equipped Slurm cluster.

Exploring, annotating and proof-reading terabyte-scale datasets requires a data architecture that ensures robustness of the original dataset over long periods of time, that facilitates sharing and access by many users from diverse locations worldwide, that streams efficiently the data users require using the minimal bandwidth necessary, that integrates well with high-performance image analysis computational routines, and that relies on an open-access data format ensuring long-term support for both datasets and metadata. OME-Zarr has been suggested as the consensus data format that can meet all these requirements: this chunked data format allows any user to explore a petascale dataset from their browser, it is optimised for distributed computing routines addressed at high performance clusters and is open-source [37, 38]. We employed a Zarr-compatible dataset architecture, WEBKNOSSOS [39] to store, explore and annotate our datasets. We release all the densely annotated training data and tagged test data in this format, enabling us and others to develop more accurate segmentation algorithms in the future.

Altogether, we have developed an image classifier that accurately segments nuclei in 3D mouse brain datasets obtained with X-ray holographic nanotomography. This classifier generalises well to other datasets of independently processed and imaged samples with minimal additional training data and operates natively on Zarr-compatible datasets stored in a WEBKNOSSOS architecture. It makes it an integral tool of a scalable image analysis pipeline for future datasets exploring mouse brain tissue with X-ray holographic nanotomography at the mm³ scale. This classifier will provide biologically relevant insights on terabyte-scale volume datasets mapping biological tissues early in the analytical pipeline, and it will pave the way to establish robust image segmentation solutions for even larger datasets.

4.

MATERIALS AND METHODS

Segmentation post-processing

We computed the segment volumes and center positions chunk-wise with subsequent aggregation. As positions are non-aggregatable, we enlarged each chunk by 20μm at each side to enable determining the positions for segments spanning multiple chunks. Segment volumes were calculated via voxel counting, segment center positions via Euclidean distance transform.

Mesh calculation

Aside from the meshes for visualisations (see Rendered visualisations), we generated meshes within our post-processing procedure with the marching cubes algorithm [51]. The resulting meshes consisted of connected vertices and the enclosed triangular faces (Fig. 2g). We applied mesh simplification using an implementation of the quadric mesh reduction algorithm [52] from the Python package pyfqmr to reduce the number of triangles by a factor of 100 while restricting the simplification error to a maximum of 16000 nm. The meshes were computed and stored in chunks, merging was performed when loading a mesh.

Sphericity

We computed the surface in each mesh by summing up the areas of all the faces of the mesh. Then, for each mesh of a segmented object, we computed its sphericity - a metric of compactness, reporting the ratio of the surface-to-volume ratio of a 3D object to the one of a sphere of the same volume [28]. Accordingly, perfectly spherical objects will have a sphericity of 1, whereas highly eccentric objects will display sphericity values close to 0. Sphericity can be obtained from the particle’s volume and surface:

Note that by generating meshes with our mesh calculation process, mesh vertices do not match exactly the borders of voxels. This mismatch leads to slight inaccuracies in the calculated sphericities because we computed the surface from the meshes and the volume as voxel count from the segments. However, those inaccuracies become only relevant for small (≤3 voxels-wide) 3D features which can then generate artifactual sphericity values larger than 1. Since nuclei are expected to present overall round structures and their diameter (10 μm) is typically >100 times larger than the voxel size (100 nm), their characteristic sphericities can be reliably retrieved (Fig. 2h).

To determine the volume and sphericity thresholds for C432 and Y391, we conducted 2-fold cross-validation over a set of 12 subvolumes spanning the different histological regions. We obtained the thresholds via grid search over the volumes and sphericities of the detected objects and selected the combination that resulted in the best F1 score. We valued higher recall over higher precision. Thus, we restricted the recall computed over all training boxes combined to be at least 0.8. Further, we performed aggregation over the folds by selecting the minimum for each threshold. For Y489, we applied the thresholds computed on C432.

Evaluation

As described in (Sphericity), we calculated the thresholds for C432 and Y391 on 12 boxes. We further kept a held-out test set of 6 additional boxes. On Y489, we used 14 boxes for classifier selection and 6 boxes as test set. All boxes covered (50 μm)³ and contained skeleton annotations of nuclei instances in the MCL, GL and EPL layer - three boxes in the centre and three boxes in the border region of each layer. For C432 and Y391, the boxes were evenly distributed among the cross-validation folds and the test set, each containing one box per region and layer. The test set for Y489 followed the same principle, the classifier selection set contained two boxes per region and layer with addition of one box from a darker region of each the MCL and the GL layer.

When calculating precision, recall, and F1 score per box, we aimed to minimise measurement artefacts due to nuclei being cut off at the box border. Therefore, we enlarged each box by 8μm in each direction beforehand, which was treated as a tolerance margin that contained no ground truth annotations but predicted nuclei in the margin were not counted as false positives. The average precision and recall values in Table 1 reflect the average score across subvolumes (Fig. 3f,l,r).