 Methodology article
 Open Access
 Published:
Skeleton optimization of neuronal morphology based on threedimensional shape restrictions
BMC Bioinformatics volume 21, Article number: 395 (2020)
Abstract
Background
Neurons are the basic structural unit of the brain, and their morphology is a key determinant of their classification. The morphology of a neuronal circuit is a fundamental component in neuron modeling. Recently, singleneuron morphologies of the whole brain have been used in many studies. The correctness and completeness of semimanually traced neuronal morphology are credible. However, there are some inaccuracies in semimanual tracing results. The distance between consecutive nodes marked by humans is very long, spanning multiple voxels. On the other hand, the nodes are marked around the centerline of the neuronal fiber, not on the centerline. Although these inaccuracies do not seriously affect the projection patterns that these studies focus on, they reduce the accuracy of the traced neuronal skeletons. These small inaccuracies will introduce deviations into subsequent studies that are based on neuronal morphology files.
Results
We propose a neuronal digital skeleton optimization method to evaluate and make fine adjustments to a digital skeleton after neuron tracing. Provided that the neuronal fiber shape is smooth and continuous, we describe its physical properties according to two shape restrictions. One restriction is designed based on the grayscale image, and the other is designed based on geometry. These two restrictions are designed to finely adjust the digital skeleton points to the neuronal fiber centerline. With this method, we design the threedimensional shape restriction workflow of neuronal skeleton adjustment computation. The performance of the proposed method has been quantitatively evaluated using synthetic and real neuronal image data. The results show that our method can reduce the difference between the traced neuronal skeleton and the centerline of the neuronal fiber. Furthermore, morphology metrics such as the neuronal fiber length and radius become more precise.
Conclusions
This method can improve the accuracy of a neuronal digital skeleton based on traced results. The greater the accuracy of the digital skeletons that are acquired, the more precise the neuronal morphologies that are analyzed will be.
Background
Neurons are the basic structural unit of the brain. Their morphology is a key determinant of neuron type classification [1,2,3,4,5]. The morphology of a neuronal circuit is a fundamental component for neuron modeling [6]. With improvements in sparse labeling and optical microscopy techniques, scientists are able to acquire threedimensional wholebrain images at a submicronresolution for mammals [7,8,9,10,11]. From these images, the techniques that accurately transform them into digital descriptions of neuron morphology are very important; this transformation process is called neuron tracing. Neuron tracing is an important technique in neuroscience that includes many automatic tracing methods, such as skeletonization [12, 13], minimum spanning trees [14], snake models [15, 16], principle curve models [17] and neural network methods [18]. These automatic tracing methods are always used in small image blocks and work well. However, when dealing with a wholebrain image dataset, researchers prefer semimanual tracing [2, 8, 19,20,21] (researchers manually mark the neuron fiber nodes according the images, and software automatically links these nodes to create a neuron morphology file).
Existing neuron tracing studies mostly focus on the correctness and completeness of neuronal morphology. For neuron tracing from wholebrain images, the results of semimanual tracing are credible, so they are chosen as the ground truth in many studies to evaluate automatic method performance [17, 22]. However, in a semimanually traced neuronal skeleton, we note some inaccuracies. The distance between consecutive nodes is very long, spanning multiple voxels. The nodes are marked around the centerline of the neuronal fiber, not on the centerline. Although these inaccuracies will not seriously affect the topology of the neuronal morphology file, they will reduce the accuracy of the traced neuron skeleton. Since the ability of people to distinguish local voxels with little difference in grayscale is not as accurate and objective as that of a computer, inaccuracies generally appear in semimanually traced neuronal skeletons and may also appear in some automatically traced neuronal skeletons. These small inaccuracies will introduce deviations into subsequent tasks based on neuronal morphology files, such as neuronal fiber radius estimation and neuronal fiber length calculation. At present, with the emergence of more longprojection neuronal skeletons, it is urgent for us to improve the accuracy of these neuronal skeletons.
Here, we propose a neuronal digital skeleton optimization method to adjust a traced neuronal skeleton to the centerline of the neuronal fiber in images. This method evaluates and finely adjusts the skeleton according to threedimensional shape restrictions, which are designed on the basis of the smooth and continuous properties of neuronal fibers. This method reduces the distance between the traced neuronal skeleton and the centerline of the neuronal fiber in images, thereby increasing the accuracy of the traced neuronal skeleton. Furthermore, based on the optimized skeletons, we calculate the neuronal morphology parameter to evaluate its effect. The results show that the more accurate the skeleton is, the more the accuracy of morphology metrics, such as neuronal fiber lengths and radiuses, are improved.
Methods
Dataset and computing environments
We used two datasets for testing. The first is provided by Wang et al. [11], and it consists of genetic labeling chandelier cells imaged by the fMOST system with a resolution of 0.2 × 0.2 × 1 μm^{3}. The second is provided by Zhang et al. [23] and Jiang et al. [24], and it consists of pyramidal neurons at the same resolution. Two mice were used, one for the chandelier cells and the other for the pyramidal neurons. The mice were anesthetized deeply using a 1% solution of sodium pentobarbital and sacrificed by injection of ketamine. The semiautomatic neuron tracing used in the quantitative evaluation was performed with Amira software (version 6.1.1; FEI, Mérignac Cedex, France); the code of this method was written in MATLAB (R2019a, MathWorks) and run on a graphics workstation (Xeon CPU E5–2687 v3).
Workflow of the proposed method
Figure 1 depicts the computational workflow, including its three stages. The skeleton optimization method is the second stage, and the first and third stages are the pre and postprocessing of the neuronal morphology, respectively.
The first phase is to extract a branch from the digital skeleton of the tree structure (Steps 1–3 in Fig. 1). Specifically, we linearly interpolate the digital skeleton and select the branch to be adjusted. The second stage involves finetuning the points derived from the digital skeleton to the centerline of the neuronal fibers from the images (Steps 4–5 in Fig. 1). Specifically, we initially adjust the digital skeleton points that are far from the centerline of the neuronal fibers based on the image grayscale. Subsequently, provided that the neuron fiber shape is smooth and continuous, we describe this physical property with two shape restrictions. One restriction is designed based on the image grayscale: the difference in the signal brightness between the local voxels around the neuronal fiber centerline must be small. The other is designed based on geometry: the angle formed from three continuous points on the neuronal fiber centerline must not be extreme. These two restrictions are designed to finely adjust the digital skeleton points to the neuronal fiber centerline. Subsequently, we join the branches into a tree structure according to the connected relationships, which comprises the third stage (Step 6 in Fig. 1). Finally, the digital skeleton yields a digital neuron morphological file.
With regard to the above three stages, the second stage is a looping structure. See the flowchart in Fig. 1 for details.
Skeleton point adjustment based on grayscale
According to the grayscale distribution of the voxels around the neuronal fibers, a preliminary adjustment of the digital skeleton points that are far away from the neuronal fiber centerline shown in the image is performed. We take advantage of the idea of using the mean shift to find the location of this centerline point. Because of anisotropy in the image stacks, we linearly interpolate the threedimensional images in the Z direction to obtain images with the same resolution in the X, Y and Z directions. For a skeleton point p_{0}(x_{0}, y_{0}, z_{0}), take an adjustment in the Z direction as an example. We keep the X and Y coordinates unchanged and calculate the convolution of the voxels in the local range centered on p_{0} and the Gaussian kernel function in the Z direction. The point that has the largest convolution value is regarded as the centerline point p_{1} of the neuronal fibers from the images.
Similarly, the X and Y coordinates of the skeleton point p_{1} are adjusted to the centerline point to finally obtain the adjusted point p_{2}.
Skeleton point adjustment based on shape restrictions
To describe the tortuous neuronal fiber morphology with more accurate digital skeletons, it is necessary to finetune the skeleton points after the preliminary adjustment. The neuronal fibers are smooth and continuous, and this neuronal fiber property is shown in the images. Considering this property, for the point p_{2}, we designed an angle restriction and brightness restriction to judge whether p_{2} needs further fine adjustment. Specifically, for p_{2}, we determined its parent point p_{f}, the parent point of p_{f}, p_{ff}, and the child point p_{c}. The point p_{2} forms an angle ∠p_{f}p_{2}p_{c} with p_{f} and p_{c}. Because the neuron fibers are smooth, the angle ∠p_{f}p_{2}p_{c} tends to be close to 180°. We set the critical value of the angle α_{std} as an obtuse angle (135° is chosen in this dataset). If ∠p_{f}p_{2}p_{c} < α_{std}, p_{2} is considered to require further adjustment. For the brightness restriction, we compared the signal strength \( {G}_{p_2} \) at point p_{2} with the signal strength threshold G_{std}, thus designing G_{std} as follows:
Here, \( {G}_{p_f} \) is the signal strength of p_{f}, taking the grayscale value of p_{f} in the same way as for \( {G}_{p_2} \) and \( {G}_{p_c} \). G_{Branch _ mean} and G_{Branch _ stdev} are the mean and standard deviation of the signal strength values of all points on the current branch. If \( {G}_{p_2}>{G}_{std} \), the signal strength of p_{2} is in accordance with that of the current branch. However, if \( {G}_{p_2}<{G}_{std} \), p_{2} is identified as a point that requires further adjustment. If its signal strength is small, this is contrary to the fact that the signal strength at the center of the neuronal fiber in images tends to be relatively strong.
According to these restrictions, we set out to find the centerline points of the neuronal fiber from the images. Prior to adjusting p_{2}, because p_{2} itself may be a biased point, we have to judge whether we need to find a new starting point other than p_{2} for subsequent computation. Here, if ∠p_{f}p_{2}p_{c} > α_{std} and \( {G}_{p_2}<{G}_{std} \), then there is only a problem with the brightness of point p_{2}. In this situation, the starting point p_{2 _ new} will be equal to p_{2}. If ∠p_{f}p_{2}p_{c} < α_{std}, this indicates that there is a serious deviation in the angle pertaining to p_{2}, and such a point is likely to have a problem with both brightness and angle. Therefore, we have to find a new starting point p_{2 _ new}. The assumption here is that \( \left\Vert \overrightarrow{p_2{p}_f}\right\Vert <\left\Vert \overrightarrow{p_2{p}_c}\right\Vert \). In the \( \overrightarrow{p_2{p}_c} \) direction, with p_{2} as the starting point, we draw \( \left\Vert \overrightarrow{p_2{p}_{rc}}\right\Vert =\left\Vert \overrightarrow{p_2{p}_f}\right\Vert \) and obtain a point p_{rc} on \( \overrightarrow{p_2{p}_c} \). The intersection of the vertical bisector of \( \overrightarrow{p_f{p}_{rc}} \) and \( \overrightarrow{p_f{p}_{rc}} \) is defined as the point p_{mid}. In the \( \overrightarrow{p_f{p}_{ff}} \) direction, we draw \( \left\Vert \overrightarrow{p_f{p}_s}\right\Vert =\left\Vert \overrightarrow{p_2{p}_f}\right\Vert \) and obtain a point p_{s} on the extension line of \( \overrightarrow{p_f{p}_{ff}} \) (the topological graph is shown in Figure S1 A of the supplementary materials). Then, p_{2 _ new} is located as:
Here, \( {G}_{p_{mid}} \) represents the signal strength of p_{mid}, taking the grayscale value of p_{mid} in the same way as for \( {G}_{p_2} \) and \( {G}_{p_s} \). With the starting point p_{2 _ new}, we use the idea of the mean drift to search for the centerline points of the neuronal fiber from the images. The drift direction \( {\overset{\rightharpoonup }{\upsilon}}_2 \) in a spherical region S_{p} centered on p_{2 _ new} is calculated as follows:
where p_{i} is the ith voxel in the spherical region S_{p}. \( {a}_{p_i} \) is the angle factor of p_{i}, \( {a}_{{\mathrm{p}}_i}=\sin {\theta_{p_i}}^2 \), where \( {\theta}_{p_i} \) is the angle between the vectors \( \overrightarrow{p_f{p}_{rc}} \) and \( \overrightarrow{p_{2\_ new}{p}_i} \), and \( {b}_{p_i} \) is the signal strength factor of p_{i}, \( {b}_{p_i}=\exp \left({\left({G}_{p_i}\mu \right)}^2/2{\sigma}^2\right)\exp \left({\left({G}_{p_{2\_ new}}\mu \right)}^2/2{\sigma}^2\right) \), where , and \( \sigma =\left\mu {G}_{p_{2\_ new}}\right \).
In the \( {\overset{\rightharpoonup }{\upsilon}}_2 \) direction, we search for the point
\( {G}_{p_i} \) is the signal strength at voxel p_{i}, taking the grayscale value of p_{i}. \( {G}_{p_{2\_ new}} \) is the signal strength of p_{2 _ new}. This means that ω_{1} and ω_{2} are the parameters that can be adjusted according to the image quality. When the parameter ω_{1} is set to a larger value, this means that we regard the grayscale value as a large influence on the points. However, when ω_{2} is set to a larger value, this means that we regard the angle deviation of the points as having a greater influence.
Quantitative evaluation on the synthetic dataset
Since the goals of most automatic neuron tracing methods are different from ours, the semimanually traced neuron skeleton, which is regarded as the ground truth in most studies, is one of our targets that needs to be optimized. Therefore, we have to find a ground truth, rather than the semimanually traced results, to evaluate the performance of our methods. Here, we design a threedimensional synthetic neuron fiber image to evaluate our methods (Fig. 2). First, we choose a helix curve (t∈(π,3π], x = 30 × sin(t) + 60, y = 30 × cos(t) + 60, z = 10 × t) as the centerline of the neuron fiber, which is the ground truth of the traced neuron skeleton. Then, in threedimensional images, we use continuous voxels to represent this helix curve. Next, we choose the Gaussian function as the point spread function (PSF) of the imaging system and convolve it with the helix curve images. The result simulates the neuron fiber images. Based on these images, we semimanually trace the neuron fibers. Following a previous procedure [8, 11, 19], the semimanually traced neuron skeleton is completed interactively in Amira visualization and data analysis software. The skeleton is checked backtoback by 3 people. The traced skeleton is regarded as the optimization target. Finally, we use our method to optimize the skeleton.
The noise of nearby neuron fibers and other fibers can affect the quality of the automatically traced neuron skeleton. To evaluate the performance of our method in such cases, we add Gaussian noise with different signaltonoise ratio levels (Fig. 2c and d). For the case of nearby neurons, we set a straight line as representing another neuron whose distance from the helix curve changes (Fig. 2e and f).
In the research on automatic neuron tracing methods, there are various parameters used to evaluate the traced results from different perspectives [22, 25, 26]. Among these parameters, the significant spatial distance (SSD) is a parameter that can intuitively evaluate the distance between two neuron skeletons, and it is suitable for us to use to quantitatively evaluate the difference between the traced neuron skeleton and the ground truth. The SSD is defined as the average distance of the reciprocal minimal spatial distances between two neuron skeleton nodes. A larger SSD means that these two neuron skeletons are more different from each other. In Peng’s work [22], if the distance between two nodes is less than two pixels, it is difficult to determine the difference visually. As a result, when calculating the SSD, the minimum reciprocal distance is set to 2 pixels. However, in our experiment, the maximum reciprocal distance between the unoptimized skeleton and the ground truth is less than 2 voxels. Therefore, we set the minimum reciprocal distance to 1 voxel. In addition, we also provide the SSDs that the minimum reciprocal distances are set to 0 and 0.5 voxels.
Process for a wholebrain threedimensional image dataset
We used image blocks from TDat [27] and converted the threedimensional image dataset, which is 32,266 × 54,600 × 10,730 voxels (more than 30 TB) with a depth of 16 bit/voxel, to the TDat data format. TDat divides a threedimensional image dataset into many small threedimensional image blocks. We then split the reconstructed digital skeleton into many fragments according to the range of each TDat image block, and the skeleton fragments in each block were calculated in turn. Finally, the fragments were connected to become a complete digital skeleton. This process was slightly different from the workflow for the local neuronal fibers. In summary, for the locally distributed neuronal fibers, we split the digital skeletons according to the branch structure. However, for the longprojection neuronal fibers, we split the digital skeletons according to the range of the blocks.
Neuronal fiber radius estimation
In a real dataset, it is difficult to obtain the centerline of the neuron fiber images. Thus, we do not have a way to directly evaluate the difference between the traced neuron skeletons and the ground truth. However, according to the traced neuron skeleton, we can calculate the neuron morphology parameters to evaluate the differences in the sizes. Here, we choose the neuron fiber length and neuron fiber radius to evaluate the performance of our method.
The estimation process for the neuron fiber radius is as follows (see Fig. 5): Based on the topology of the digital skeleton, we calculate a crosssectional plane of a skeleton point. Thereafter, with a series of image processing operations, this plane is used to segment the area of the neuronal fiber tubular structure from the crosssectional image, and we estimate the neuronal fiber radius. Specifically, the skeleton point p_{3} and its parent point p_{f} and child point p_{c} can be approximately considered to be in the same circle with a very large radius. Moreover, the normal vector \( {\overset{\rightharpoonup }{\upsilon}}_3 \) of this circle at p_{3} can be regarded as parallel to \( \overrightarrow{p_f{p}_{rc}} \) (the topological graph is shown in the supplementary materials, Figure S1 C). Subsequently, the crosssectional plane at p_{3} is defined by the normal vector \( {\overset{\rightharpoonup }{\upsilon}}_3 \), named Area_{3}. With regard to Area_{3}, we segment the neuronal fiber tubular structure with a series of image processing operations, including linear grayscale transformations, adaptive thresholding [28], open operations, and connected component operations on the crosssectional image. Finally, we obtain the inner part of the crosssection of the neuronal fiber structure at p_{3}, named Area_{t}. Here, we estimate this segmented result with a circle and estimate the radius according to the formula of a circle \( r=\sqrt{S/\pi } \), where S = Area_{t}.
Statistical analysis
The statistical significance analysis was performed using MATLAB, and statistical comparisons were performed using a twotailed ttest. All measurements are listed in the form mean ± std. The confidence level was set to 0.05 (*p < 0.05, **p < 0.01, ***p < 0.0001).
Results
The performance of the method on synthetic neuron fiber images
To verify that our method works as intended, we apply our method to synthetic images in different situations (Fig. 2). We verified three cases, including an ideal neuron fiber without any interference, noisy images, and a case where other neuron fibers are near the target neuron.
According to the results of the first case (Fig. 2b), the SSD of the optimized neuron skeleton is smaller than that of the original neuron skeleton. After making adjustments with the fiber property, the minimum reciprocal distance is not larger than 1 voxel. When the minimum reciprocal distance is set to 0 and 0.5 voxels, the SSD of the skeleton adjusted by the fiber property is slightly smaller than the SSD of the neuron skeleton adjusted by the grayscale value. This means that after optimization, the difference between the traced skeleton and the ground truth is reduced.
In wholebrain threedimensional images, there is considerable noise. To simulate this case, we test our method under different signaltonoise ratios (SNRs, Fig. 2c and d). The results show that the SSD increases with decreasing SNR. When the SNR is larger than 5.2 dB, the SSD of the optimized neuron skeleton is smaller than that of the original skeleton, and this trend changes slowly. However, when the SNR is lower than 5.2 dB, the SSD of the optimized neuron skeleton suddenly increases. These results suggest that our method can work well with SNRs larger than 5.2 dB. The images from the reported wholebrain imaging system [9, 10] are sufficient for our method to work.
In neuron tracing research, distinguishing two neuron fibers that are close to each other is a very difficult situation, and we test our method in such simulations. The results (Fig. 2e and f) show that the SSD of the optimized skeleton first increases and then decreases, and the SSD reaches a peak at a distance of 6 voxels between two neuron fibers. From the crosssections of the neuron fibers in Fig. 2e, if the distance between the two neuron fibers is less than 6 voxels, it is difficult to visually distinguish the boundary of the two neuron fibers. In this case, our method performs worse than semimanual tracing. If the distance between two neuron fibers is greater than 6 voxels, the boundary can be roughly distinguished by humans. In this case, the SSD of the optimized skeleton is similar to that of the semimanually traced skeleton. This result suggests that our method does not perform very well in distinguishing adjacent neurons. If humans cannot distinguish the boundaries of nearby neuron fibers, our method will be severely affected by adjacent neuron fibers. If humans can distinguish the boundaries of adjacent neuron fibers, the accuracy of our methods is not much different from the accuracy of semimanual tracing. Although our method is not better than semimanual tracing when dealing with adjacent neuron fibers, the skeleton processed by the machine is more objective and repeatable.
In real applications, adjacent neuronal fibers often have different signal intensities compared with the target neuronal fiber. We set different grayscales for adjacent neuronal fibers in the synthetic data and tested our method on these data as shown in Figure S2 in the additional file. The results show that when the difference in signal strength between two adjacent neuronal fibers is great, the SSD score of our method is low, which indicates that our method performs well. In Fig. 2f, the two neuronal fibers have the same signal intensity, which is an extreme situation for our method.
Applicability to the fibers of pyramidal neurons
For pyramidal neurons with longprojection axonal fibers, researchers prefer semimanual tracing. Due to the great fiber length of this type of neuron, slight inaccuracies accumulate and increase, which may cause large deviations.
In applying this method to dendritic and axonal fibers, the image dataset comes from the fMOST imaging system. Figure 3 shows the axonal and dendritic fiber adjustment results for a pyramidal neuron in the mouse motor cortex. From Fig. 3A and B, the original skeletons describe the general direction of the neuron fibers (Fig. 3, yellow line). The enlarged images (Fig. 3A1A8 and B1B8) show that it is necessary to improve the accuracy of the traced skeletons to describe the tortuous shape of these neuron fibers. After the linear interpolation and adjustment based on grayscale, the number of skeleton points increases, and the tortuous morphology becomes more accurate (Fig. 3, magenta line). However, after grayscale adjustment, a jagged line appears in the skeleton (magenta line, Fig. 3, B4, B7 and B8). After the adjustment, based on the shape restrictions, the skeletons become smoother and match the centerline of the neuron fibers in the images more closely (Fig. 3, green line).
For longprojection neuron fibers in wholebrain images, Fig. 5A shows the results. In Fig. 5A2A4, the optimized skeleton matches the centerline of the fibers from the images more closely.
Applicability to the fibers of chandelier cells
Chandelier cells are GABAergic interneurons whose axonal fibers are slender and tortuous. Because their axon arbors are extremely dense and complex, researchers prefer semimanual tracing to ensure the correctness of their topology. We choose fiber images of these neurons as a representative case to test our method. Figure 4 shows the axonal and dendritic fiber optimized results for a chandelier cell. Compared with the pyramidal neuron fibers in Fig. 3, the axonal fibers of the chandelier cell are thinner and more tortuous. Although in Fig. 4A and B, the axonal fibers seem very dense, after we remove the voxels near the axonal fibers from the threedimensional image (Fig. 4B1), the axonal fibers look less difficult for machines to identify.
In Fig. 4B1 and C1, the traced skeletons describe the general direction of the fibers (Fig. 4, yellow line). However, the enlarged images (Fig. 4B2B4 and C2C4) show that it is necessary to improve the accuracy of the digital skeleton to describe the tortuous shape of these dendritic and axonal fibers. After the optimization, the number of skeleton points increases, and the tortuous morphology becomes more accurate (Fig. 4, green line). Compared with the original traced skeleton, the optimized skeleton matches the centerline of the dendritic and axonal fibers from the images more closely (Fig. 4, B2B4 and C2C4).
Influence of the optimized skeleton on the neuron morphological parameters
Changes in the digital neuron skeleton will inevitably have an impact on the calculation results of morphological parameters. Here, we select the neuron fiber radius, which is always calculated based on the digital neuron skeleton, and the neuron fiber path length to evaluate the performance of our method.
Figure 5CE show the information of the pyramidal neuron fiber radius. According to the geometric relationship of the continuous nodes in the neuron skeleton, we extract the cross sections of the dendritic and axonal fibers (Fig. 5C; additional crosssections are shown in the supplementary materials, Figure S3). These crosssections show that the dendritic fiber radius is larger than the axonal fiber radius. The shape of the neuron fiber crosssections is not a circle but an irregular ellipse. We select several fibers from different parts of the neuron, colorcoded in Fig. 5D, and calculate their average fiber radius. The resulting statistic in Fig. 5D2 shows that the dendritic radius is significantly larger than the axonal radius. The radius of the major axonal fibers is significantly different from the radius of the terminal axonal branch fibers. Furthermore, we show the neuron fiber radius distribution with a heat map in Fig. 5E. These heat maps (Fig. 5E1 and E2) show that the radius of the major axonal fiber seems larger than those of the fibers branching from the major axonal fiber, and the size of the dendritic fiber radius does not show an obvious relationship with the branch order. These results are consistent with previous studies [29, 30]
The path lengths of the neuron fibers are listed in Table 1. Overall, the optimized skeletons are longer than the original skeletons. This may be caused by the optimized skeleton depicting the tortuous shape of neuron fiber in the images more accurately. In the simulation experiment, the original traced skeleton is 5.67% shorter than the ground truth. The optimized skeleton is 3.80% longer than the ground truth. This suggests that the path length of the optimized skeleton is closer to the ground truth.
Discussion
The goal of our method is to reduce the difference between the traced neuron skeleton and the centerline of the neuron fiber in images rather than correctly reconstructing the topology of the neuron skeleton. In most studies, the topology of the neuron fiber skeleton carefully traced by researchers is delicate enough to enable the projection patterns to be analyzed. Therefore, the optimized skeleton improves the accuracy of the morphological parameters. In this study, the original skeletons are semimanually traced because inaccuracy generally appears in the semimanually traced results of the longprojection neuron morphology. However, if automatically traced neuron skeletons have the same problem, it is possible to use our method to optimize the traced neuron skeleton or integrate our method into the process of neuron reconstruction to obtain a more accurate neuron fiber skeleton.
In fact, neuron fibers are not as ideal as we have seen in synthetic data. From the images of real neuron fibers in Fig. 5, the crosssections of fibers are not circular (Fig. 5C), and the spine and bouton grow on the fibers (Fig. 5A2 and A3). These factors influence the recognition of the neuron fiber centerlines. For example, if the spine is considered to be part of a neuron fiber, the centerline may be biased toward the side where the spine exists. If not, we may have to use algorithms [31, 32] to identify these tiny structures first to avoid their influence. In this article, we optimize the skeleton without distinguishing these tiny structures. However, in the future, distinguishing these tiny structures may further improve the accuracy of the skeleton.
In our method, we use the mean shift to find the nodes on the centerline. Without the limitation of the search range, the results of the mean shift will be easily affected by noise and nearby neuron fibers. To avoid these effects, we set the search range according to the full width at half maximum of the gray distribution of a truncated line in the maximumprojection images of the neuron fibers. Therefore, when our method is used for different neuron fibers, the range must be set according to the specific situation. For example, the ranges of the axonal fibers of chandelier cells are smaller than the dendritic fibers of pyramidal neurons. On the other hand, the radius of neuron fibers is very important information for our method, especially in areas with dense fibers. After skeleton optimization, we compare the optimized skeleton with the original skeleton. If the optimized skeleton is very strange, we estimate the radius of the neuron fibers based on the original skeleton. According to the radius information, the program can optimize the traced skeleton again.
In the part of skeleton adjustment that uses shape restrictions, the reason we choose the property of neuron fibers to finetune the neuron skeleton is that neuron fiber shapes are variable, but this property is stable locally. In the finetuning step, there are two main thresholds: the angle threshold and grayscale threshold. The values of these two thresholds are adjustable. In the method section, the angle threshold is set to 135°. However, for tortuous neuron fibers such as the axonal fibers of chandelier cells, this threshold has to be set to a smaller value. For the grayscale threshold, if the local signal intensity varies greatly, the ranges of the gray degree (formula 2) have to be set to larger values. We use the angle threshold and grayscale threshold to limit the adjusting range of the skeleton nodes in a more detailed way. According to different images of the neuron fibers, these thresholds may have to be adjusted. In the methods section, we give the reference values of these thresholds, which does not mean that the values are suitable for all images.
In the radius estimation, we use only a crosssection to calculate the fiber radius. Although we use the adaptive method [28] to set the threshold for image binarization in different crosssections, the impacts of noise and signal intensity variation cannot be avoided completely. In addition, in previous studies, images of the same neuron type labeled with different methods showed different neuron fiber radii [30]. Even if the imaging system and labeling method are the same, the tiny structure will affect the radius estimation results. For example, if spines are counted as part of the dendritic fibers, the estimated radius will be larger than the radius not including the spines. In the future, it may be more appropriate to use a threedimensional image instead of a crosssection to estimate the radius, and more information may contribute to the robustness to noise or the signal intensity variation. In contrast, in this study, we used a circle to approximate the neuronal fiber tubular structure from crosssectional images. However, from the crosssectional images shown in the results, the shape of the tubular fiber is more like an irregular ellipse shape. In the future, a more complex model must be considered in order to describe neuronal fiber tubular structures.
On the other hand, in our experiments, the adjustmentbased fiber property step helps to reduce the jagged lines introduced by adjustments based on grayscale. In our method, we take voxels as the minimum calculation unit; as a result, the optimized coordinates are integers. These integers may also introduce jagged lines. The jagged lines can be smoothed by performing Gaussian smoothing on the coordinates of the ordered skeleton nodes. When the nodes are evenly distributed around the ground truth, the performance of Gaussian smoothing is very good. However, when the nodes are distributed on one side of the line, Gaussian smoothing will pull the line away from the ground truth. At present, we have not completely clarified the situations when Gaussian smoothing is not effective. To ensure the stability of the algorithm, we have not yet added Gaussian smoothing in skeleton optimization. In the future, we believe that after clarifying these uncertain situations, adding Gaussian smoothing to neuron skeleton optimization will further improve the accuracy of the skeletons.
Conclusion
In conclusion, this method focuses on improving the accuracy of the traced skeletons of neuronal fibers after neuron tracing instead of improving the completeness and correctness of neuron tracing. In the skeleton optimization stage, we take advantage of the neuronal fiber property shown in images and design shape restrictions and an evaluation index for the adjustment computation. With this method, the accuracy of the digital skeleton is improved. Compared with the morphology obtained without optimization, the optimized neuronal fiber is longer. Furthermore, based on the more accurate neuronal morphology, we analyze the neuronal fiber radius distribution and find that the major branches are thicker than the fibers split from major branches. In the future, we believe that with this method, more accurate digital skeleton data will improve the results of neuron modeling.
Availability of data and materials
All the image datasets supporting the conclusions of the article are available, and can be obtained via contacting the corresponding author A. Li (aali@mail.hust.edu.cn) or found at http://atlas.brainsmatics.org/a/jiang2004. The codes of the method are also available via above webpage.
Abbreviations
 VIS:

Visual cortex
 MC:

Motor cortex
 fMOST:

Fluorescence microoptical sectioning tomography
 SNR:

Signal to noise ratio
 SSD:

Significant spatial distance
 MD:

Minimal distance
 PSF:

Point spread function
References
 1.
Zeng H, Sanes JR. Neuronal celltype classification: challenges, opportunities and the path forward. Nat Rev Neurosci. 2017;18(9):530–46.
 2.
Sun Q, Li X, Ren M, Zhao M, Zhong Q, Ren Y, Luo P, Ni H, Zhang X, Zhang C, et al. A wholebrain map of longrange inputs to GABAergic interneurons in the mouse medial prefrontal cortex. Nat Neurosci. 2019;22(8):1357–70.
 3.
Economo MN, Viswanathan S, Tasic B, Bas E, Winnubst J, Menon V, Graybuck LT, Nguyen TN, Smith KA, Yao Z, et al. Distinct descending motor cortex pathways and their roles in movement. Nature. 2018;563(7729):79–84.
 4.
Petilla Interneuron Nomenclature Grou (PING). Petilla terminology: nomenclature of features of GABAergic interneurons of the cerebral cortex. Nat Rev Neurosci. 2008;9(7):557.
 5.
DeFelipe J, LopezCruz PL, BenavidesPiccione R, Bielza C, Larranaga P, Anderson S, Burkhalter A, Cauli B, Fairen A, Feldmeyer D, et al. New insights into the classification and nomenclature of cortical GABAergic interneurons. Nat Rev Neurosci. 2013;14(3):202–16.
 6.
Markram H, Muller E, Ramaswamy S, Reimann MW, Abdellah M, Sanchez CA, Ailamaki A, AlonsoNanclares L, Antille N, Arsever S, et al. Reconstruction and simulation of neocortical microcircuitry. Cell. 2015;163(2):456–92.
 7.
Economo MN, Winnubst J, Bas E, Ferreira TA, Chandrashekar J. Singleneuron axonal reconstruction: the search for a wiring diagram of the brain. J Comp Neurol. 2019;527(13):2190–9.
 8.
Lin R, Wang R, Yuan J, Feng Q, Zhou Y, Zeng S, Ren M, Jiang S, Ni H, Zhou C, et al. Celltypespecific and projectionspecific brainwide reconstruction of single neurons. Nat Methods. 2018;15(12):1033–6.
 9.
Economo MN, Clack NG, Lavis LD, Gerfen CR, Svoboda K, Myers EW, Chandrashekar J. A platform for brainwide imaging and reconstruction of individual neurons. Elife. 2016;5:e10566.
 10.
Gong H, Xu D, Yuan J, Li X, Guo C, Peng J, Li Y, Schwarz LA, Li A, Hu B, et al. Highthroughput dualcolour precision imaging for brainwide connectome with cytoarchitectonic landmarks at the cellular level. Nat Commun. 2016;7:12142.
 11.
Wang X, Tucciarone J, Jiang S, Yin F, Wang BS, Wang D, Jia Y, Jia X, Li Y, Yang T, et al. Genetic single neuron anatomy reveals fine granularity of cortical axoaxonic cells. Cell Rep. 2019;26(11):3145–59 e3145.
 12.
He W, Hamilton TA, Cohen AR, Holmes TJ, Pace C, Szarowski DH, Turner JN, Roysam B. Automated threedimensional tracing of neurons in confocal and brightfield images. Microsc Microanal. 2003;9(4):296–310.
 13.
Turetken E, Benmansour F, Fua P. Automated reconstruction of tree structures using path classifiers and mixed integer programming. Proc Cvpr Ieee. 2012:566–73.
 14.
Yang J, Hao M, Liu X, Wan Z, Zhong N, Peng H. FMST: an automatic neuron tracing method based on fast marching and minimum spanning tree. Neuroinform. 2019;17(2):185–96.
 15.
Wang Y, Narayanaswamy A, Tsai CL, Roysam B. A broadly applicable 3D neuron tracing method based on opencurve snake. Neuroinform. 2011;9(2–3):193–217.
 16.
Kass M, Witkin A, Terzopoulos D. Snakes: active contour models. Int J Comput Vis. 1988;1(4):321–31.
 17.
Quan T, Zhou H, Li J, Li S, Li A, Li Y, Lv X, Luo Q, Gong H, Zeng S. NeuroGPStree: automatic reconstruction of largescale neuronal populations with dense neurites. Nat Methods. 2016;13(1):51–4.
 18.
Zhou Z, Kuo HC, Peng H, Long F. DeepNeuron: an open deep learning toolbox for neuron tracing. Brain Inform. 2018;5(2):3.
 19.
Li X, Yu B, Sun Q, Zhang Y, Ren M, Zhang X, Li A, Yuan J, Madisen L, Luo Q, et al. Generation of a wholebrain atlas for the cholinergic system and mesoscopic projectome analysis of basal forebrain cholinergic neurons. Proc Natl Acad Sci U S A. 2018;115(2):415–20.
 20.
Han Y, Kebschull JM, Campbell RAA, Cowan D, Imhof F, Zador AM, MrsicFlogel TD. The logic of singlecell projections from visual cortex. Nature. 2018;556(7699):51–6.
 21.
Gerfen CR, Economo MN, Chandrashekar J. Long distance projections of cortical pyramidal neurons. J Neurosci Res. 2018;96(9):146775.
 22.
Peng H, Ruan Z, Atasoy D, Sternson S. Automatic reconstruction of 3D neuron structures using a graphaugmented deformable model. Bioinformatics. 2010;26(12):i38–46.
 23.
Zhang Y, Jiang S, Xu Z, Gong H, Li A, Luo Q, Ren M, Li X, Wu H, Yuan J, et al. Pinpointing morphology and projection of excitatory neurons in mouse visual cortex. Front Neurosci. 2019;13:912.
 24.
Jiang S, Guan Y, Chen S, Jia X, Ni H, Zhang Y, Han Y, Peng X, Zhou C, Li A, et al. Anatomically revealed morphological patterns of pyramidal neurons in layer 5 of the motor cortex. Sci Rep. 2020;10(1):7916.
 25.
Wan Z, He Y, Hao M, Yang J, Zhong N. MAMST: an automatic 3D neuron tracing method based on mean shift and adapted minimum spanning tree. BMC Bioinformatics. 2017;18(1):197.
 26.
Peng H, Ruan Z, Long F, Simpson JH, Myers EW. V3D enables realtime 3D visualization and quantitative analysis of largescale biological image data sets. Nat Biotechnol. 2010;28(4):348–53.
 27.
Li Y, Gong H, Yang X, Yuan J, Jiang T, Li X, Sun Q, Zhu D, Wang Z, Luo Q, et al. TDat: an efficient platform for processing petabytescale wholebrain volumetric images. Front Neural Circuits. 2017;11:51.
 28.
Yen JC, Chang FJ, Chang S. A new criterion for automatic multilevel thresholding. IEEE T Image process. 1995;4(3):370–8.
 29.
Lei W, Jiao Y, Del Mar N, Reiner A. Evidence for differential cortical input to direct pathway versus indirect pathway striatal projection neurons in rats. J Neurosci. 2004;24(38):8289–99.
 30.
Blackman AV, Grabuschnig S, Legenstein R, Sjostrom PJ. A comparison of manual neuronal reconstruction from biocytin histology or 2photon imaging: morphometry and computer modeling. Front Neuroanat. 2014;8:65.
 31.
Cheng S, Wang X, Liu Y, Su L, Quan T, Li N, Yin F, Xiong F, Liu X, Luo Q, et al. DeepBouton: automated identification of singleneuron axonal Boutons at the brainwide scale. Front Neuroinform. 2019;13:25.
 32.
Meijering E. Neuron tracing in perspective. Cytom Part A. 2010;77(7):693–704.
Acknowledgements
We thank the MOST group members of Britton Chance Center for Biomedical Photonics for fMOST dataset and experimental assistance.
Funding
This work was financially supported through the National Natural Science Foundation of China (Grant Nos. 81827901, 61890954), the Science Fund for Creative Research Group of China (Grant No.61721092), the director fund of the WNLO.
Author information
Affiliations
Contributions
A.L. and S.C. originally conceived the study. S.J. and Z.P. proposed this method. S.J., A.L. and Y.G. draft the manuscript. Q.L., H.G., M.R. and Z.D. provide mouse brain dataset. S.J. and Z.F. performed image preprocessing. All authors discussed the conceptual and practical implications of the methods. All authors have read and approved the manuscript.
Corresponding author
Ethics declarations
Ethics approval and consent to participate
All the procedures followed in the animal experiments were approved by the Institutional Animal Ethics Committee of Huazhong University of Science and Technology. The animal care and use were done in accordance with the guidelines of the Administration Committee of Affairs Concerning Experimental Animals in Hubei Province of China.
Consent for publication
Not applicable.
Competing interests
The authors declare no financial or nonfinancial competing interests.
Additional information
Publisher’s Note
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Supplementary information
Additional file 1: Figure S1.
The topological graph of skeleton points. Figure S2. Quantitative performance of skeleton optimization method on synthetic datasets with different signal intensity. Figure S3. The cross sections of axonal and dendritic fibers of pyramidal neuron.
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, S., Pan, Z., Feng, Z. et al. Skeleton optimization of neuronal morphology based on threedimensional shape restrictions. BMC Bioinformatics 21, 395 (2020). https://doi.org/10.1186/s1285902003714z
Received:
Accepted:
Published:
Keywords
 Neuronal morphology
 Neuron tracing
 Neuronal skeleton optimization
 Shape restriction