A first (and intuitive) attempt for the data-driven clustering of spinal tracts consisted in inputting a 2D axial slice (corresponding to one spinal level) for each of the five morphometrics. Hence 31 clusterings were performed independently (one for each spinal level).

Only voxels inside the white matter mask were considered for the clustering. Each morphometric map was normalized between 0 and 1 before entering the clustering.

Clustering of the morphometrics data was performed using scikit-learn’s agglomerative clustering algorithm. The number of clusters was varied between 5 and 11 (with an increment of 1) to explore the impact of this parameter on the generated clusters. A connectivity matrix was defined to enforce pixel-to-pixel connection within the axial plane. Clusters were merged by minimizing the variance of the Euclidean distance between cluster candidates (linkage = ward, affinity = Euclidean).

While the slice-wise approach described above seemed the most faithful to the data, it produced unsatisfactory results (see results section), which we attributed to noise and to the relatively low number of samples. To obtain a more robust clustering, we enforced additional spatial priors, namely: we aggregated slices within each of the four spinal regions (cervical, thoracic, lumbar and sacral), and we non-linearly co-registered slices together within the same region to be able to define the connectivity matrix along the superior-inferior axis in addition to the axial plane. The approach is described hereafter.

Building on the anatomical knowledge that the spatial distribution of tracts is relatively similar within a given region (e.g., cervical, thoracic), we ran the clustering by aggregating adjacent axial slices within a given region. But before doing that, to account for the slight displacement of funiculi across slices, we co-registered each axial slice to the slice located at the center of a region. This co-registration was performed using ANTs’ bspline-regularized SyN transformation on the axon density metric map. To ensure a smooth and robust transformation, instead of directly registering each slice to the center slice, we instead performed a step-by-step registration where the center-slice r is registered with slice r + 1, then slice r + 1 is registered with slice r + 2, and so on. Transformations were then concatenated and applied to the corresponding slice. This algorithm is illustrated in Figure 3. Following co-registration, morphometric data were averaged along the superior-inferior axis, producing 2D slices for each region (n = 4) and each morphometric (n = 5), so 20 in total.

For convenient comparison with the Watson atlas, the estimated warping fields were also applied to the digitized Watson atlas. Then, as done for the morphometrics data, the region-specific registered Watson slices were averaged together, resulting in a single 2D axial slice per region. Technical note: each tract was defined in a 4th dimension, so that averaging across the 3rd dimension did not result in an “overlap” of the tracts, i.e., each z-averaged tract ranged from 0 to 1 in order to estimate partial volume information, which was subsequently used to produce colored maps of the atlas for qualitative assessment.

A quantitative validation of the clustering results was tricky to implement because there was no clear way to justify what cluster to consider for computing a distance (e.g., Hausdorff) or overlap (e.g., Dice) metric between the clusters and tracts from the Watson atlas. Therefore, we opted for a qualitative (visual) evaluation, as done in Assaf et al. (2008). Moreover, the qualitative evaluation was augmented by quantitatively selected colormaps and intensity: we generated 2D images showing the clusters next to the Watson atlas, with matching color-coding based on the level of overlap. When multiple clusters were found within a given tract, the intensity level would be varied between 0.2 and 1, scaled with the amount of overlap with the Watson tract.

Figure 4 shows, for each spinal level, the clustering results (left hemi-section) and the Watson atlas (right hemi-section). The color coding was topologically matched for convenient comparison. Here, the number of clusters was set to 8, which corresponds to the number of tracts in the Watson atlas. We observe a somewhat good topological correspondence between the generated clusters and that from the Watson atlas, notably for the cuneate fasciculus (red) and the dorsal corticospinal tract (cyan). Some smaller tracts (gracile, post synaptic dorsal column, lateral spinal nucleus, rubrospinal) are not well captured by the clustering. The ventral funiculus, which is a large collection of ventral tracts from the Watson atlas, appears parcellated on the generated clusters, and this parcellation is not consistent across slices. This lack of consistency is possibly related to noise from the histology data and to the ill-posed problem of clustering.

The inconsistencies observed with the slice-wise clustering were tackled by aggregating slices within the same region (cervical, thoracic, lumbar, sacral). To achieve this, we first non-linearly registered slices pertaining to the same group. Figure 5 shows satisfactory results of this inter-slice registration for the cervical region. Registration went equally well for the other three regions.

The clustering algorithm was applied using the registered slices as inputs for all the metric maps. Figure 6 shows results of the clustering for n = 8 and n = 10 clusters. Like for the slice-wise clustering, the generated clusters are shown on the left section of the spinal cord while the Watson atlas is shown on the right section. Overall, the same two clusters (cuneate fasciculus, dorsal corticospinal tract) are clearly delineated across all four regions, with a particularly good topological correspondence for the cervical and thoracic regions. The Gracile Fasciculus (black) and the Post Synaptic Dorsal Column (green) appear to be merged into the same cluster (light red) for the cervical and thoracic regions. This merging could be caused by similar morphometrics features between these two tracts, that the automatic clustering was not able to distinguish. The sacral slice seems a bit “noisy,” which could be due to the fact that tracts are smaller in the lower spinal cord region.

To get a sense for how the fixed number of clusters affect the results, Figure 7 shows clustering results when the number of clusters was varied from 7 to 30. This result is presented for a representative slice (C7 level). Visually, a better correspondence is observed with a higher number of clusters. More specifically, the gracile fasciculus (black) and the lateral spinal nucleus (yellow) start to appear more consistently as the number of clusters increase. However, even with a relatively high number of clusters, the rubrospinal tract (purple) is hard to distinguish. Interestingly, at 13 clusters and higher, the dorsal corticospinal tract (cyan) is horizontally split between two clusters, suggesting two different fiber populations within this tract.

The slice-wise clustering generated two clusters that were relatively stable across slices and that corresponded well to the topology of the Watson atlas, namely the cuneate fasciculus and the dorsal corticospinal tract (See Figure 4). Other (smaller) clusters were less consistent. We hypothesized that this lack of consistency was mainly due to the noise in the histology maps, and that more robust clustering would be obtained by aggregating more slices together (see additional considerations in the discussion below). The region-wise clustering did indeed produce clusters that were more consistent across regions, and with more clusters matching the topology of the Watson atlas. Despite that, some smaller tracts (gracile, postsynaptic dorsal column, lateral spinal nucleus, rubrospinal) were not well captured by the clustering. The ventral funiculus, which is a large collection of ventral tracts from the Watson atlas, appears parcellated on the generated clusters, and this parcellation is not consistent across slices.

This lack of consistency is possibly related to noise from the histology data and/or to the ill-posed nature of clustering (Rajapakshage and Pensky, 2020). Indeed, clustering assumes that (i) a solution exists, (ii) the solution is unique and (iii) the data provide sufficient information to reach the solution. In this project, a core assumption was that tracts could be identified from each other only based on their morphometric signature. This immediately questions the choice of morphometrics that were selected for this clustering work. Here, we included the following metrics: axon density, axon diameter, axon volume fraction, g-ratio and myelin thickness. This decision was based on the qualitative and quantitative observation of the variability of each of these metrics across spinal tracts. However, it is possible that this choice was suboptimal, and that clustering would benefit from removing some metrics or adding others. For example, the g-ratio is relatively flat across tracts, so it could potentially be removed from the analysis (hence reducing the noise). To investigate the effect of selecting less metric maps for the clustering, we rerun the slicewise clustering with only the axon density and the axon volume fraction maps (see results in Supplementary Figure 2). Interestingly, the clustering result is fairly similar to that when using all five metric maps, suggesting that some of the metrics used in the initial clustering analysis are mostly contributing. Further work could focus on the optimal choice of axon/myelin morphometrics for clustering spinal tracts.

Another (strong) assumption is that each tract is considered ‘homogeneous’, i.e., the axons that compose them are assumed to all have the same morphometrics (e.g., axon diameter of 800 nm, myelin thickness of 500 nm). We know very well this assumption is wrong, because axons that compose a given white matter tract exhibit a relatively wide range of morphometrics (Saliani et al., 2017, 2019; Duval et al., 2019). For example, Saliani et al. (2019) showed that the distributions of axonal density clearly overlap across most tracts, as illustrated in Figure 8. Interestingly, the dorsal corticospinal tract has a very different distribution from the other tracts (much denser axons), which likely explains the relatively good stability of that cluster in the present study. Similarly, axon diameter distributions also overlap across multiple tracts. The implication, from a clustering point-of-view, is that it makes the morphometric-based separation of individual tracts particularly challenging, because of the ill-posed nature of the problem.

The decision to perform region-wise clustering was twofold. Firstly, without the region-wise aggregation, clustering results appeared somewhat noisy and unreliable across adjacent slices, hence motivating the aggregation of pixels sharing similar anatomical features, even though we acknowledge this may be an oversimplification in some parts of the anatomy, especially in small tracts which shape and location is changing along the rostro-caudal axis. The medium and bigger tracts, however, show some levels of spatial overlap within a given region as can be qualitatively observed in the Watson atlas. Secondly, it has been quite common in the neuroanatomy community to represent the topology of the spinal cord white matter within the cervical, thoracic, lumbar and sacral. For example, the Gray’s anatomy atlas of the human spinal cord only shows one level per region (Standring, 2008). The work of Nathan et al. (1990), which focuses on the corticospinal tract in humans, does report that “the arrangement [of the lateral CST] is more constant throughout the thoracic segments,” although, to be fair, more topological variability is observed within the cervical spinal cord (notably around the cervical enlargement). The idea behind using non-linear bspline-regularized registration based on the shape of the gray matter and spinal cord, however, is somewhat supported, in the same article by Nathan et al., by the observation that the “location of the lateral corticospinal tract is influenced by the size and position of the anterior and posterior horns as well as by the shape of the cord itself.” However, this observation only concerns the CST, and other tracts might be displaced differently than the CST. To sum up, while we acknowledge that the proposed within-region clustering is an oversimplification, especially for smaller tracts, the slowly varying topology of the white matter tracts along the rostro-caudal direction is not a new idea, and while the region-wise aggregation is possibly too aggressive, the proposed methodology could be adapted to apply the rostro-caudal regularization between fewer slices (e.g., 3–5 adjacent slices, instead of within a full region).

The clustering algorithm asked for a fixed number of clusters to be generated. We started with 8, based on the number of tracts identified in the Watson atlas. Increasing the number of clusters (up to 30) produced better delineation of some tracts that were otherwise not captured with a smaller number of clusters. This is not surprising, given that more clusters enable adjacent regions with subtle morphometric differences to be separated. However, fixing a high number of clusters also defeats the purpose of automatically finding one and only one cluster associated with an individual tract.

The previous considerations naturally bring up the question of validation. As mentioned in the methods section, a quantitative validation of the clustering results would have been convoluted due to the exploratory nature of the present work. For example, computing a spatial similarity metric between the clusters and tracts from the Watson atlas would be inherently biased toward the arbitrary pairing between a given cluster and a given tract. As we observed, some tracts are represented by several clusters (especially in the ventral funiculi), therefore it is unclear how to best compute a spatial similarity in those instances. Instead, we opted for a visual interpretation of the produced clusters, aided by a color-coding and intensity leveling computed from the overlap between the produced clusters and each of the tracts from the Watson atlas.

Future work could also explore inputting spatial prior to the clustering algorithm for making the clusters more stable and reliable. For example, spatial priors based on the existing Watson atlas (and/or other atlases) would enforce the identification of tracts that are difficult to identify solely based on their morphometrics.

There exist different histology techniques with their pros and cons. Optical microscopy can provide a full coverage of the slices but at a limited resolution (∼200 nm). The other extreme, transmit electron microscopy, can provide nanometric resolution, but only on small spatial windows (about 50 μm × 50 μm) taken sparsely across the tissue, hence it does not provide a full picture. Blockface transmit electron microscopy can mitigate this covering issue, but at the cost of a more complex and expensive imaging setup. Scanning electron microscopy falls in the middle, with the possibility to achieve 50–100 nm resolution while covering the entire tissue. In this project we opted for scanning electron microscopy. While we were able to obtain morphometric maps of the entire spinal cord cross section across all spinal levels, the resolution was somewhat limiting in our ability to precisely and accurately delineate axon and myelin tissue, biasing quantities such as myelin thickness and g-ratio (Saliani et al., 2019). Also, it is very likely that very small myelinated axons (with internal diameter smaller than 200–300 nm) were missed by the segmentation algorithm and hence not accounted for in all metrics. Moreover, histology sections can be hampered by all sorts of artifacts such as improper fixation (inducing degradation of the myelin sheath), poor penetration of osmium, improper polishing, intensity bias across scanning electron microscopy sub-images and bad focus (Cohen-Adad, 2018).

The segmentation of individual axons and surrounding myelin sheath originated from a previous study, which relied on the AxonSeg software (Zaimi et al., 2016). We acknowledge the limitation of this fully automated histology segmentation tool, where a non-negligible amount of axons were poorly segmented (false negatives, false positives, “leaky” segmentation or sub-segmentation). However, the number of segmented instances is very high: with about 250 axons per 2500 μm² (see Figure 8), and about 50,000 voxels in the white matter of the rat-atlas (Saliani et al., 2019), it sums up to 125,000,000 segmented axons and the same amount of segmented myelin, which realistically cannot be manually verified/corrected. While the AxonSeg software has since been superseded by more performant methods based on deep learning (Zaimi et al., 2018), every fully automatic has some levels of failures that will unfortunately bias the produced morphometric maps.

In addition to segmentation issues, other processing steps could have introduced undesired bias. Notably, the various filtering, inter-rat co-registration, downsampling to 50 μm resolution, and pixel-wise metric aggregation methods performed when building the metric maps which is described in Saliani et al. (2019), the co-registration between adjacent slices for the within-region clustering, and the right-left symmetrization to maximize shared information during clustering. In particular, the choice of 50 μm resolution for the metric maps (Saliani et al., 2019) results in a compromise between having sufficient axon count to obtain a reliable axon and myelin volume fraction estimates, and having sufficiently small pixel sizes to be able to resolve fine anatomical features. In the context of this clustering work, where the goal was to retrieve white matter bundles that are at least several hundreds of μm² (see Saliani et al., 2017 for a review on axon morphometry in the rat spinal cord), this working resolution was likely sufficient. However, for other work looking at subtle variations within bundles, finer resolution might be desirable.

The ground truth used to validate the clustering results was the rat atlas of the spinal cord (Watson et al., 2009) (a.k.a. the “Watson” atlas), which is not without flaws. Firstly, this atlas was created from a single rat, hence does not represent the possible inter-specimen variability. Secondly, the white matter tracts were identified and manually drawn, based on multiple evidence from previous studies that relied on cellular injection of staining agents (Kayalioglu, 2009; Watson and Harvey, 2009). Moreover, some tracts are missing, as stated by the authors:

The “missing” tracts from the Watson atlas include the ventral and lateral spinothalamic tracts and the spinocerebellar tracts (ascending), the medial and lateral vestibulospinal tracts and the reticulospinal tract (descending). It is possible that our clustering method did identify some of these tracts (or at least in part), although without a proper ground truth one can only speculate. Another – possibly more interesting – approach, would be to turn the problem around, and instead of aiming at validating the clustering method based on hard-to-obtain cellular tracking techniques, the idea would be to expand our knowledge of the white matter tracts distribution by combining the traditional with the cluster-based techniques. There is an undeniable lack of knowledge about how white matter tracts are organized and distributed in the spinal cord across species. Including histology-driven clustering is potentially a useful technique to increase this knowledge. The open-source analysis pipeline published with this study could be useful to further develop and study white matter anatomy in rats and other species. Although this pipeline was applied to microscopic images, it can equally be used on other imaging modalities, including MRI, which offers quantitative measures of axons and myelin (Seiberlich et al., 2020).