 Research
 Open access
 Published:
A robust approach to 3D neuron shape representation for quantification and classification
BMC Bioinformatics volume 24, Article number: 366 (2023)
Abstract
We consider the problem of finding an accurate representation of neuron shapes, extracting subcellular features, and classifying neurons based on neuron shapes. In neuroscience research, the skeleton representation is often used as a compact and abstract representation of neuron shapes. However, existing methods are limited to getting and analyzing “curve” skeletons which can only be applied for tubular shapes. This paper presents a 3D neuron morphology analysis method for more general and complex neuron shapes. First, we introduce the concept of skeleton mesh to represent general neuron shapes and propose a novel method for computing mesh representations from 3D surface point clouds. A skeleton graph is then obtained from skeleton mesh and is used to extract subcellular features. Finally, an unsupervised learning method is used to embed the skeleton graph for neuron classification. Extensive experiment results are provided and demonstrate the robustness of our method to analyze neuron morphology.
Introduction
The importance of neuronal morphology has been recognized from the early days of neuroscience [1]. There are three obstacles in automatic neuron morphology analysis. First, we need to have a good shape representation of each neuron. Skeleton representations are widely used in neuroscience [2,3,4,5] as they provide a compact and abstract shape representation. Mathematically, skeletonization or medial axis transform (MAT) has a rigorous definition for arbitrary shapes. The skeleton of a shape is defined as a collection of interior points that have at least two closest points on the surface of the shape. We refer to those interior points as skeleton points and each skeleton point is associated with a radius. Figure 1 shows an example of MAT. However, in reality, it is not an easy task to get skeleton representation directly from images. Most automatic or manual segmentation methods output a cloud of surface points. Thus, we need to compute the 3D neuron skeleton from 3D surface point clouds. The skeleton representation further enables computing subcellular features such as length and number of branches of neurons as well as classification of neurons.
The main contribution of this paper is a robust and efficient method for computing a skeleton representation from a set of 3D surface points. This 3D skeleton representation can be used for a quantitative analysis of neuronal cell structures, including subcellular feature calculations and for neuron type classification based on 3D shapes.
There is an extensive literature on neuron skeleton extraction [3, 6, 7]. In [3], the skeleton representation is computed from a 3D mesh by using a traditional morphological thinning algorithm [8]. This method has two main drawbacks. First, the thinning algorithm is sensitive to noise of 3D mesh. Second, in reality, we usually get discrete 3D surface points of neurons from the segmentation step, and constructing 3D mesh from those discrete 3D surface points will introduce additional noise. To make the skeleton extraction model more robust, Ref. [7] proposes to use deep learning network to learn skeleton representation. The main idea of the paper is to use the deep learning network to predict skeletons from features of multiple spatial scale layers. This model still takes a continuous surface as input, as opposed to discrete surface points. Further, this is a supervised method and it requires training samples. In [6], they propose extracting skeleton representations directly from discrete surface points by using a 3D discrete distance transform. However, this only works well for curve skeletons and only tubular structures have curve skeletons. General 3D shapes will result in surface skeletons as shown in Fig. 2. In practice, skeleton mesh is used to represent the surface skeleton. There is no existing method to extract skeleton mesh from surface point clouds for neuron morphology analysis.
There are also methods to analyze neuron shapes using the skeleton representation. For skeleton classification task, [3, 9,10,11,12,13] use predefined handcrafted features to represent neuron morphology and classify neurons based on those representations. [3] proposes to compare similarity of skeletons by using local skeleton features. It breaks a neuron skeleton into short segments and characterize segments by location and direction of segments. Similarly, compared to [3, 11] provides a more robust method to find branches of neuron which do not have many local fluctuations. Reference [9, 10] extract other cellular and subcellular features, such as length of neuron, the surface area of soma, dendritic length from “curve” skeletons to represent neuron morphology. Reference [12] models the skeleton as a graph and uses paths in the graph as the feature to represent neuron morphology. Reference [13] defines a specific descriptor function to capture global and local information from the skeleton. The main drawback of the handcrafted feature is its limited representation capability. To solve the problem of handcrafted feature, several learning based methods [14,15,16] have been proposed recently to classify neuron skeletons. References [14, 15] project 3D skeletons back to 2D images and use convolution neural networks to get skeleton representation from 2D images. However, 3D shape information is lost when projecting onto 2D images even with multiview projections like [14]. To avoid projecting 3D skeletons into 2D images, Ref. [16] models the neuron skeleton as a graph directly and proposes a contrastive graph neural network (GNN) learning framework to represent the neuron. Similar to all above mentioned methods, Ref. [16] only works for “curve” skeletons but not for surface skeletons. There are fundamental differences between curve skeletons and skeleton meshes. For example, there are no cycles in curve skeletons but that is not the case in skeleton meshes. Ref. [17] can handle skeleton meshes. However, they only consider features from skeleton points for the classification and it does not fully utilize the skeleton mesh information.
To analyze general neuron shapes, this paper presents a robust 3D neuron morphology analysis framework based on the surface skeleton representation of neurons. In [18], the authors propose an unsupervised deep learning skeleton mesh extraction method. However, this method does not work well when neurons have concave shapes. Our skeleton mesh extraction method is built upon [18], and by using estimated surface norm of point clouds as part of the optimization function, we address this drawback. Next the skeleton mesh is converted to an undirected graph called skeleton graph. Inspired by [19], we embed the skeleton graph by maximizing mutual information, and then classify neurons based on the embedding of each skeleton. To compute cellular/subcellular features of neurons from the skeleton representation, we also utilize the skeleton graph. A simple but effective recursive algorithm is proposed to get number of branches and length of neurons.
We apply our neuron morphology analysis method to classify Ciona neurons. The Ciona sea squirt is one of the widely studied tunicates in neuroscience [20]. Its brain is closely related to vertebrates with a much simpler neuronal structure. In a single Ciona larva, it has only about 187 neurons with about 6600 synapses [20]. Studying the Ciona brain in depth can reveal the general principles behind the mechanism of how vertebrate brains work [21]. We also present our results on the NeuroMorpho[22] dataset. In summary, the main contribution of our paper include

A robust and efficient skeleton mesh extraction method with novel cost function by using properties of MAT. To the best of our knowledge, this is the first one to use skeleton mesh instead of the curve skeleton to analyze neuronal shapes

A 3D Ciona neuron dataset that can be used for neuron morphology analysis.
Method
Figure 3 illustrates our overall neuron morphology analysis method. Given a set of surface point clouds as input, we introduce an unsupervised deep learning method to get the skeleton mesh representation of each neuron. This is achieved by using the properties of the traditional medial axis transform (MAT). The skeleton mesh representation includes skeleton points with radii as well as the connection of those skeleton points as shown in Fig. 3. Second, the skeleton representation of each neuron is transformed into a skeleton graph. Each node in the skeleton graph represents a skeleton point. If there is an edge between two nodes, it means those two skeleton points are connected. The weight of the edge represents distance between the two skeleton points. Radii as well as the location of each skeleton point are attributes of each node. Next, length and number of branches of neurons are computed from the skeleton graph. To compare different shapes of neurons, a graph level representation learning method is used to embed the skeleton graph. The representation learning method is an unsupervised method that maximizes mutual information of the skeleton graph.
Skeleton mesh from 3D surface point cloud
Our unsupervised 3D neuron skeleton extraction method is built upon the method in [18] as illustrated in Fig. 4. Given a 3D surface point cloud as input, PointNet++ [23] is used as the encoder to obtain the sampled surface points with features. Next, a multilayerperceptron (MLP) is used to learn the geometric transformation to predict the skeleton points with their radii using linear combination of the MLP input points. Compared to [18], we propose to use properties of skeleton points as the prior knowledge to make the geometric transformation learning process more robust to general shapes. After getting skeleton points with radii, a graph autoencoder is used to predict links between skeleton points.
Skeleton points prediction
Mathematically, given a set of 3D surface points \({\textbf {P}}\in \mathbb {R}^{M\times 3}\) where M is a number of points, we want to predict N skeleton spheres \({\textbf {s}}_i=<{\textbf {c}}_i,r_i>\) where \({\textbf {c}}_i\) is the center of each sphere and \(r_i\) is the radius of each sphere.
As illustrated in Fig. 4, we first use PointNet++ [23] as the encoder to obtain a set of sampled input points \({\textbf {P}}'\in \mathbb {R}^{M'\times 3}\) and their associated contextual features \({\textbf {F}}\in \mathbb {R}^{M'\times D}\). \(M'(M'<M)\) is the number of feature points after PointNet++ and D is the feature dimension of each feature point. PointNet++ groups points and extract point features hierarchically. It contains a number of set abstraction levels. For each set abstraction level, there are three layers: Sampling layer, Grouping layer, and PointNet layer. For the first set abstraction level, the input is \({\textbf {P}}\), a set of M number of 3D surface points. Next, Sampling layer is applied. The iterative farthest point sampling (FPS) [23] is used to get \(M'\) sampled points. In the grouping layer, those \(M'\) sampled points are used as the centroid points. Then for each centroid, all M points within a radius are viewed as neighbor points and are grouped into that local region. Therefore, each centroid has K neighbors and K can vary for different groups. After Sampling and Grouping layers, PointNet [24] layer is used to extract features for each local region. The sampling layer, grouping layer, and PointNet layer consists one set abstraction level, and we stack such abstraction levels to form a hierarchical architecture to get features at different spatial scales. Next, multiscale grouping is applied to concatenate the features from different spatial scales.
A MultiLayer Perceptron (MLP) is used to estimate the geometrical transformation to get a set of skeleton spheres’ center points \({\textbf {C}}\), \(\{ {\textbf {c}}_1,{\textbf {c}}_2,...,{\textbf {c}}_N \}\). The geometric transformation we use is a convex combination of input points \({\textbf {P}}'\). MLP with softmax function is used to estimate the weight \({\textbf {W}}\in \mathbb {R}^{M'\times N}\) of the convex combination in Eq. 1.
As shown in [18], the same weight matrix W can be used to compute \(r({\textbf {c}})\in {\textbf {R}}\) using Eq. 2
where \(D\in \mathbb {R}^{M'\times 1}\) is a vector of \(d({\textbf {p}}',{\textbf {C}})\). \(d({\textbf {p}}',{\textbf {C}})\) is the closest distance of one surface point \({\textbf {p}}'\) to all skeleton points \({\textbf {C}}\) and is defined as \(d({\textbf {p}}',{\textbf {C}})=\min _{{\textbf {c}}\in {\textbf {C}}}{\textbf {p}}'{\textbf {c}}_2\)
A set of loss functions are defined in [18] to train the MLP. The loss function includes sampling loss \(L_s\), pointtosphere loss \(L_p\), and radius regularizer loss \(L_r\). The first two losses are based on the recoverability of skeleton representation. The last loss term is to encourage larger radii to avoid instability induced by surface noise.
For the sampling loss \(L_s\), we sample points on the surface of each skeleton sphere and measure the Chamfer Distance (CD) between the set of sampled sphere points \({\textbf {T}}\) and the set of sampled surface points from PointNet++ \({\textbf {P}}'\) as in Eq. 3:
Pointtosphere loss \(L_p\) measures the reconstruction error by explicitly optimizing the coordinate of skeleton points and their radii:
where \({\textbf {C}}\) is a set of predicted skeleton points, r(c) is a radius of skeleton point \({\textbf {c}}\), and \({\textbf {c}}_{p'}^{min}\) is the cloest skeleton points to a point \({\textbf {p}}'\).
Radius regularizer loss \(L_r\) is defined in Eq. 5 where \(r({\textbf {c}})\) is a radius of the skeleton point \({\textbf {c}}\). By minimizing this loss, it encourages large radii of skeleton points to make the skeleton points prediction more stable.
However, based on above three losses, predicted skeleton points can be outside of a shape if the shape is concave. Therefore, we introduce the skeletontosurface norm loss \(L_n\) to encourage the skeleton points to be inside the shape. \(L_n\) is a term to measure the reconstruction error by utilizing the property of spoke direction in MAT. Figure 5 illustrates a spoke of a skeleton point in MAT. The length of a spoke of the skeleton point \({\textbf {c}}\) is \(r({\textbf {c}})\). This is also one of our main contribution compared to [18] for skeleton points prediction. In theory, the direction of the spoke should be perpendicular to the object surface at the surface point [17]. Also the spoke direction should be pointing outside of a shape. To capture this property, we define \(L_n\):
\({\textbf {p}}_c^{'min}\) is the closest surface point to the skeleton point c and \({\textbf {n}}_{{\textbf {p}}_c^{'min}}\) is the estimated surface norm of the surface point. The “\(\cdot \)” denotes the dot product between two vectors. To estimate the norm of each point in the 3D surface points, the adjacent points are found first and then principal axis of the adjacent points using covariance analysis are calculated. More details of the norm estimation of each surface point are described in [25]. \(L_n\) encourages the skeleton points to be located within a shape even the shape is concave.
Links prediction
After predicting skeleton points, our target is to predict links to connect skeleton points to form the skeleton mesh. In theory, for any pair of skeleton points, if all points that are on the line connecting the two skeleton points are also on the skeleton surface, there should be links between those two points. We adapt the graph auto encoder (GAE) as used in [18] to predict links between skeleton points. GAE takes input the initialized adjacency matrix \(A_{ini}\) of the skeleton mesh graph \(G_{mesh}\) and the skeleton points’ features. The skeleton points’ features is concatenation of \({\textbf {C}}\), \({\textbf {R}}\), and \({\textbf {W}}^T{\textbf {F}}\). \({\textbf {C}}\) are coordinates of skeleton points, \({\textbf {R}}\) are radii of skeleton points, \({\textbf {W}}\) is the learned weights from the MLP, and \({\textbf {F}}\) is the contextual features of the surface points from PointNet++. GAE outputs the estimated adjacency matrix \(\hat{A}\) of \(G_{mesh}\). The loss function is a Masked Balanced CrossEntropy (MBCE) loss as proposed in [26].
Subcellular feature extraction from the skeleton graph
Skeleton model can be represented as the skeleton graph G(V, E) where V represents all skeleton points and E represents connection between skeleton points. The weight of the edge represents distance between skeleton points.
Neuron length computation from skeleton graph
We formulate the neuron length computation problem as finding the longest simple path in the skeleton graph G. A simple path in the graph is a path that does not have repeat nodes. In G, the length of the path is the sum of all edges’ weights along the path. Note that our skeleton graph can have loops. To solve this problem, for each node, we find the longest simple path from that node and denote as path(i) for node i. To avoid getting stuck during the loop, we mark any node when we visits as shown in the algorithm below. Next, we find \(i^*\) that maximizes path. We use the following recurrent algorithm to find the longest simple path from one node in the skeleton graph.
Neuron branch calculation
After finding the longest simple path, we are able to identify a set of nodes on that path. Those nodes are possible branching nodes. We name a set containing all possible branching nodes as \({\textbf {B}}\) For each node \(i\in B\), we find the longest simple path from that node i which does not contain any other nodes in B. Therefore, the branch is identified as the longest simple path.
Skeleton model comparison
We cluster neuron morphology by comparing different skeleton graphs. Specifically, we embed the skeleton graphs and then cluster neurons based on their embeddings. We embed the skeleton graph based on InfoGraph [19] as illustrated in Fig. 6. The embedding process is in an unsupervised manner.
First, graph convolutional layers are used to generate node features which we refer to as patch representation \({\textbf {h}}_i^j\) (i is the skeleton graph index and j is the node index of the skeleton graph i). Then graphlevel pooling is used on all patch representations to get the graph level representation (global representation) \({\textbf {H}}_i\). The mutual information (MI) estimator on globalpatch pairs over the given graph dataset \({\textbf {G}}\) is defined as:
where K is the total number of graphs in the dataset and \(G_i\in {\textbf {G}}\).
MI is the mutual information estimator modeled by the discriminator T. The JensenShannon MI estimator proposed in [19] is
where \(\mathbb {E}\) is the expectation (here it is just average operation) and \(sp(z)=log(1+e^z)\). i and \(i'\) denote two graph samples from the dataset \({\textbf {G}}\). The discriminator T estimates globalpatch representation pairs by passing two representations to different nonlinear transformations and then takes the dot product of the two transformed representations. Both nonlinear transformations consist 3 linear layers with ReLU activation functions. Therefore, the discriminator will output a score between \([0,\infty )\) to represent whether the input patch (node) is from the input graph. If the input global/patch pairs are from the same graph, we refer to them as positive samples, otherwise negative samples. We randomly sample pairs as input to the discriminator.
Dataset
Ciona neuron EM dataset
The first dataset (Dataset 1) contains two Ciona larva 3D TEM images [20]. The section thickness for TEM images varies between 35 nm and 100 nm. For each section, xy resolution is 3.85\(\times \)3.85 nm. Animal 1 contains 7671 sections and animal 2 contains about 8000 sections. In each Ciona larva, 187 neurons are annotated. Those 187 neurons can be grouped into 31 types. For animal 2, Ciona neuron skeletons are traced using TrackEM2 [27], an ImageJ [28] plugin. This dataset is summarized in Table 1, and we refer to [20] for more details.
C.elegans Neuron dataset from NeuroMorpho
NeuronMorpho [22] is a publicly available dataset that is used for neuron morphology research. It has dozens of different animals’ neurons. So far, it is the largest neuron skeletons dataset with associated metadata. In this paper, we take a subset of C.elegans dataset (Dataset 2) from the whole NeuroMorpho dataset to verify our subcellular feature extraction method and neuron skeleton comparison method. Dataset 2 consists of 1240 neuron skeletons (with radii) and it is classified into 3 different types. Each neuron with detailed metadata information such as number of branches and length of neuron.
Major type neurons from NeuroMorpho
We collect 5 major cell types from NeuromMorpho dataset (Dataset 3) including 20614 pyramidal neurons, 956 ganglion neurons, 2674 granule neurons, 1617 medium spiny neurons, and 320 motoneurons. These neurons come from 3 species, human, rat, and mouse.
In additional to the above datasets, we also use the benmark ShapeNet [29] and detailed experimental results are provided in the supplemental material.
Results
Skeleton model from 3D surface point cloud
We apply our method on animal 1 neurons from Dataset 1 for the purpose of building a shape model to analyze neuron morphology. To get the fixed number of 3D surface input points, we use the sampling strategy described in [30]. The main idea of the sampling strategy is to give each point a weight based on its distance to neighbor points. Then we sample points based on the weights until we reach the number of desired points. Details of defining the neighbor points and computing the weights are described in [30].
To evaluate our skeleton extraction method on Dataset 1, we carefully repair surface mesh using screened poisson surface reconstruction method [31] with spherical harmonics to smooth the surface. Figure 7 shows the qualitative comparison of our methods and other stateoftheart methods [18, 32]. Our method has better visual results. Reference [32] can generate unstructured skeleton points but it lacks topological constraint. We sample the number of points using [30] to be the same with the other methods for a fair comparison. It performs well when neuron has tube like structure but it is not good when neuron has a more circular shape. Compared with [18], our method can capture more detailed structures which are important for subcellular feature extraction, such as branches. For a quantitative evaluation of our method on Dataset 1, we compute the strictly defined MAT and use the handcrafted method in [18] to remove insignificant spikes to get the simplified MAT. We sample points on the simplified MAT as the ground truth skeleton points. Then we compute Chamfer Distance (CD) and Hausdorff distance (HD) between computed skeleton points and ground truth skeleton points. We refer to them as CDskel and HDskel, respectively. To compute CDskel and HDskel, we shift and rescale each skeleton so that the skeleton center is located at (0,0,0) and their x,y,z coordinates are all between −1 and 1 for every skeleton. We also measure the difference between the shapes reconstructed from the skeletons and ground truth surface points using CD and HD. We refer to them as CDrecon and HDrecon, respectively. Similarly, we also shift and rescale the ground truth points so that each neuron is centered at (0,0,0) and each neuron’s surface coordinates are between −1 and 1. Other than those four aforementioned evaluation metrics, we also use the reconstructed neuron volume difference as the evaluation metric, considering neuron volume as one of the important property of a neuron. We denote it as \(volpct\). Mathematically, it is defined as \(volpct=\frac{vreconvg}{vg}\) where vrecon is the volume of the reconstructed neurons from the skeleton model, and vg is the ground truth volume. Table 2 gives the detailed results of different methods. It shows our method has the best performance compared to other methods on Dataset 1 in terms of all of the above evaluation metrics.
Subcellular feature extraction from skeleton model
We apply our subcellular feature extraction method on Dataset 1, Dataset 2, and Dataset 3. We define the length difference percentage (lenpct) and number of branches difference percentage (numpct) to measure the neuron length error and number of branches error of the computation methods. We compare our subcellular feature extraction method with the stateoftheart subcellular feature extraction method proposed in [3]. Table 3 gives details of subcellular feature computation results. lenpct and numpct for Dataset 1 is for both animals. Our method provides the better subcellular feature extraction results in most cases and the percentage error is no more than 8 percent. Also, as we see, our method is comparable to [3] on “curve” skeletons but has better performance on the skeleton meshes (“curve” skeleton is a special case of the skeleton mesh).
For Dataset 1, we also analyze the relationships between length and branches of neurons computed from our method. For animal 1 in Dataset 1, we first use our skeleton extraction method to get the skeleton representation of each neuron. Next, we use our subcellular feature extraction method to get length of the neuron and number of branches of the neuron. Figure 8 shows the relationships between length and branches of neurons. Blue dots represent neurons from animal 1 and red dots represent neurons from animal 2. Overall, longer neurons tend to have more branches. This relationship between neuron length and number of branches is consistent between animals 1 and 2.
Skeleton model comparison
We apply our skeleton model comparison method on Dataset 1 for the purpose of analyzing Ciona neuron morphology of different neuron types. Given the skeleton graph, we embed it into a 100dimension vector. Then we use Kmeans++ to cluster vector representations of skeleton graphs. For Kmeans++, the number of clusters is set to be the same as number of neuron classes. After Kmeans++, the cluster label is given by using the majority vote of neuron types within the cluster. We use animal 1 neurons to train the Kmeans++ model and get the cluster (neuron type) centers. Next, we use animal 2 as the test set. For each neuron in the test set, we assign the label based on its closest cluster center. The distance metric we use is the euclidean distance in the vector embedding space. Table 4 shows the comparison of clustering (classification) accuracy on both training and test sets using different neuron classification method. The neuron classification methods include graph spectrum method, graph2vec [33], srep [17] and our graph level representation method. The graph spectrum method uses the eigen values of the graph’s adjacency matrix to form the vector representation. Similar to our method, graph2vec method is another way to convert the skeleton graph to the graph level vector representation. For the srep method, it uses the skeleton points’ features such as spoke direction, spoke length, and skeleton points’ locations to classify neurons. From Table 4, we observe that grouping neurons by our graph embedding provide the best classification results on both train and test sets. It shows that neuron types are closely related to its morphology. Also, our method is a better way to represent skeleton graphs in terms of clustering accuracy.
Based on previous observations, we do further morphology analysis based on our graph level representation results on Dataset 1. After we get the vector representation of each graph, we compute euclidean distance between each pair of vectors. Then we compute the inter class and intra class distance based on pairwise neuron distance as Fig. 9 shows. Diagonal entries tend to be smaller than other values, confirming a strong correlation between structure and function. More specifically, neurons within a neuron type tend to have a smaller morphological distance than neurons between different groups. Also, two animals inter and intra distance look very similar.
Based on this inter and intra class distance, we do hierarchical clustering as shown in Fig. 10. The hierarchical clustering results show that BVIN and prBTN RN have larger morphology distances from other neuron types. The BVIN neurons are a broad group of intrinsic interneurons located in the brain vesicle of Ciona. The main role of this group is to connect the sensory neurons to other groups within the brain vesicle. The BVIN neurons have partial subclassification based on sensory input [20]. Receiving specific sensory information is an indication of functional role, therefore, the BVIN can be further subdivided into different groups based on the sensory group(s) from which they receive input. Using the sensory input as criteria, the entire group was split up into four groups: prIN if receiving photoreceptor input, antIN if receiving antenna cell input, prant IN if receiving from both, and BVIN if not receiving from either. The prBTN RN only have two neurons and their functions are mostly unknown. According to the connectome [20], they receive input from both the photoreceptors and the BTN neurons (neurons involved in processing touch stimuli in the tail), so it’s possible they play a role in integrating the two inputs. Any functional differences that may exist between the two are currently unknown, however, the hierarchical clustering suggests that this may be the case.
To show the generality of our skeleton model comparison method, we apply it on Dataset 2 and Dataset 3 which both contain “curve” skeleton neurons for classification. For Dataset 2, we try to classfiy C.elegans neurons into 3 predefined cell types (somatic, motoneuron, and interneuron). For Dataset 3, we do two experiments. First, we classify each type of neurons into classes based on which species they belong to. Second, we classify neurons within same species to different neuron types. For all three experiments, we embed each skeleton mesh into a 100dimension vector and then use Kmeans++ for the classification, which is the same process for Dataset 1. We compare the performance of our method with available stateoftheart methods on NeuroMorpho dataset, LMeasure [9], TMD [11], and NeuroPath2Path [12]. We use 10fold cross validation for our classification task. Tables 5, 6 and 7 show our method has the best performance in most cases. Therefore, in general, our experiment shows that although our method is specifically designed for neuron skeleton mesh representation, it also has comparable performance in terms of “curve” skeleton neuron classification.
Conclusion
In this paper, we propose a novel neuron morphology analysis pipeline. It mainly includes three parts. First, we propose a robust shapre representation using skeleton mesh. Next, we compute subcellular features from the skeleton mesh. Finally, we compare different neuron shapes using skeleton mesh. To the best of knowledge, this is the first time that such an approach is used to represent and classify neuronal shapes.
The introduction of the estimated surface norm penalty results in a robust mesh representation that achieves the stateoftheart performance using well defined metrics.
Based on skeleton graph, we formulate subcellular feature computation problem as a longest simple path problem that can be easily computed. To compare different neuron morphology, we use a novel unsupervised graph level representation method to get the vector representation of each skeleton graphs. We provide detailed experimental results to demonstrate the effectiveness of our method. Specifically, our analysis of the Ciona dataset demonstrates that shape could be used as a significant feature for classifying neuronal types.
Availability of data and materials
The shape analysis code is available https://github.com/matlabbaltam/NeuronMorphology. Dataset 1 is available https://drive.google.com/drive/u/1/folders/1jCGek7U9e7dddsmDz0JKCkWfU0esaiJR and Dataset 2 and 3 are available https://neuromorpho.org/.
References
Halavi M, Hamilton KA, Parekh R, Ascoli GA. Digital reconstructions of neuronal morphology: three decades of research trends. Front Neurosci. 2012;6:49.
Jiang S, Pan Z, Feng Z, Guan Y, Ren M, Ding Z, Chen S, Gong H, Luo Q, Li A. Skeleton optimization of neuronal morphology based on threedimensional shape restrictions. BMC Bioinf. 2020;21(1):1–16.
Costa M, Manton JD, Ostrovsky AD, Prohaska S, Jefferis GS. Nblast: rapid, sensitive comparison of neuronal structure and construction of neuron family databases. Neuron. 2016;91(2):293–311.
Abdellah M, Hernando J, Eilemann S, Lapere S, Antille N, Markram H, Schürmann F. Neuromorphovis: a collaborative framework for analysis and visualization of neuronal morphology skeletons reconstructed from microscopy stacks. Bioinformatics. 2018;34(13):574–82.
Li S, Quan T, Xu C, Huang Q, Kang H, Chen Y, Li A, Fu L, Luo Q, Gong H, et al. Optimization of traced neuron skeleton using lassobased model. Front Neuroanat. 2019;13:18.
Saha PK, Xu Y, Duan H, Heiner A, Liang G. Volumetric topological analysis: a novel approach for trabecular bone classification on the continuum between plates and rods. IEEE Trans Med Imag. 2010;29(11):1821–38.
Panichev O, Voloshyna A. Unet based convolutional neural network for skeleton extraction. In: Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition Workshops; 2019
Lee TC, Kashyap RL, Chu CN. Building skeleton models via 3d medial surface axis thinning algorithms. CVGIP: Gr Models Image Process. 1994;56(6):462–78.
Scorcioni R, Polavaram S, Ascoli GA. Lmeasure: a webaccessible tool for the analysis, comparison and search of digital reconstructions of neuronal morphologies. Nat Protoc. 2008;3(5):866–76.
Wan Y, Long F, Qu L, Xiao H, Hawrylycz M, Myers EW, Peng H. Blastneuron for automated comparison, retrieval and clustering of 3d neuron morphologies. Neuroinformatics. 2015;13:487–99.
Kanari L, Dłotko P, Scolamiero M, Levi R, Shillcock J, Hess K, Markram H. A topological representation of branching neuronal morphologies. Neuroinformatics. 2018;16:3–13.
Batabyal T, Condron B, Acton ST. Neuropath2path: classification and elastic morphing between neuronal arbors using pathwise similarity. Neuroinformatics. 2020;18:479–508.
Li Y, Wang D, Ascoli GA, Mitra P, Wang Y. Metrics for comparing neuronal tree shapes based on persistent homology. PLoS ONE. 2017;12(8):0182184.
Schubert PJ, Dorkenwald S, Januszewski M, Jain V, Kornfeld J. Learning cellular morphology with neural networks. Nat Commun. 2019;10(1):2736.
Li Z, Fan X, Shang Z, Zhang L, Zhen H, Fang C. Towards computational analytics of 3d neuron images using deep adversarial learning. Neurocomputing. 2021;438:323–33.
Zhao J, Chen X, Xiong Z, Zha ZJ, Wu F. Graph representation learning for largescale neuronal morphological analysis. IEEE Trans Neural Netw Learn Syst 2022.
Pizer SM, Hong J, Vicory J, Liu Z, Marron J, Choi Hy, Damon J, Jung S, Paniagua B, Schulz J, et al. Object shape representation via skeletal models (sreps) and statistical analysis. Riemannian Geomet Stat Med Image Anal. 2020;1:233–71.
Lin C, Li C, Liu Y, Chen N, Choi YK, Wang W. Point2skeleton: learning skeletal representations from point clouds. In: Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition 2021;4277–4286
Sun FY, Hoffman J, Verma V, Tang J. Infograph: Unsupervised and semisupervised graphlevel representation learning via mutual information maximization. In: International Conference on Learning Representations 2019.
Ryan K, Lu Z, Meinertzhagen IA. The cns connectome of a tadpole larva of ciona intestinalis (l.) highlights sidedness in the brain of a chordate sibling. Elife. 2016;5:16962.
Scheffer LK, Xu CS, Januszewski M, Lu Z, Takemura Sy, Hayworth KJ, Huang GB, Shinomiya K, MaitlinShepard J, Berg S, et al. A connectome and analysis of the adult drosophila central brain. Elife. 2020;9:57443.
Ascoli GA, Donohue DE, Halavi M. Neuromorpho.org: a central resource for neuronal morphologies. J Neurosci. 2007;27(35):9247–51.
Qi CR, Yi L, Su H, Guibas LJ. Pointnet++: deep hierarchical feature learning on point sets in a metric space. Adv Neural Inf Process Syst. 2017;30:1.
Qi CR, Su H, Mo K, Guibas LJ. Pointnet: Deep learning on point sets for 3d classification and segmentation. Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition. 2017;652–60.
Zhou QY, Park J, Koltun V. Open3d: a modern library for 3d data processing. 2018. arXiv preprint arXiv:1801.09847.
Tran PV. Learning to make predictions on graphs with autoencoders. In: 2018 IEEE 5th International Conference on Data Science and Advanced Analytics (DSAA), 2018;237–245. IEEE
Cardona A, Saalfeld S, Schindelin J, ArgandaCarreras I, Preibisch S, Longair M, Tomancak P, Hartenstein V, Douglas RJ. Trakem2 software for neural circuit reconstruction. PLoS ONE. 2012;7(6):38011.
Collins TJ. Imagej for microscopy. Biotechniques. 2007;43(S1):25–30.
Chang AX, Funkhouser T, Guibas L, Hanrahan P, Huang Q, Li Z, Savarese S, Savva M, Song S, Su H, et al. Shapenet: An informationrich 3d model repository. 2015. arXiv preprint arXiv:1512.03012.
Yuksel C. Sample elimination for generating poisson disk sample sets. Comput Gr Forum. 2015;34:25–32.
Kazhdan M, Hoppe H. Screened poisson surface reconstruction. ACM Trans Gr (ToG). 2013;32(3):1–13.
Wu S, Huang H, Gong M, Zwicker M, CohenOr D. Deep points consolidation. ACM Trans Gr (ToG). 2015;34(6):1–13.
Narayanan A, Chandramohan M, Venkatesan R, Chen L, Liu Y, Jaiswal S. graph2vec: Learning distributed representations of graphs. 2017. arXiv preprint arXiv:1707.05005.
Acknowledgements
The authors would like to thank Kerrianne Ryan for providing the Ciona EM dataset as well as the great annotations on the dataset.
Funding
This work was supported in part by grants from NSFSSI award # 1664172, NIH # 5R01NS10377404, and Department of Energy Office of Science grant # DESC002197.
Author information
Authors and Affiliations
Contributions
JJ designed the overall neuron morphology analysis pipeline, carried out experiments on different datasets, and drafted manuscript. MG helped on the clustering/classification part of the experiments and helped draft the manuscript. CB provided the Ciona neuron class information and helped draft the paper. WS coordinated Ciona dataset acquisation, helped analyzed the Ciona data, and reviewed the manuscript. BSM coordinated the overall design, development and evaluations of the neuron morphology analysis method, and helped prepare the manuscript.
Corresponding authors
Ethics declarations
Ethics approval and consent to participate
Ethical approval was not required as confirmed by the license attached with the openaccess data.
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.
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.
About this article
Cite this article
Jiang, J., Goebel, M., Borba, C. et al. A robust approach to 3D neuron shape representation for quantification and classification. BMC Bioinformatics 24, 366 (2023). https://doi.org/10.1186/s1285902305482y
Received:
Accepted:
Published:
DOI: https://doi.org/10.1186/s1285902305482y