Automatic and accurate reconstruction of long-range axonal projections of single-neuron in mouse brain

eLife Assessment

This important paper takes a novel approach to the problem of automatically reconstructing long-range axonal projections from stacks of images. The key innovation is to separate the identification of sections of an axon from the statistical rules used to constrain global structure. The authors provide compelling evidence that their method is a significant improvement over existing measures in circumstances where the labelling of axons and dendrites is relatively dense.

https://doi.org/10.7554/eLife.102840.3.sa0

Significance of the findings:

Important: Findings that have theoretical or practical implications beyond a single subfield

Landmark
Fundamental
Important
Valuable
Useful

Strength of evidence:

Compelling: Evidence that features methods, data and analyses more rigorous than the current state-of-the-art

Exceptional
Compelling
Convincing
Solid
Incomplete
Inadequate

During the peer-review process the editor and reviewers write an eLife Assessment that summarises the significance of the findings reported in the article (on a scale ranging from landmark to useful) and the strength of the evidence (on a scale ranging from exceptional to inadequate). Learn more about eLife Assessments

Abstract
Introduction
Results
Discussion
Materials and methods
Data availability
References
Article and author information
Metrics

Abstract

Single-neuron axonal projections reveal the route map of neuron output and provide a key cue for understanding how information flows across the brain. Reconstruction of single-neuron axonal projections requires intensive manual operations in tens of terabytes of brain imaging data and is highly time-consuming and labor-intensive. The main issue lies in the need for precise reconstruction algorithms to avoid reconstruction errors, yet current methods struggle with densely distributed axons, focusing mainly on skeleton extraction. To overcome this, we introduce a point assignment-based method that uses cylindrical point sets to accurately represent axons and a minimal information flow tree model to suppress the snowball effect of reconstruction errors. Our method successfully reconstructs single-neuron axonal projections across hundreds of GBs (Gigabytes) images within a mouse brain with an average of 80% f1-score, while current methods only provide less than 40% f1-score reconstructions from a few hundred MBs (Megabytes) images. This huge improvement is helpful for high-throughput mapping of neuron projections.

Introduction

Neuronal axons in general project to different brain regions, and their projection distribution is an essential cue for neuron type identification, neuronal circuit construction, and deeper insight into how information flows in the brain (Huang and Luo, 2015; Meijering, 2010; Parekh and Ascoli, 2013; Zingg et al., 2014). Advances in optical imaging and molecular labeling techniques (Cai et al., 2019; Chung and Deisseroth, 2013; Çiçek et al., 2016; Kim and Schnitzer, 2022; Li et al., 2010; Osten and Margrie, 2013) have allowed us to observe the entire mouse brain at single-axon resolution and provided the database for the study of neuronal projection patterns (Foster et al., 2021; Gao et al., 2022; Muñoz-Castañeda et al., 2021; Peng et al., 2021; Qiu et al., 2024; Sun et al., 2019; Xu et al., 2021; Zeng, 2022). However, the reconstruction of these long-range projected axons still requires extensive manual annotation in tens of TBs volumetric images (Çiçek et al., 2016; Friedmann et al., 2020; Wang et al., 2019; Winnubst et al., 2019; Zhou et al., 2021), this labor-intensive process creates a major bottleneck for high-throughput mapping of neuronal projections (Zeng and Sanes, 2017).

The difficulties in reconstructing the long-range projections of neurons are as follows. On the one hand, while molecular labeling techniques can shed light on a very small fraction of neurons, a significant fraction of neuronal axons is still densely distributed due to the morphological complexity of neurons. The identification of densely distributed axons is considered an open problem in the field (Li et al., 2019; Lichtman and Denk, 2011; Zeng and Sanes, 2017), which still has no good solution. On the other hand, during neuron reconstruction, reconstruction errors accumulate, and a single reconstruction error can result in an entire branch being connected erroneously to other neurons or missing (Helmstaedter, 2013). Therefore, effective large-scale reconstruction of neurons requires extremely high identification accuracy of dense axons. The contradictions between these two aspects seem hard to reconcile.

The current neuron reconstruction frameworks focus on how to accurately extract skeletons of neurites and establish the connections between skeletons (Meijering, 2010; Peng et al., 2015). The BigNeuron project (Manubens-Gil et al., 2023) conducts a systematic evaluation of 35 automatic neuron reconstruction algorithms, all of which are based on tracing neurite skeletons and can be divided into two categories: local and global approaches. In the local approach (Choromanska et al., 2012; Li et al., 2020; Peng et al., 2011; Yang et al., 2013), the localization of the next skeleton point requires computation of the signal anisotropy of the image region near the current skeleton point. Localization errors typically occur when this image region contains other neurite signals. The global approach (Li et al., 2019; Türetken et al., 2011; Xiao and Peng, 2013) first generates multiple seed points that are commonly located at the neurite centerline and then establishes connections between these seed points for generating the neurite skeleton. This connection relies mainly on spatial location information, resulting in densely distributed neurites being connected to each other erroneously. While deep learning is widely used in neuron reconstruction (Huang et al., 2020; Li and Shen, 2020; Liu et al., 2022; Zhou et al., 2018), - mainly for neuronal image segmentation and signal intensity enhancement to reduce reconstruction errors - even ideal segmentation with all neurite centers identified and their signal enhanced still exhibits significant reconstruction errors with skeleton-based methods (Figure 1—figure supplement 1).

To address the problem of error accumulation during neuron reconstruction, it is common practice to utilize statistical information of neuron morphology, such as the angle between two neurites, to identify and remove spurious connections between the reconstructed neurites. This strategy (Li et al., 2019; Quan et al., 2016) achieves 80% reconstruction accuracy from GB-scale images under two critical constraints: (1) precise identification of neurite terminals and branch points is required for accurate angle computation and morphological analysis, and (2) somatic locations are required as critical information to remove some links between the reconstructed neurites to ensure that each cell body can be mapped to the root node of a single tree structure. However, for long-range axonal reconstruction across hundreds of GB-scale images, the strategy is not effective to eliminate the accumulation of errors due to factors such as the position of the axon at a distance from the soma and slight morphological differences between axon junction and termination. Consequently, current long-range projection reconstruction methods are semi-automatic and require substantial human intervention (Gao et al., 2023; Wang et al., 2019; Winnubst et al., 2019; Zhou et al., 2021).

Here, we propose a new neuron reconstruction method called PointTree, which aims at how to assign foreground points in neuronal images to their own neurons. In the workflow, we design a constrained Gaussian clustering method to partition the foreground region of a neuronal image into a series of columnar regions whose centerline belongs to only a single neurite. This operation essentially eliminates the interference of different neurites in the dense reconstruction. In addition, each columnar region is characterized by a minimal envelope ellipsoid for constructing connections between columnar regions, which forms the neurite shapes. Based on the reconstructed shapes, we design a minimal information flow tree model to suppress the cumulative reconstruction error. Using the proposed method, we successfully achieve accurate reconstruction of long-range projections of neurons across hundreds of gigabytes of volumetric image.

Results

The architecture and principles of PointTree

In the design of PointTree, we have developed a series of optimization problems to assign foreground points in data blocks to their respective neurites. Firstly, the segment network is utilized for each data block to obtain foreground points. Subsequently, we apply a constrained Gaussian clustering method (Reynolds, 2009) to partition the foreground points into columnar regions and determine their geometrical parameters by solving the minimum-volume covering ellipsoids problem (Sun and Freund, 2004). Using these geometrical parameters, we construct a 0–1 assignment problem (Volgenant, 1996) to establish links between these columnar regions. Finally, skeletons are extracted from these linked columnar regions to reduce data redundancy by using region growing (Harris, 2011). The key procedures for neuron reconstruction are presented in Figure 1A.

Figure 1 with 4 supplements see all

Download asset Open asset

Summary and principle of PointTree.

(A) The reconstruction procedure of PointTree involves the generation, clustering, and connection of foreground points (the first row). Within this procedure, three optimization problems are designed to allocate the foreground points into their respective neurites (the second row). (B) Schematic diagram of information flow score calculation. In a neurite branch with a fixed root node (green circle), the information flow score is calculated based on the assumption that a neurite has few directional changes. The assumption determines the neurite directly connecting to the root node (red), resulting in two branch angles used to calculate the information flow score. (C) Statistical analysis of the consistency between the minimum information flow and the real situation. For 208 neurite branches, the information flow scores are calculated as ground truth according to their manually determined skeletons and root nodes. These scores are then displayed in ascending order. The root nodes of neurite branches are changed to generate both maximum and minimum information flow scores. (D) One neurite branch is decomposed into two by minimizing the total information flow scores. (E) Performance of different methods on separating closely paralleled neurites. In PointTree, a single neurite is represented by a series of ellipsoids whose centerlines are not simultaneously located within different neurites. They are connected using an ellipsoid shape, which results in perfect reconstruction (Left). However, skeleton-based methods fail to separate two closely paralleled neurites due to interference from other signals (Red circle in middle) or connections being interfered with by another neighboring skeleton point (Red circle in right).

In addition, PointTree employed the statistical prior information to reduce the reconstruction errors. At the branching point (node) of the neurites, it can be divided into three segments of neurite skeletons. The segment entering the node forms two angles with the other two segments exiting the node respectively. The node angle is defined as the smaller angle between the entering segment and each exiting segment (Figure 1B). With node angle, we can identify the single complete neurite and its corresponding node angles. The skeleton of the neurite is generally smooth, with very few sudden directional changes and even fewer at the nodes. So, the node angles should be as small as possible. For neuronal branches, the node angles are uniquely determined when the root node is given, and the sum of the negative cosine of these node angles expressed by information flow value is small when the root node is correctly identified. This rule is defined as a minimal information flow tree (MIFT).

In image blocks of densely distributed neurites, we used semi-automatic software (Zhou et al., 2021) extracting 208 neuronal branches and identifying their root nodes. For each branch, we calculated their information flow values as the ground-truth information flow values (Figure 1C). To validate MIFT, we looped through all possible structures of these branches by changing the root node in order to compute the maximum and minimum information flow values (Figure 1C). It is evident that, for most neuronal branches (195/208), the ground-truth values of the information flow achieve the minimum value, suggesting that MIFT rule is reasonable. We utilized MIFT to modify skeleton structure and remove spurious connections between reconstructed neurites (Figure 1D and Figure 1—figure supplement 2), both for reconstructions within individual blocks and for the fused reconstruction in adjacent blocks.

PointTree has the capability to separate densely distributed neurites. When dealing with two parallel neurites in close proximity to each other, their shapes can be represented by a series of columnar regions (the left panels of Figure 1E). We have modified the Gaussian clustering algorithm by constraining the estimated mean and covariance parameters so that the cluster shape approaches a columnar shape. Additionally, foreground points within the same cluster are connected to each other. These two features ensure that the central line in the columnar region belongs to only a single neurite, which is crucial for separating densely packed neurites. Furthermore, we utilize the minimum volume covering ellipsoid to extract shape information of the columnar regions for constructing their connections. These designs enable PointTree to successfully reconstruct packed neurites. In contrast, skeleton-based local methods rely on determining the position of the next skeleton point based on the shape anisotropy of the region. This often leads to localization errors when there are two neurite image signals within a region (the middle panels of Figure 1E). When it comes to skeleton-based global methods, although seed points can be located at individual neurite centers, accurately constructing connections between these seed points proves challenging due to the reliance on distance between points and susceptibility to interference from densely distributed neurites (the right panels of Figure 1E).

The merits of PointTree in dense reconstruction

In dense reconstruction, one of the main concerns is how well to separate densely distributed neurites that behave as crossover and closely paralleled neurites. These neurites can be manually identified by visualization with different view angles (Figure 2—figure supplement 1). We compared PointTree with several skeleton-based methods such as neuTube (Feng et al., 2015), PHDF (Radojevic and Meijering, 2017), NGPST (Quan et al., 2016), and MOST (Wu et al., 2014) in performing this task. We manually labeled the locations where neurites are crossover or closely parallel from five 256×256 × 256 image blocks. For a fair comparison, all methods are performed on segmented images derived from the segmentation network. Figure 2A illustrates the process of PointTree’s separation of crossover and closely paralleled neurites. PointTree can successfully separate the densely distributed neurites in a range of 71.4% and 91.7%, while these skeleton-based methods only separate 25.0% densely distributed neurites (Figure 2B) at most. We also present the comparison of PointTree and other methods on some reconstruction examples in which multi-crossover neurites (Figure 2C) and closely paralleled neurites are involved. PointTree provides the perfect reconstruction while other methods fail to reconstruct these neurites.

Figure 2 with 1 supplement see all

Download asset Open asset

Performance of PointTree on crossover and closely paralleled neurites.

(A) The reconstruction process of crossover and closely paralleled neurites. (B) Quantitative evaluation of PointTree and several skeleton-based methods on identifying closely distributed neurites. The box plots present the statistical information (n=5) in which the horizontal line in the box, the lower and upper borders of the box represent the median value, the first quartile (Q1), and the third quartile (Q3), respectively. The vertical black lines indicate 1.5 × IQR. (C) Three reconstruction examples derived from PointTree and several skeleton-based methods.

Furthermore, we present the quantitative results derived from PointTree and five widely used skeleton-based reconstruction methods, including APP2, neuTube, NGPST, PHDF, and MOST. Eight 256×256 × 256 image blocks that include many densely distributed neurites are of the testing dataset. All reconstruction algorithms are performed on the segmentation images of these testing datasets. We give the intuitive reconstruction comparisons (Figure 3A). PointTree provides the reconstruction close to the ground truth. The skeleton-based methods generate lots of reconstruction errors and incorrectly combine multi-neurites into a single branch. The quantitative reconstructions suggest that PointTree is far superior to skeleton-based methods (Figure 3B). For PointTree, the average precision is above 90%, both recall and f1-score are above 85%. The skeleton-based methods cannot provide a good solution to separate the densely packed neurites. The f1-score of these reconstructions ranges from 30% to 40%, which indicates the ineffective reconstructions.

Figure 3

Download asset Open asset

Comparison of reconstruction methods on image blocks containing densely distributed neurites.

(A) Comparison of reconstruction performance among six methods, including PointTree, NGPST, neuTube, APP2, PHDF, and MOST. Individual neurite branches are delineated in different colors. (B) Quantitative evaluation of reconstruction performance using precision, recall, and f1-score. The box plots display these three evaluation indexes (n=8). In the box, the horizontal line represents the median value. The box shows the interquartile range (IQR) from the first quartile (Q1) to the third quartile (Q3). The vertical lines indicate 1.5×IQR.

Reconstruction of data with different signal-to-noise ratios

In the field of neuronal reconstruction, data acquired by different imaging systems often exhibit varying signal-to-noise ratio (SNR) characteristics. For some low-SNR datasets, severe noise interference makes it difficult even for human observers to accurately identify neurite structures. To systematically evaluate PointTree’s reconstruction performance across datasets with different SNRs, we selected and analyzed data from three imaging systems: light sheet microscopy (Stelzer et al., 2021) (LSM), fluorescent micro-optical sectioning tomography (Wang et al., 2021) (fMOST), and high-definition fluorescent micro-optical sectioning tomography (Zhong et al., 2021) (HD-fMOST), with SNR ranges of 2–7, 6–12, and 9–14, respectively (Figure 4A).

Figure 4

Download asset Open asset

Reconstruction performance of PointTree across data with different signal-to-noise ratios.

(A) Data blocks from light sheet microscopy (LSM), fluorescent micro-optical sectioning tomography (fMOST), and high-definition fluorescent micro-optical sectioning tomography (HD-fMOST) are selected. SNR and corresponding reconstruction scores with PointTree are drawn with line charts. Each dataset is of sample size n=25 and each data block size of 128×128 × 128. (B) shows reconstruction performance of PointTree on different datasets. (C) The zoomed-in view displays the region marked by white box in the first column of (B), with 25 foreground points and 25 background points sampled respectively. The signal intensities of both the foreground points and background points are plotted in the adjacent line charts.

Experimental results demonstrate that, thanks to the powerful feature extraction capability of the deep learning network, the trained neural network achieves satisfactory segmentation performance (third row in Figure 4B) even on low-SNR data (first two columns in Figure 4B, top row), laying a solid foundation for subsequent accurate reconstruction (bottom row in Figure 4B). Quantitative analysis reveals that PointTree delivers stable reconstruction performance across all SNR levels. Specifically: for LSM data (sample size n=25, mean SNR = 5.01), average precision = 96.0%, recall = 88.7%, and f1-score=91.0%; for fMOST data (sample size n=25, mean SNR = 8.68), average precision = 95.8%, recall = 87.3%, and f1-score=90.0%; for HD-fMOST data (sample size n=25, mean SNR = 11.4), average precision = 98.1%, recall = 91.0%, and f1-score=93.3% (Figure 4A).

Notably, in low-SNR LSM data, background regions contain more artifactual signals (first panel in Figure 4C) due to similar intensity distributions between background and foreground points. In contrast, high-SNR datasets (fMOST and HD-fMOST) exhibit cleaner background features with distinct intensity separation between background noise and neurite signals (second and third panel in Figure 4C). This observation highlights the critical impact of SNR on reconstruction quality while simultaneously validating the robustness of PointTree, which is aided by the segmentation network, across diverse SNR conditions.

Restrain error accumulation in the reconstruction

In order to achieve accurate axon reconstruction, it is essential to effectively suppress the snowballing accumulation of reconstruction errors. The performance of the minimal information flow tree (MIFT) in retraining the reconstruction errors is evaluated in this study. Figure 5A presents six 512×512 × 512 image blocks and their reconstructions using PointTree in the first column. The reconstruction fusing procedure is then performed on these axonal reconstructions (Figure 5A). By employing MIFT to revise the reconstructions and remove false connections between axons, reasonable reconstructions are achieved. In contrast, when the same fusion procedure is conducted without MIFT to revise the reconstruction, almost all axons are incorrectly connected together (bottom-right panel in Figure 5A).

Figure 5

Download asset Open asset

Minimal information flow tree effectively restrains the accumulation of reconstruction errors.

(A) Reconstruction comparisons in the fusion process with MIFT and without MIFT are shown. Both image blocks and neurite reconstructions are displayed using maximum projection along the z-direction. Two fusion procedures are performed, and the final fusion reconstructions are presented in the third column. (B) The variation in reconstruction accuracy during the fusion process with MIFT and without MIFT is illustrated. Blue points represent the initial reconstruction accuracy from six image blocks, while green points and red points denote the merged reconstruction accuracy with MIFT and without MIFT, respectively. The squares represent the mean values of the evaluation indexes. (C) The skeletons of three neurite branches from the final merged reconstructions with MIFT are shown. Additionally, corresponding ground-truth reconstructions and reconstruction evaluations are also presented.

We furthermore measure the enhancement in the reconstruction accuracy achieved by MIFT (Figure 5B). For the initial reconstructions from six image blocks, the average of f1-score is about 0.86. By using MIFT, the average of f1-score is above 0.8 for the reconstructions from two image blocks which are generated with the first fusion. In the second fusion (top-right panel in Figure 5A), the f1-score still keeps 0.79. In contrast, without MIFT, the first fusion leads to a drop of about f1-score of 0.3. After the second fusion, the f1-score is less than 0.2. We also present some reconstruction examples after two fusions in Figure 5C, which are close to the ground truth. These results suggest that the MIFT model takes consideration of the proper structure of axons and thus can restrain the error communications in the reconstruction fusion process.

Long-range axonal projections reconstruction

We applied PointTree for long-range axon reconstruction. The testing image block has the size of 11226×8791 × 1486 voxels and includes axons from eight neurons (Figure 6A). We also used GTree to manually reconstruct these neurons as the ground-truth reconstruction (Figure 6B). Except for the labeling of training data for segmentation network and of the axon starting points of a single neuron, the whole reconstruction process is totally automatic. The results show PointTree successfully recovered the axonal morphology of these eight neurons without manual interference (Figure 6C and Videos 1 and 2), and we compared these reconstructions with ground truth (Figure 6—figure supplement 1). The average precision is above 85% and the average recall and f1-score are above 80% (Figure 6E). In addition, we presented the axon reconstructions from two image blocks (Figure 6C1 and C2) which include a large number of densely distributed axons. This reconstruction performance suggests that the point assignment and the minimal information flow tree mode, as the two key strategies in PointTree, perform well in long-range axonal reconstruction.

Figure 6 with 2 supplements see all

Download asset Open asset

Long-range axonal reconstruction using PointTree.

(A) The image block contains eight neurons in the ventral posteromedial thalamic region. The projection of these neurons includes a large number of densely distributed axons, which are enlarged in A₁ and A₂. (B) The reconstruction of the eight neurons is achieved by annotators with semi-automatic software GTree, serving as ground-truth reconstruction to evaluate automatic algorithms. The reconstructions B₁ and B₂ correspond to the image blocks A₁ and A₂. (C) Automatic reconstruction with PointTree results in reconstructions of the densely distributed axons, which are enlarged in C₁ and C₂. (D) A comparison between automatic reconstruction and ground-truth reconstruction of axonal projection for one neuron is shown. Green indicates consistent reconstruction, blue indicates missed branches, and red denotes branches from other neurons. (E) Quantitative analysis of long-range projections for these neurons is presented. Statistical information is displayed in boxes (n=8), the horizontal line in the box, the lower and upper borders of the box represent the median value, the first quartile (Q1) and the third quartile (Q3) respectively, the vertical black lines indicate 1.5 × IQR, while black points represent the accuracy of the reconstructions for these neurons.

Video 1

Download asset

posterframe for video — Reconstructed long-range axonal projections and raw image data shown in Figure 6, individual axonal projections are delineated in different colors.

Video 2

Download asset

We also applied PointTree to process another 10739×11226 × 3921 image blocks collected with HD-fMOST system (Zhong et al., 2021). The high signal-to-noise ratio in this optical system results in a significantly extended dynamic range of the signal. PointTree can effectively deal with this case, and all 14 long-range projections are successfully reconstructed (Figure 6—figure supplement 2). The quantitative results suggest that the average f1-score is above 90% (Table 1).

Table 1

Quantitative metrics comparing ground truth and reconstructed neurons are presented in Figure 6—figure supplement 2.

ID	Precision	Recall	F1-Score
1	1.00	0.92	0.95
2	1.00	1.00	1.00
3	0.98	0.76	0.86
4	1.00	0.82	0.90
5	1.00	0.77	0.87
6	1.00	0.92	0.96
7	0.96	0.75	0.84
8	1.00	0.87	0.93
9	1.00	0.82	0.90
10	1.00	0.96	0.98
11	1.00	0.99	0.99
12	1.00	0.77	0.87
13	1.00	0.90	0.95
14	0.99	0.87	0.93

Despite the need to solve multiple large-scale optimization problems, the reconstruction speed using PointTree is generally faster than the imaging speed. For instance, in a typical scenario involving 254 image blocks with 512×512 × 512 voxels, the total time required for reconstruction is approximately 44 min. Even for a larger dataset comprising 821 image blocks with 512×512 × 512 voxels and including a significant number of sparsely distributed neurites, the total time cost amounts to about 60 min (Table 2). It should be noted that the time cost does not increase linearly as data volume increases due to the influence of neurite density on overall reconstruction time. In summary, PointTree demonstrates remarkable speed in reconstructing long-range axons (Video 3).

Table 2

Time cost of three modules in the entire reconstruction for two testing datasets shown in Figure 6, Figure 6—figure supplement 2.

block number(size: 512×512 × 512)	Points clustering(mins)	Clusters connection(mins)	Reconstruction merging (mins)
254	23	18	3
821	22	35	3

Video 3

Download asset

Discussion

We have presented an automated method for reconstructing the long-range projections of neurons. In this study, we address the problem of mutual interference among densely distributed neurites and the cumulative error during reconstruction by designing a reconstruction method based on point set assignment and the minimal information flow tree, respectively. As a result, our approach enables accurate reconstruction of long-range neuron projections from hundreds of gigabytes of data. This advance significantly enhances the efficiency of whole-brain-scale neuron reconstruction, bridging the substantial gap between factory-level generation of whole-brain-scale neuronal imaging data and tens of hours required to reconstruct one neuron.

Our approach is performed on image foregrounds where the segmented neurites have a fixed radius approximately equal to the total size of the three voxels. In this case, we can estimate the total number of foreground points (voxels) and set a suitable number of columnar regions for ensuring the anisotropy of each columnar region, which is based on the fact that the union of columnar regions equals the foreground region. The anisotropy of the columnar regions will reduce the difficulty in establishing their connection. The requirement that all segmented neurites have a relatively fixed radius can be fulfilled. For all neurites, the value of their voxels decreases as these voxels deviate from the nearest centerline. The deep learning network is able to grasp this feature and segment only the neurite centerline and its neighborhood. Typically, in reconstructions of neurons whose projections are distributed over hundreds to thousands of GBs of data, less than GB-sized images with labels are needed as training data. The labeling process takes a few hours, which is negligible for semi-automatic reconstruction of all neurons in the whole volume images.

We propose a new reconstruction mode centered on point set assignment instead of the current reconstruction mode focused on skeleton extraction. In the current reconstruction paradigm, most deep networks are used to enhance the signal-to-noise ratio of neuronal images and do not address the issue of signal interference during skeleton extraction. In contrast, our reconstruction approach is based on directly processing the foreground points generated by the deep learning network. With continued advances in deep learning techniques, the generality and accuracy of image segmentation will be continuously enhanced, thereby significantly boosting the application scope of our method in various scenarios. Essentially, our method can be applied to any skeleton tracking-based application scenario and effectively eliminate dense signal interference.

Our method still generates a few reconstruction errors. This is due to the following three aspects. First, our method directly handles image foregrounds, which leads to reconstruction errors when some neurites with weak image intensities are not identified. Second, relying solely on foreground point information and rule-based judgment methods may generate some connection errors when establishing connections between neurites. Finally, the minimal information flow tree’s fundamental assumption, that axons should be as smooth as possible, does not always hold true. In fact, real axons can take quite sharp turns (Figure 1—figure supplement 3) leading the algorithm to erroneously separate a single continuous axon into disjoint fibers (Figure 1—figure supplement 3). Therefore, for the automatic reconstruction of neurons on a brain-wide scale, further work is needed to enhance the imaging intensity and incorporate soma shapes and raw image signals for neurites connection recognition.

Materials and methods

Data collections

Request a detailed protocol

All animal experiments followed procedures approved by the Institutional Animal Ethics Committee of the Huazhong University of Science and Technology. The test datasets are collected through the preparation of two kinds of samples. For one C57BL/6 male mouse, 100 nl AAV-Cre virus and 100 nl of AAV-EF1α-DIO-EYFP virus were injected into the VPM nucleus at the same time. 21 days later, the chemical sectioning fluorescence tomography (CSFT) system (Wang et al., 2021) was used to acquire imaging data (Figures 1—6), more details can be seen in the reference (Zhang et al., 2021). For one C57BL/6 J male mouse, 100 nl of AAV-YFP was injected into the motor area. 21 days later, high-definition fluorescent micro-optical sectioning tomography (HD-fMOST) was used to acquire imaging data (Zhong et al., 2021; Figure 6—figure supplement 2).

Generation of foreground points

Request a detailed protocol

Our reconstruction method performs on the image foregrounds. Here, we used UNet3D (Çiçek et al., 2016) for image stacks segmentation without network structure modification. The detailed information about UNet3D can be found in the reference (Çiçek et al., 2016). Considering the requirement that the network output, the segmented neurites, have the relatively fixed radius, we calculate the distance field of the neurite’s skeleton as the ground truth for supervising the network. Initially, the semi-automatic software GTree was utilized to extract the neurite skeleton and subsequently interpolate the skeleton points. The interpolation operation ensured that the distance between any skeleton point and its nearest point was less than 1 μm. Subsequently, the interpolated skeleton points were used as centers to mark spherical regions with a radius of 5 voxels. These spherical regions served as candidate areas for foreground. Within these candidate areas, the distance from each point to its nearest interpolated skeleton point was calculated. Finally, the distances are mapped into Gaussian kernel distances, which form the Gaussian density map. This map normalized by maximum value leads to the distance field map to supervise UNet3D output.

In the training stage, Adam optimizer is used with an initial learning rate at 3e-4. The input image size is 128×128 × 128. Batch size is set to 1, the L1-norm is used as loss function to train the network. We presented the reconstructions from two kinds of fMOST datasets. One is from the reference (Zhang et al., 2021) and the other is from the reference (Zhong et al., 2021). Therefore, we created two sets of training data, each consisting of 20 512×512 × 512 image blocks (each divided into 64 image blocks of size 128×128 × 128). In each set, 10 image blocks contain densely distributed neurites, while the other 10 blocks contain sparsely distributed neurites. In the predicting stage, we applied the threshold operation to the distance field image. The voxels whose values are more than 0.5 are regarded as the foreground points.

Neuron Reconstruction based on Points assignment

Request a detailed protocol

For the image stack, we allocated the foreground points to their respective neurites and established connections between neurites by constructing three optimization models: (1) the constrained Gaussian mixture model divides the foreground points into a set of points, each of which has a column shape; (2) the minimum-volume covering ellipsoids model extracts the features of the column-shaped point set; (3) the 0–1 assignment optimization model establishes connections between the column-shaped point sets, resulting in the shapes of individual neurites, and then builds connections between the reconstructed neurites.

Constrained Gaussian mixture model

Request a detailed protocol

The three-dimensional Gaussian function exhibits an ellipsoidal shape in space, which we have utilized to approximate the columnar shape of local neurites. In this study, Gaussian distribution mixture functions with $K$ components are employed to approximate the shape of all neurites in an image block. The component number $K$ is obtained by point density and will be discussed later. Given the foreground points $x_{1}, x_{2}, \dots, x_{n}$ , for each foreground points $x_{i}$ , the probability density function $P (x_{i})$ is calculated as follows:

P (x_{i}) = \sum_{j = 1}^{K} π_{j} N (x_{i} | μ_{j}, Σ_{j})

Here, $N (x_{i} | μ_{j}, Σ_{j})$ is the Gaussian density function with mean value $μ_{j}$ and covariance matrix $Σ_{j}$ . Weight $π_{j}$ is the regularization parameter. $N (x_{i} | μ_{j}, Σ_{j})$ is given by the formula:

N (x_{i} | μ_{j}, Σ_{j}) = \frac{1}{2 π^{3 / 2} | Σ_{j} |^{1 / 2}} e^{- \frac{1}{2} {(x_{i} - μ_{j})}^{T} Σ_{j}^{- 1} (x_{i} - μ_{j})}

Based on probability density function, the conditional probability can be computed as:

p_{i, j} = P (x_{i} | c l u s t e r_{j}) = \frac{π_{j} N (x_{i} | μ_{j}, Σ_{j})}{\sum_{j = 1}^{K} π_{j} N (x_{i} | μ_{j}, Σ_{j})} j = (1, 2, . . ., K)

Here, $p_{i, j}$ is the conditional probability for $x_{i}$ to assign to the j-th cluster. If $p_{i, k}$ is the maximum value among ${p_{i, 1}, . . . p_{i, K}}$ , the foreground point $x_{i}$ will be assigned to the k-th cluster. All the points assigned to the k-th cluster form a columnar region. Considering that both the number of foreground points and component number are large, we have added some constrained conditions for Gaussian mixture model as follows:

\sum_{j = 1}^{K} π_{j} = 1

I (μ_{j}) \geq ε_{0}, | Σ_{j} | \leq ε_{1}

$\sum_{j = 1}^{K} π_{j} = 1$ refers to the fact that the total probability distribution normalizes to 1. $I (\cdot)$ represents the signal intensity from segment image, $ε_{0}$ is the minimum signal intensity of foreground points and is set to 128 in the algorithm. $I (μ_{i}) \geq ε_{0}$ restrain the center of the Gaussian distribution to be a foreground point. $| Σ_{j} | \leq ε_{1}$ restrain the determinant of the covariance matrix which controls the suitable number of foreground points for each columnar region. $ε_{1}$ is set to the cube of three times the average diameter of neurite.

Maximum likelihood is employed to estimate the parameters of Gaussian mixture model and the final optimization problem is formed as follows:

{(π_{j}^{*}, μ_{j}^{*}, Σ_{j}^{*})}_{j = 1, 2, \dots, K} = \arg max \prod_{i = 1}^{n} P (x_{i}) = \arg max \prod_{i = 1}^{n} (\sum_{j = 1}^{K} π_{j} N (x_{i} | μ_{j}, Σ_{j}))

s . t . \sum_{j = 1}^{K} π_{j} = 1, I (μ_{j}) \geq ε_{0}, | Σ_{j} | \leq ε_{1}

In solving this optimization problem, we employ peak density algorithm (Wei et al., 2023) to compute density for each foreground points and sort them in descending order. We first select a point as a seed point, and the foreground points within a radius of 5 centered on it will be excluded. Then we continue selecting seed points until all foreground points are either selected or excluded. The selected $K$ seed points represent the initial $K$ components. We select signal points from the median (based on density) to both sides as seed points, which can decrease the situations that seed points lie in the center of a crossover or the edge of neurites. This strategy can make the generated columnar regions be more reasonable. The positions of the $K$ seed points are set to the initial $(μ_{1}, μ_{2}, \dots, μ_{K})$ . The initial setting of the covariance matrix is the identity matrix. The constrained Gaussian mixture model was solved by the EM algorithm (McLachlan and Krishnan, 2007), the EM algorithm is divided into two steps:

E-step: For each point $x_{i}$ , compute its probability within each Gaussian distribution using the probability density function:

p_{i, j} = \frac{π_{j} N (x_{i} | μ_{j}, Σ_{j})}{\sum_{j = 1}^{K} π_{j} N (x_{i} | μ_{j}, Σ_{j})}

M-step: Update the mean value, covariance matrices, and weight vectors.

π_{j} = \frac{\sum_{i = 1}^{n} p_{i, j}}{n}

μ_{j} = \frac{\sum_{i = 1}^{n} p_{i, j} x_{i}}{\sum_{i = 1}^{n} p_{i, j}}

Σ_{j} = \frac{{\sum_{i = 1}^{n} p_{i, j} (x_{i} - μ_{j}) (x_{i} - μ_{j})}^{T}}{\sum_{i = 1}^{n} p_{i, j}}

Besides, the constrained Gaussian mixture model possesses additional constraints: $I (μ_{j}) \geq ε_{0}$ and $| Σ_{j} | \leq ε_{1}$ . After finishing the M-step, $μ_{j}$ with $I (μ_{j}) < ε_{0}$ are selected. Eigenvalue decomposition is applied on $Σ_{j}$ and obtains eigenvalues $(γ_{1}, γ_{2}, γ_{3})$ in descending order and eigenvectors $(v_{1}, v_{2}, v_{3})$ . $μ_{j}$ is updated along $v_{1}$ and $- v_{1}$ to generate two new clusters with mean value and covariance matrices $(u_{j, 1}, Σ_{j, 1})$ and $(u_{j, 2}, Σ_{j, 2})$ as follows:

u_{j, 1} = u_{j} + v_{1} \cdot \frac{γ}{2}

u_{j, 2} = u_{j} - v_{1} \cdot \frac{γ}{2}

Σ_{j, 1} = \frac{\sum_{i = 1}^{n} p_{i, j} (x_{i} - μ_{j, 1}) {(x_{i} - μ_{j, 1})}^{T}}{\sum_{i = 1}^{n} p_{i, j}}

Σ_{j, 2} = \frac{\sum_{i = 1}^{n} p_{i, j} (x_{i} - μ_{j, 2}) {(x_{i} - μ_{j, 2})}^{T}}{\sum_{i = 1}^{n} p_{i, j}}

For $Σ_{j} > ε_{1}$ , it will be updated as follows:

Σ_{j}^{'} = \frac{ε_{1}}{Σ_{j}} Σ_{j}

Iteration of E-step and M-step will continue until the k-th result ${μ^{k}, Σ^{k}}$ and (k-1)-th result satisfy the stopping criteria:

‖ \frac{u^{k} - u^{k - 1}}{u^{k - 1}} ‖ < ε a n d ‖ \frac{Σ^{k} - Σ^{k - 1}}{Σ^{k - 1}} ‖ < ε

Here the division represents element-wise division and $‖ \cdot ‖$ denotes $L_{2}$ -norm and $ε$ is set to 0.01.

Shape characterization of columnar regions

Request a detailed protocol

After deriving the columnar regions through solving the constrained Gaussian mixture model, it is imperative to characterize their geometric shape (terminals and centerlines). For this purpose, we calculate the minimum-volume ellipsoids that can fully encompass each individual columnar region. For $c \in R^{3}$ , $Q \in S_{+ +}^{3}$ , a three-dimensional ellipsoid can be defined as follows Sun and Freund, 2004:

E_{c, Q} : = {x \in R^{3} | {(x - c)}^{T} Q (x - c) \leq 1}

Here, $c$ is the center of ellipsoid, $Q$ represents the geometric shape, $S_{+ +}^{3}$ denotes the convex cone of 3×3 symmetric positive definite matrices. The volume of $E_{c, Q}$ is given by the formula:

V o l u m e (E_{c, Q}) = \frac{π^{3 / 2}}{Γ (3 / 2 + 1)} \frac{1}{\sqrt{d e t (Q)}}

Here, $Γ (\cdot)$ is the standard gamma function of calculus, $d e t (Q)$ means the determinant of matrix Q. Minimizing the volume of $E_{c, Q}$ is equivalent to minimizing $\det (Q^{- 1 / 2})$ . Therefore, for a columnar region with foreground points $P {x_{1}, x_{2}, \dots x_{m}}$ , we define the target function as follows:

P 1 : (c^{*}, Q^{*}) = {\arg min}_{c, Q} d e t (Q^{- 1 / 2})

s . t . {(x_{i} - c)}^{T} Q (x_{i} - c) \leq 1, i = 1, 2 . . . m

c \in C H u l l (P), Q \in S_{+ +}^{3}

Here $c \in C H u l l (P_{i})$ restrain the solved center of ellipsoid to locate within the smallest convex hull formed by the clustering points. To solve this problem, a variable substitution $A = Q^{1 / 2}$ and $y = Q^{1 / 2} c$ were applied to Equation 20 and Equation 21, the original problem P1 can be transformed into a convex optimization problem as follows:

P 2 : (A^{*}, y^{*}) = \underset{A, y}{\arg min} - l n d e t (A)

s . t . {(A x_{i} - y)}^{T} (A x_{i} - y) \leq 1, i = 1, 2, . . ., m

A \in S_{+ +}^{3}

Through adding the logarithmic barrier function, we can obtain the following formula:

P 3 : (A^{*}, y^{*}, θ^{*}) = \underset{A, y, θ}{\arg min} - l n d e t (A) - θ \sum_{i = 1}^{m} \ln (z_{i})

s . t . {(A x_{i} - y)}^{T} (A x_{i} - y) + z_{i} = 1, i = 1, 2, . . ., m

A \in S_{+ +}^{3}, z_{i} > 0

As $θ$ varies in the interval $(0, \infty)$ , the solution of $P 3$ changes. When $θ$ approaches 0, the optimal solution of $P 3$ tends to the optimal solution of $P 2$ . By adding the dual multipliers $d_{i}$ which satisfies $d_{i} \cdot z_{i} = θ$ , the optimality conditions can be written as:

\sum_{i = 1}^{m} d_{i} [(A x_{i} - y) x_{i}^{T} + x_{i} {(A x_{i} - y)}^{T}] - A^{- 1} = 0

\sum_{i = 1}^{m} d_{i} (y - A x_{i}) = 0

{(A x_{i} - y)}^{T} (A x_{i} - y) + z_{i} = 1 i = 1, 2, . . ., m

\sum_{i = 1}^{m} d_{i} \cdot z_{i} = θ, i = 1, 2, . . ., m

d_{i}, z_{i} \geq 0

At this point, the error between the solution of the system of equations and the optimal solution of $P 3$ is less than $d^{T} z$ . Through Equation 30, the explicit expression for solving $y$ can be obtained as follows:

y = \frac{A X d}{e^{T} d}

Here, $X$ stands for a $3 \times m$ matrix $[x_{1} | x_{2} | ... | x_{m}]$ , $e$ stands for vector of ones ${(1, 1, ..., 1)}_{1 \times m}^{T}$ and $d$ stands for ${(d_{1}, d_{2}, ..., d_{m})}_{1 \times m}^{T}$ . Substitute Equation 34 into Equation 29, the equation for matrix $A$ can be obtained by:

(X D X^{T} - \frac{X d d^{T} X^{T}}{e^{T} d}) A + A (X D X^{T} - \frac{X d d^{T} X^{T}}{e^{T} d}) = A^{- 1}

Here, $D$ stands for a $m \times m$ diagonal matrix $D i a g (d_{1}, d_{2}, ..., d_{m})$ . And the explicit expression for $A$ is formed as

A = A (d) = {[2 (X D X^{T} - \frac{X d d^{T} X^{T}}{e^{T} d})]}^{- 1 / 2}

And explicit expression for $y$ :

y = \frac{{[2 (X D X^{T} - \frac{X d d^{T} X^{T}}{e^{T} d})]}^{- 1 / 2} X d}{e^{T} d}

Through substituting the above two equations to the system of Equations 29-33, variables A and y are eliminated. The following system of equations with only variables d and z can be obtained:

f (d) + z - e = 0

D z - θ e = 0

d_{i}, z_{i} \geq 0

Here, $f (d)$ is nonlinear function of variable $d$ :

f_{i} (d) = (x_{i} - \frac{X d}{e^{T} d}) {[2 (X D X^{T} - \frac{X d d^{T} X^{T}}{e^{T} d})]}^{- 1} \cdot (x_{i} - \frac{X d}{e^{T} d}) i = 1, 2, . . ., m

For a fixed barrier parameter $θ$ , we employ Newton’s method to solve the system of equations. We use $\nabla_{d} f (d)$ to represent the Jacobian matrix of $f (d)$ . Thus, the Jacobian matrix of the system of equations can be computed as follows:

[\begin{array}{cc} \nabla_{d} f (d) & I \\ Z & D \end{array}]

And the Newton’s direction is written as:

Δ (d) = {(\nabla_{d} f (d) - D^{- 1} Z)}^{- 1} (h_{1} - D^{- 1} h_{2})

Δ (z) = D^{- 1} h_{2} - D^{- 1} Z {(\nabla_{d} f (d) - D^{- 1} Z)}^{- 1} (h_{1} - D^{- 1} h_{2})

h_{1} = e - z - f (d), h_{2} = θ e - D z

With initial $(d_{0}, z_{0})$ , iterate with $(d_{n}, z_{n}) = (d_{n - 1}, z_{n - 1}) + \tilde{β} (Δ (d_{n - 1}), Δ (z_{n - 1}))$ to obtain the final optimal solution, $\tilde{β}$ represents the Newton’s step. Detailed process can see the pseudo code as follows:

Algorithm 1. Compute Newton’s direction.
Input: $(d, z, θ)$ satisfying $d, z > 0$ , $θ \geq 0$ 1. $A^{- 2} (d) = [2 (X D X^{T} - \frac{X d d^{T} X^{T}}{e^{T} d})]$ 2. $Σ (d) = (X - \frac{X d e^{T}}{e^{T} d}) A^{2} (d) (X - \frac{X d e^{T}}{e^{T} d})$ 3 $\nabla_{d} f (d) = - 2 (\frac{Σ (d)}{e^{T} d} + Σ (d) \circ Σ (d))$ 4. $(Δ (d), Δ (z)) = ({(\nabla_{d} f (d) - D^{- 1} Z)}^{- 1} (h_{1} - D^{- 1} h_{2}), D^{- 1} h_{2} - D^{- 1} Z Δ (d))$ Output: $(Δ (d), Δ (z))$

Algorithm 1. Compute Newton’s direction.

Input:

(d, z, θ)

satisfying

d, z > 0

θ \geq 0

A^{- 2} (d) = [2 (X D X^{T} - \frac{X d d^{T} X^{T}}{e^{T} d})]

Σ (d) = (X - \frac{X d e^{T}}{e^{T} d}) A^{2} (d) (X - \frac{X d e^{T}}{e^{T} d})

\nabla_{d} f (d) = - 2 (\frac{Σ (d)}{e^{T} d} + Σ (d) \circ Σ (d))

(Δ (d), Δ (z)) = ({(\nabla_{d} f (d) - D^{- 1} Z)}^{- 1} (h_{1} - D^{- 1} h_{2}), D^{- 1} h_{2} - D^{- 1} Z Δ (d))

Output:

(Δ (d), Δ (z))

Algorithm 2. Process of solving P2.
Input: ${x_{1}, x_{2}, . . ., x_{m}}$ 1. $r = 0.99$ , $(d_{0}, z_{0}) = (\frac{3}{2 m} e, e - f (d_{0}))$ 2. $E = - det (A (d))$ 3. while ( $\| e - f (d) - z \| > ε_{1}$ or $\frac{d^{T} z}{E} > ε_{2}$ ) 4. $θ = \frac{d^{T} z}{10 m}$ 5. $(Δ (d), Δ (z))$ = Compute_Newton_direction $(d, z)$ 6. $\bar{β} = max {β \| (d, z) + β (Δ (d), Δ (z) \geq 0)}$ 7. $\tilde{β} = min (\bar{r β}, 1)$ 8. $(d, z) = (d, z) + \tilde{β} (Δ (d), Δ (z))$ 9. $E = - det (A (d))$ Output: $Q = A {(d)}^{2}, c = A {(d)}^{- 1} y (d)$

Algorithm 2. Process of solving P2.

Input:

{x_{1}, x_{2}, . . ., x_{m}}

r = 0.99

(d_{0}, z_{0}) = (\frac{3}{2 m} e, e - f (d_{0}))

E = - det (A (d))

3. while (

| e - f (d) - z | > ε_{1}

\frac{d^{T} z}{E} > ε_{2}

)

4.

θ = \frac{d^{T} z}{10 m}

(Δ (d), Δ (z))

= Compute_Newton_direction

(d, z)

\bar{β} = max {β | (d, z) + β (Δ (d), Δ (z) \geq 0)}

\tilde{β} = min (\bar{r β}, 1)

(d, z) = (d, z) + \tilde{β} (Δ (d), Δ (z))

E = - det (A (d))

Output:

Q = A {(d)}^{2}, c = A {(d)}^{- 1} y (d)

With the solved optimal solution of $(Q, c)$ , we then check whether $c$ is located within the convex hull of the input point set ${x_{1}, x_{2}, . . ., x_{m}}$ . If it is not, a constrained Gaussian mixture model will be applied to partition it into two subsets and solve the minimum-volume covering ellipsoids problem again in the two subsets. Through solving the above minimum-volume covering ellipsoids problem, we can characterize the columnar regions more accurately.

Note that from constrained GMM, each cluster has the corresponding mean and covariance matrix of points in the cluster. These two values essentially describe the shape of the cluster. However, if these two values directly replace $c^{*}$ and $Q^{*}$ , the exported ellipsoid may only encompass a part of points in the cluster. For covering all points in the cluster, all elements in the covariance matrix are needed to be proportionally enlarged, but the volume of the corresponding ellipsoid is not minimum. These two cases will reduce the accuracy of the connections between clusters, that is columnar regions. So, we introduce the minimum-volume covering ellipsoid model to extract the shape of columnar region.

Skeleton generation using 0-1 assignment model

Request a detailed protocol

The 0–1 assignment model (Volgenant, 1996) can robustly and accurately establish connections between particles in live-cell imaging (Jaqaman et al., 2008). It is particularly effective in handling cases where particles are densely distributed, merged, or split. We analogize column regions to particles and apply the 0–1 assignment model to build the connections between column regions. For the i-th columnar region, the center and the two endpoints of the longest axis of its minimum-volume covering ellipsoid are denoted by $c_{i}, t_{i, 0}, t_{i, 1}$ . The direction refers to the pointing of the center point towards $t_{i, k}$ , k equal to 0 or 1. According to the direction and the endpoints, we design the cost matrix for building the 0–1 assignment model.

C = {[\begin{array}{cccccc} c (t_{1, 0}, t_{1, 0}) & c (t_{1, 0}, t_{1, 1}) & \dots & c (t_{1, 0}, t_{n, 1}) \\ c (t_{1, 1}, t_{1, 0}) & c (t_{1, 1}, t_{1, 1}) & \dots & c (t_{1, 1}, t_{n, 1}) & D \\ ⋮ & ⋮ & ⋮ & ⋮ \\ c (t_{n, 1}, t_{1, 0}) & c (t_{n, 1}, t_{1, 1}) & \dots & c (t_{n, 1}, t_{n, 1}) \\ D & D \end{array}]}_{4 n \times 4 n}

c (t_{i, i 0}, t_{j, j 0}) = {\begin{cases} 100 & i f (i = j) \\ \frac{n o r m (t_{i, i 0}, t_{j, j 0})}{{(0.5 \times (\frac{θ (t_{i, i 0}, t_{j, j 0})}{3} + 1.001))}^{4}} & i f (i \neq j) \end{cases}

\begin{aligned} θ (t_{i, i 0}, t_{j, j 0}) = & ⟨ d i r (c_{i}, t_{i, i 0}), d i r (c_{i}, t_{j, j 0}) ⟩ + ⟨ d i r (c_{j}, t_{j, j 0}), d i r (c_{j}, t_{i, i 0}) ⟩ \\ - ⟨ d i r (c_{i}, t_{i, i 0}), d i r (c_{j}, t_{j, j 0}) ⟩ \end{aligned}

Here, D is 2n×2n auxiliary matrix all elements of which are all set 100. Both $i 0$ and $j 0$ in Equation 47 are equal to 0 or 1, labeling the two endpoints of the longest axis of the ellipsoid. $n o r m (t_{i, i 0}, t_{j, j 0})$ represents the Euclidean distance between $t_{i, i 0}$ and $t_{j, j 0}$ . $θ (t_{i, i 0}, t_{j, j 0})$ describes the angle between two ellipsoids, that is two columnar regions. $d i r (c_{i}, t_{i, i 0})$ represents the line from point $c_{i}$ to $t_{i, i 0}$ . $〈 d i r (c_{i}, t_{i, i 0}), d i r (c_{i}, t_{j, j 0}) 〉$ represents cosine angle between the two lines. The threshold of 100 in D in Equation 46 and Equation 47 is an experimental value designed to ensure that the terminal points of neurites do not connect to more than one other terminal point.

After setting the cost matrix, the 0–1 assignment problem is defined as follows:

A = {a r g m i n}_{A} \sum_{i = 1}^{4 n} \sum_{j = 1}^{4 n} A_{i j} C_{i j}

s . t . \sum_{i = 1}^{4 n} A_{i, j} = 1 (j = 1, 2, \dots, 4 n)

\sum_{j = 1}^{4 n} A_{i, j} = 1 (i = 1, 2, \dots, 4 n)

Here, $A$ represents the connectivity matrix between different terminals of columnar regions: if $A_{i, j} = 1$ , then establish connection between terminal $i$ and terminal $j$ , if $A_{i, j} = 0$ , then establish no connection between terminal $i$ and terminal $j$ . $\sum_{i = 1}^{4 n} A_{i, j} = 1 (j = 1, 2, \dots, 4 n)$ and $\sum_{j = 1}^{4 n} A_{i, j} = 1 (i = 1, 2, \dots, 4 n)$ restrain each terminal from establishing connection with at most one other terminal. The Lapjv algorithm (Volgenant, 1996) is utilized to solve this optimization problem and the shapes of individual neurites in block images are formed. Furthermore, we employ the region growing method to generate skeletons from the reconstructed shape, achieving the neurites reconstruction from individual image blocks.

Minimal information flow tree for revising the reconstruction

The minimal information flow tree model is designed to modify the topology of skeletons, eliminate incorrect connections, and decompose them into multiple branches. When given an input skeleton file such as the swc file (Cannon et al., 1998), we convert it into a binary tree structure with the following steps.

Step 1

Request a detailed protocol

select the neurite skeleton $S_{1}$ . $S_{1}$ has the largest length in the neurite skeletons that connect with each other. One of its terminal nodes is recorded as the head node $n_{1}$ .

Step 2

Request a detailed protocol

generate the initial tree structure. Starting at head node $n_{1}$ , search the linking nodes along the skeleton $S_{1}$ , denoted by $n_{1}^{s_{1}}, n_{2}^{s_{1}}, \dots, n_{k_{1}}^{s_{1}}$ . The topology structure is $n_{i} \to l e f t n o d e = n_{i + 1}^{s_{1}}$ .

Step 3

Request a detailed protocol

generate new structure induced by the linking node $n_{1}^{s_{1}}$ . $n_{1}^{s_{1}}$ is regarded as the head node and its corresponding neurite skeleton is denoted by $S_{2}$ . Let $n_{1}^{s_{2}}, n_{2}^{s_{2}}, \dots, n_{k_{2}}^{s_{2}}$ represent the linking nodes in skeleton $S_{2}$ . The corresponding topology structure is $n_{1}^{s_{1}} \to r i g h t n o d e = n_{1}^{s_{2}}$ , $n_{i}^{s_{2}} \to l e f t n o d e = n_{i + 1}^{s_{2}}$ .

Step 4

Request a detailed protocol

repeat the operation in Step 3 for dealing with the linking nodes $n_{2}^{s_{1}}, \dots, n_{k_{1}}^{s_{1}}$ . The corresponding topology structures are added into the total tree structure. After obtaining the tree structures induced by linking nodes in $S_{1}$ , use the operation in Step 3 to generate the tree structures induced by linking nodes in $S_{2}$ . Continue in this manner until all linking nodes have been processed.

To gain a better understanding of the above process, we have provided a demonstration of how to generate the corresponding binary tree from the skeletons of neurites (Figure 1—figure supplement 4).

For the skeletons of neurites in an image block, the corresponding number of binary tree structures will be generated. We use the MIFT model to merge or split these binary structures. Suppose that an image stack contains $m$ skeletons all of which have K nodes, denoted by $n_{1}, \dots, n_{K - 1}, n_{K}$ . The connections among these nodes are stored in a matrix $W$ with $K \times K$ elements. $W_{i, j} = 0$ indicates that there is no connection between node $i$ and node $j$ . $W_{i, j} = - 1$ indicates that $j \to h e a d n o d e = i$ , $W_{i, j} = - 2$ indicates that $j \to l e f t n o d e = i$ , $W_{i, j} = - 3$ indicates that $j \to r i g h t n o d e = i$ .

The information flow can be computed as follows:

W^{*} = {a r g m i n}_{W} \sum_{i = 1}^{K} f (W, n_{i})

f (W, n_{i}) = c o s (θ (n_{i} \to h e a d n o d e, n_{i}, n_{i} \to l e f t n o d e))

Here, the optimization objective function in Equation 53 is called information flow. $θ (\cdot)$ is the angle between flow from $n_{i} \to h e a d n o d e$ to $n_{i}$ and flow from $n_{i}$ to $n_{i} \to l e f t n o d e$ . To minimize the optimization problem while ensuring that the topology matrix $W$ does not exhibit abnormal values, we adopt the strategy of dynamic programming to update the topology matrix $W$ . Briefly, we calculate the other two possible angles $θ (n_{i} \to h e a d n o d e, n_{i}, n_{i} \to r i g h t n o d e)$ and $θ (n_{i} \to l e f t n o d e, n_{i}, n_{i} \to r i g h t n o d e)$ at the first linking node $n_{i}$ . The minimum information flow is selected, and $W$ is updated. Following the updated $W$ , the next branching node is found and information flow and $W$ is updated. The updating process iterates until all nodes are updated. The final root nodes ${r_{1}, r_{2}, . . ., r_{m}}$ are obtained (node satisfies $W (r_{t}, i) = 0 o r - 1 (i = 1, . . . n)$ is set root node). The pseudo-code for solving the optimization problem is provided below:

Algorithm 3. Generation of Minimal Information Flow Tree.
# Graph defines tree topology of the nodes, t_node->left represents the left child node of t_node, t_node->right represents the right child node of t_node, t_node->head represents the head node of t_node. Input: N: ${N_{0}, N_{1}, ..., N_{k}}$ , Graph head: ${N_{0}}$ $S e t = {N_{0}}$ While $\| S e t \| > 0$ : $t_n o d e = S e t [0]$ # calculate three possible information flow $r e s = c a l c_t h r e e_d i r e c t i o n s (t_n o d e)$ if $(r e s = = 0)$ : # maintain original structure. $S e t [0] = t_n o d e - > l e f t$ $S e t . p u s h_b a c k (t_n o d e - > r i g h t)$ if $(r e s = = 1)$ : # change the position of t_node’s two child nodes. $E x c h a n g e_c h i l d (t_n o d e)$ $S e t [0] = t_n o d e - > l e f t$ $S e t . p u s h_b a c k (t_n o d e - > r i g h t)$ if $(r e s == 2)$ : # Information flows from t_node->left to t_node->right, update the structure along t_node->left and t_node->head, generate new head if possible. $N e w_n o d e = R e v e r s e_h e a d (t_n o d e)$ $S e t [0] = N e w_n o d e$ Output: N: ${N_{0}, N_{1}, . . ., N_{k}}$ , Graph head: ${N_{0}^{'}, N_{1}^{'}, . . ., N_{m}^{'}}$ .

Algorithm 3. Generation of Minimal Information Flow Tree.

# Graph defines tree topology of the nodes, t_node->left represents the left child node of t_node, t_node->right represents the right child node of t_node, t_node->head represents the head node of t_node.
Input: N:

{N_{0}, N_{1}, ..., N_{k}}

, Graph head:

{N_{0}}

S e t = {N_{0}}

While

| S e t | > 0

t_n o d e = S e t [0]

# calculate three possible information flow

r e s = c a l c_t h r e e_d i r e c t i o n s (t_n o d e)

(r e s = = 0)

:
# maintain original structure.

S e t [0] = t_n o d e - > l e f t

S e t . p u s h_b a c k (t_n o d e - > r i g h t)

(r e s = = 1)

:
# change the position of t_node’s two child nodes.

E x c h a n g e_c h i l d (t_n o d e)

S e t [0] = t_n o d e - > l e f t

S e t . p u s h_b a c k (t_n o d e - > r i g h t)

(r e s == 2)

:
# Information flows from t_node->left to t_node->right, update the structure along t_node->left and t_node->head, generate new head if possible.

N e w_n o d e = R e v e r s e_h e a d (t_n o d e)

S e t [0] = N e w_n o d e

Output: N:

{N_{0}, N_{1}, . . ., N_{k}}

, Graph head:

{N_{0}^{'}, N_{1}^{'}, . . ., N_{m}^{'}}

Please note that the model has the capability to merge binary trees. When two branches of neurites have identifiable root nodes, and one root node is in close proximity to the skeleton points on the other branch of neurites, the root node does not contribute to the calculation of information flow without fusion. However, after fusion, the root node becomes a linking node in the other branch of neurites, resulting in an additional negative information flow value. In this merging process, a threshold is required to be set. When the minimum distance between the root node of a branch of neurites and the skeleton point of the other branch of neurites is less than 8 for individual image blocks or less than 8,12,16 for fused image blocks respectively, these two branches are merged. When splitting a branch of neurites, the minimal information flow tree model is also applied to both individual and fused image blocks.

The fusion of neurites reconstruction

Request a detailed protocol

By using the MIFT model to revise the neurites reconstruction in individual image blocks, the root nodes and leaf nodes of a branch of neurites can be extracted directly. Here, we use a 0–1 assignment model to merge the reconstructions between two adjacent image blocks. For two adjacent image blocks $P$ and $Q$ , the neurite skeleton nodes which locate near the common boundary are extracted as ${p_{1}, p_{2}, ... p_{m}}$ , ${q_{1}, q_{2}, ... q_{n}}$ and the cost matrix is constructed as follows:

C = {[\begin{array}{cc} \begin{array}{ccc} c (p_{1}, q_{1}) & \dots & c (p_{1}, q_{n}) \\ ⋮ & ⋱ & ⋮ \\ c (p_{m}, q_{1}) & \dots & c (p_{m}, q_{n}) \end{array} & D_{m \times m} \\ D_{n \times n} & D_{n \times m} \end{array}]}_{(m + n) \times (m + n)}

c (p_{i}, q_{j}) = d (p_{i}, q_{j}) \times (2 - θ (L (p_{i}), L (q_{j})))

Here, $D_{m \times m}$ , $D_{n \times n}$ , $D_{n \times m}$ are auxiliary matrix which the values are all set 20. $d (p_{i}, q_{j})$ represents the Euclidean distance between terminal $p_{i}$ and $q_{j}$ . $L (p_{i})$ and $L (q_{j})$ are fitted lines from the skeleton points near $p_{i}$ and $q_{j}$ . $θ (L (p_{i}), L (q_{j}))$ represents the cosine value of their angle. Thus, the 0–1 assignment problem is formed as follows:

A = {a r g m i n}_{A} \sum_{i = 1}^{m + n} \sum_{j = 1}^{m + n} A_{i, j} \cdot C_{i, j}

s . t . \sum_{i = 1}^{m + n} A_{i, j} = 1 (j = 1, 2, \dots m + n)

\sum_{j = 1}^{m + n} A_{i, j} = 1 (i = 1, 2, \dots m + n)

Here, $A$ represents the connectivity relationship between nodes, if $A_{i, j} = 1$ , there is connection between block $P$ ’s node $i$ and block $Q$ ’s node $j$ , if $A_{i, j} = 0$ , there is no connection between block $P$ ’s node $i$ and block $Q$ ’s node $j$ . $\sum_{i = 1}^{m + n} A_{i, j} = 1 (j = 1, 2, \dots m + n)$ and $\sum_{j = 1}^{m + n} A_{i, j} = 1 (i = 1, 2, \dots m + n)$ restrict each node to connect to one other node at most. With the solved matrix $A$ , the neurite skeletons of adjacent blocks can be merged and fused skeleton structures can be obtained.

Statistical analysis

Request a detailed protocol

In this study, three commonly used metrics defined in Quan et al., 2016 were used, including precision, recall, and f1-score, which are computed to measure the fidelity between the reconstruction results and the ground truth. They are defined as follows:

p r e c i s i o n (R, G) = \frac{| R \cap G |}{| R |} = \frac{| T P |}{| R |}

r e c a l l (R, G) = \frac{| R \cap G |}{| G |} = \frac{| T P |}{| G |}

f 1 - s c o r e (R, G) = 2 \cdot \frac{p r e c i s i o n \times r e c a l l}{p r e c i s i o n + r e c a l l}

$R$ represents the point set of reconstructed neurons, $G$ represents the point set of the ground truth, $| \cdot |$ represents the number of points of a set. The three metrics are first computed on each individual neuron and then averaged by weighting each neuron with its point number of its ground truth neuritis.

We also calculated the signal-to-ratio (SNR) of the data using the following method: For a given data block $B$ and its corresponding ground-truth skeleton $S$ , we first densify the skeleton $S$ by using linear interpolation to ensure that the Euclidean distance between adjacent skeleton points is less than 1 voxel. Next, we expand each skeleton point in the densified skeleton $S^{`}$ into a spherical mask with a radius of 3 voxels. The resulting region serves as the foreground $m a s k$ . Finally, SNR is computed with mean intensity of foreground points and standard deviation of background points as follows:

M e a n_{f o r e g r o u n d} = \sum_{x \in B} I (x) \times σ_{1} (x) / \sum_{x \in B} σ_{1} (x)

M e a n_{b a c k g r o u n d} = \sum_{x \in B} I (x) \times σ_{2} (x) / \sum_{x \in B} σ_{2} (x)

S t d_{b a c k g r o u n d} = \sqrt{\sum_{x \in B} {(I (x) - M e a n_{b a c k g r o u n d})}^{2} \times σ_{2} (x) / \sum_{x \in B} σ_{2} (x)}

σ_{1} (x) = {\begin{array}{cc} 1 & i f (x \in S^{'}) \\ 0 & i f (x \notin S^{'}) \end{array}

σ_{2} (x) = {\begin{array}{cc} 1 & i f (x \notin m a s k) \\ 0 & i f (x \in m a s k) \end{array}

Here, $I (x)$ represents the signal intensity of the voxel at position $x$ , the SNR is calculated by $M e a n_{f o r e g r o u n d}$ and $S t d_{b a c k g r o u n d}$ by the following formula:

S N R = 10 l o g_{10} (M e a n_{f o r e g r o u n d} / S t d_{b a c k g r o u n d})

Data availability

The data for Figure 1C, Figure 2, Figure 3, Figure 4, Figure 5, Figure 6 and Figure 6—figure supplement 2 is available in https://zenodo.org/records/15589145. The training code of the segmentation network is available on Github: https://github.com/FateUBW0227/Seg_Net (copy archived at Cai, 2025). The software of PointTree and its user guideline are available at https://zenodo.org/records/15589145.

The following data sets were generated

1. Lin C
(2025) Zenodo
Dataset for PointTree: Automatic and accurate reconstruction of long-range axonal projections of single-neuron.

https://doi.org/10.5281/zenodo.15589145

References

1. Cai R
2. Pan C
3. Ghasemigharagoz A
4. Todorov MI
5. Förstera B
6. Zhao S
7. Bhatia HS
8. Parra-Damas A
9. Mrowka L
10. Theodorou D
11. Rempfler M
12. Xavier ALR
13. Kress BT
14. Benakis C
15. Steinke H
16. Liebscher S
17. Bechmann I
18. Liesz A
19. Menze B
20. Kerschensteiner M
21. Nedergaard M
22. Ertürk A
(2019) Panoptic imaging of transparent mice reveals whole-body neuronal projections and skull-meninges connections
Nature Neuroscience 22:317–327.

https://doi.org/10.1038/s41593-018-0301-3
- PubMed
- Google Scholar
Software
1. Cai L
(2025) Seg_Net, version swh:1:rev:5d29c54b2b41f44cc87252c265c82998ccc494c1
Software Heritage.

https://archive.softwareheritage.org/swh:1:dir:4a6b444598ac0f8f80b664fa96504d5c9a57b4fa;origin=https://github.com/FateUBW0227/Seg_Net;visit=swh:1:snp:320efc5d379135155542b22fe0f36233033c94e6;anchor=swh:1:rev:5d29c54b2b41f44cc87252c265c82998ccc494c1
(1998) An on-line archive of reconstructed hippocampal neurons
Journal of Neuroscience Methods 84:49–54.

https://doi.org/10.1016/s0165-0270(98)00091-0
- PubMed
- Google Scholar
(2012) Automatic reconstruction of neural morphologies with multi-scale tracking
Frontiers in Neural Circuits 6:25.

https://doi.org/10.3389/fncir.2012.00025
- PubMed
- Google Scholar
1. Chung K
2. Deisseroth K
(2013) CLARITY for mapping the nervous system
Nature Methods 10:508–513.

https://doi.org/10.1038/nmeth.2481
- PubMed
- Google Scholar
Conference
(2016)
3D U-Net: learning dense volumetric segmentation from sparse annotation

Medical Image Computing and Computer-Assisted Intervention–MICCAI 2016: 19th International Conference.
- Google Scholar
1. Feng L
2. Zhao T
3. Kim J
(2015) neuTube 1.0: a new design for efficient neuron reconstruction software based on the SWC Format
eNeuro 2:ENEURO.0049-14.2014.

https://doi.org/10.1523/ENEURO.0049-14.2014
- PubMed
- Google Scholar
1. Foster NN
2. Barry J
3. Korobkova L
4. Garcia L
5. Gao L
6. Becerra M
7. Sherafat Y
8. Peng B
9. Li X
10. Choi J-H
11. Gou L
12. Zingg B
13. Azam S
14. Lo D
15. Khanjani N
16. Zhang B
17. Stanis J
18. Bowman I
19. Cotter K
20. Cao C
21. Yamashita S
22. Tugangui A
23. Li A
24. Jiang T
25. Jia X
26. Feng Z
27. Aquino S
28. Mun H-S
29. Zhu M
30. Santarelli A
31. Benavidez NL
32. Song M
33. Dan G
34. Fayzullina M
35. Ustrell S
36. Boesen T
37. Johnson DL
38. Xu H
39. Bienkowski MS
40. Yang XW
41. Gong H
42. Levine MS
43. Wickersham I
44. Luo Q
45. Hahn JD
46. Lim BK
47. Zhang LI
48. Cepeda C
49. Hintiryan H
50. Dong H-W
(2021) The mouse cortico-basal ganglia-thalamic network
Nature 598:188–194.

https://doi.org/10.1038/s41586-021-03993-3
- PubMed
- Google Scholar
1. Friedmann D
2. Pun A
3. Adams EL
4. Lui JH
5. Kebschull JM
6. Grutzner SM
7. Castagnola C
8. Tessier-Lavigne M
9. Luo L
(2020) Mapping mesoscale axonal projections in the mouse brain using a 3D convolutional network
PNAS 117:11068–11075.

https://doi.org/10.1073/pnas.1918465117
- Google Scholar
1. Gao Le
2. Liu S
3. Gou L
4. Hu Y
5. Liu Y
6. Deng L
7. Ma D
8. Wang H
9. Yang Q
10. Chen Z
11. Liu D
12. Qiu S
13. Wang X
14. Wang D
15. Wang X
16. Ren B
17. Liu Q
18. Chen T
19. Shi X
20. Yao H
21. Xu C
22. Li CT
23. Sun Y
24. Li A
25. Luo Q
26. Gong H
27. Xu N
28. Yan J
(2022) Single-neuron projectome of mouse prefrontal cortex
Nature Neuroscience 25:515–529.

https://doi.org/10.1038/s41593-022-01041-5
- PubMed
- Google Scholar
1. Gao L
2. Liu S
3. Wang Y
4. Wu Q
5. Gou L
6. Yan J
(2023) Single-neuron analysis of dendrites and axons reveals the network organization in mouse prefrontal cortex
Nature Neuroscience 26:1111–1126.

https://doi.org/10.1038/s41593-023-01339-y
- PubMed
- Google Scholar
1. Harris R
(2011) Models of regional growth: past, present and future
Journal of Economic Surveys 25:913–951.

https://doi.org/10.1111/j.1467-6419.2010.00630.x
- Google Scholar
1. Helmstaedter M
(2013) Cellular-resolution connectomics: challenges of dense neural circuit reconstruction
Nature Methods 10:501–507.

https://doi.org/10.1038/nmeth.2476
- PubMed
- Google Scholar
1. Huang ZJ
2. Luo L
(2015) Neuroscience: It takes the world to understand the brain
Science 350:42–44.

https://doi.org/10.1126/science.aad4120
- PubMed
- Google Scholar
1. Huang Q
2. Chen Y
3. Liu S
4. Xu C
5. Cao T
6. Xu Y
7. Wang X
8. Rao G
9. Li A
10. Zeng S
11. Quan T
(2020) Weakly supervised learning of 3D deep network for neuron reconstruction
Frontiers in Neuroanatomy 14:38.

https://doi.org/10.3389/fnana.2020.00038
- PubMed
- Google Scholar
1. Jaqaman K
2. Loerke D
3. Mettlen M
4. Kuwata H
5. Grinstein S
6. Schmid SL
7. Danuser G
(2008) Robust single-particle tracking in live-cell time-lapse sequences
Nature Methods 5:695–702.

https://doi.org/10.1038/nmeth.1237
- PubMed
- Google Scholar
1. Kim TH
2. Schnitzer MJ
(2022) Fluorescence imaging of large-scale neural ensemble dynamics
Cell 185:9–41.

https://doi.org/10.1016/j.cell.2021.12.007
- PubMed
- Google Scholar
1. Li A
2. Gong H
3. Zhang B
4. Wang Q
5. Yan C
6. Wu J
7. Liu Q
8. Zeng S
9. Luo Q
(2010) Micro-optical sectioning tomography to obtain a high-resolution atlas of the mouse brain
Science 330:1404–1408.

https://doi.org/10.1126/science.1191776
- PubMed
- Google Scholar
1. Li R
2. Zhu M
3. Li J
4. Bienkowski MS
5. Foster NN
6. Xu H
7. Ard T
8. Bowman I
9. Zhou C
10. Veldman MB
11. Yang XW
12. Hintiryan H
13. Zhang J
14. Dong H-W
(2019) Precise segmentation of densely interweaving neuron clusters using G-Cut
Nature Communications 10:1549.

https://doi.org/10.1038/s41467-019-09515-0
- Google Scholar
1. Li S
2. Quan T
3. Zhou H
4. Huang Q
5. Guan T
6. Chen Y
7. Xu C
8. Kang H
9. Li A
10. Fu L
11. Luo Q
12. Gong H
13. Zeng S
(2020) Brain-wide shape reconstruction of a traced neuron using the convex image segmentation method
Neuroinformatics 18:199–218.

https://doi.org/10.1007/s12021-019-09434-x
- PubMed
- Google Scholar
1. Li Q
2. Shen L
(2020) 3D neuron reconstruction in tangled neuronal image with deep networks
IEEE Transactions on Medical Imaging 39:425–435.

https://doi.org/10.1109/TMI.2019.2926568
- PubMed
- Google Scholar
1. Lichtman JW
2. Denk W
(2011) The big and the small: challenges of imaging the brain’s circuits
Science 334:618–623.

https://doi.org/10.1126/science.1209168
- PubMed
- Google Scholar
1. Liu Y
2. Wang G
3. Ascoli GA
4. Zhou J
5. Liu L
(2022) Neuron tracing from light microscopy images: automation, deep learning and bench testing
Bioinformatics 38:5329–5339.

https://doi.org/10.1093/bioinformatics/btac712
- Google Scholar
1. Manubens-Gil L
2. Zhou Z
3. Chen H
4. Ramanathan A
5. Liu X
6. Liu Y
7. Bria A
8. Gillette T
9. Ruan Z
10. Yang J
11. Radojević M
12. Zhao T
13. Cheng L
14. Qu L
15. Liu S
16. Bouchard KE
17. Gu L
18. Cai W
19. Ji S
20. Roysam B
21. Wang C-W
22. Yu H
23. Sironi A
24. Iascone DM
25. Zhou J
26. Bas E
27. Conde-Sousa E
28. Aguiar P
29. Li X
30. Li Y
31. Nanda S
32. Wang Y
33. Muresan L
34. Fua P
35. Ye B
36. He H-Y
37. Staiger JF
38. Peter M
39. Cox DN
40. Simonneau M
41. Oberlaender M
42. Jefferis G
43. Ito K
44. Gonzalez-Bellido P
45. Kim J
46. Rubel E
47. Cline HT
48. Zeng H
49. Nern A
50. Chiang A-S
51. Yao J
52. Roskams J
53. Livesey R
54. Stevens J
55. Liu T
56. Dang C
57. Guo Y
58. Zhong N
59. Tourassi G
60. Hill S
61. Hawrylycz M
62. Koch C
63. Meijering E
64. Ascoli GA
65. Peng H
(2023) BigNeuron: a resource to benchmark and predict performance of algorithms for automated tracing of neurons in light microscopy datasets
Nature Methods 20:824–835.

https://doi.org/10.1038/s41592-023-01848-5
- PubMed
- Google Scholar
Book
1. McLachlan GJ
2. Krishnan T
(2007) The EM Algorithm and Extensions
John Wiley & Sons.

https://doi.org/10.1002/9780470191613
- Google Scholar
1. Meijering E
(2010) Neuron tracing in perspective
Cytometry. Part A 77:693–704.

https://doi.org/10.1002/cyto.a.20895
- PubMed
- Google Scholar
1. Muñoz-Castañeda R
2. Zingg B
3. Matho KS
4. Chen X
5. Wang Q
6. Foster NN
7. Li A
8. Narasimhan A
9. Hirokawa KE
10. Huo B
11. Bannerjee S
12. Korobkova L
13. Park CS
14. Park Y-G
15. Bienkowski MS
16. Chon U
17. Wheeler DW
18. Li X
19. Wang Y
20. Naeemi M
21. Xie P
22. Liu L
23. Kelly K
24. An X
25. Attili SM
26. Bowman I
27. Bludova A
28. Cetin A
29. Ding L
30. Drewes R
31. D’Orazi F
32. Elowsky C
33. Fischer S
34. Galbavy W
35. Gao L
36. Gillis J
37. Groblewski PA
38. Gou L
39. Hahn JD
40. Hatfield JT
41. Hintiryan H
42. Huang JJ
43. Kondo H
44. Kuang X
45. Lesnar P
46. Li X
47. Li Y
48. Lin M
49. Lo D
50. Mizrachi J
51. Mok S
52. Nicovich PR
53. Palaniswamy R
54. Palmer J
55. Qi X
56. Shen E
57. Sun Y-C
58. Tao HW
59. Wakemen W
60. Wang Y
61. Yao S
62. Yuan J
63. Zhan H
64. Zhu M
65. Ng L
66. Zhang LI
67. Lim BK
68. Hawrylycz M
69. Gong H
70. Gee JC
71. Kim Y
72. Chung K
73. Yang XW
74. Peng H
75. Luo Q
76. Mitra PP
77. Zador AM
78. Zeng H
79. Ascoli GA
80. Josh Huang Z
81. Osten P
82. Harris JA
83. Dong H-W
(2021) Cellular anatomy of the mouse primary motor cortex
Nature 598:159–166.

https://doi.org/10.1038/s41586-021-03970-w
- PubMed
- Google Scholar
1. Osten P
2. Margrie TW
(2013) Mapping brain circuitry with a light microscope
Nature Methods 10:515–523.

https://doi.org/10.1038/nmeth.2477
- PubMed
- Google Scholar
1. Parekh R
2. Ascoli GA
(2013) Neuronal morphology goes digital: a research hub for cellular and system neuroscience
Neuron 77:1017–1038.

https://doi.org/10.1016/j.neuron.2013.03.008
- PubMed
- Google Scholar
1. Peng H
2. Long F
3. Myers G
(2011) Automatic 3D neuron tracing using all-path pruning
Bioinformatics 27:i239–i247.

https://doi.org/10.1093/bioinformatics/btr237
- PubMed
- Google Scholar
1. Peng H
2. Hawrylycz M
3. Roskams J
4. Hill S
5. Spruston N
6. Meijering E
7. Ascoli GA
(2015) BigNeuron: large-scale 3D neuron reconstruction from optical microscopy images
Neuron 87:252–256.

https://doi.org/10.1016/j.neuron.2015.06.036
- PubMed
- Google Scholar
1. Peng H
2. Xie P
3. Liu L
4. Kuang X
5. Wang Y
6. Qu L
7. Gong H
8. Jiang S
9. Li A
10. Ruan Z
11. Ding L
12. Yao Z
13. Chen C
14. Chen M
15. Daigle TL
16. Dalley R
17. Ding Z
18. Duan Y
19. Feiner A
20. He P
21. Hill C
22. Hirokawa KE
23. Hong G
24. Huang L
25. Kebede S
26. Kuo H-C
27. Larsen R
28. Lesnar P
29. Li L
30. Li Q
31. Li X
32. Li Y
33. Li Y
34. Liu A
35. Lu D
36. Mok S
37. Ng L
38. Nguyen TN
39. Ouyang Q
40. Pan J
41. Shen E
42. Song Y
43. Sunkin SM
44. Tasic B
45. Veldman MB
46. Wakeman W
47. Wan W
48. Wang P
49. Wang Q
50. Wang T
51. Wang Y
52. Xiong F
53. Xiong W
54. Xu W
55. Ye M
56. Yin L
57. Yu Y
58. Yuan J
59. Yuan J
60. Yun Z
61. Zeng S
62. Zhang S
63. Zhao S
64. Zhao Z
65. Zhou Z
66. Huang ZJ
67. Esposito L
68. Hawrylycz MJ
69. Sorensen SA
70. Yang XW
71. Zheng Y
72. Gu Z
73. Xie W
74. Koch C
75. Luo Q
76. Harris JA
77. Wang Y
78. Zeng H
(2021) Morphological diversity of single neurons in molecularly defined cell types
Nature 598:174–181.

https://doi.org/10.1038/s41586-021-03941-1
- PubMed
- Google Scholar
1. Qiu S
2. Hu Y
3. Huang Y
4. Gao T
5. Wang X
6. Wang D
7. Ren B
8. Shi X
9. Chen Y
10. Wang X
11. Wang D
12. Han L
13. Liang Y
14. Liu D
15. Liu Q
16. Deng L
17. Chen Z
18. Zhan L
19. Chen T
20. Huang Y
21. Wu Q
22. Xie T
23. Qian L
24. Jin C
25. Huang J
26. Deng W
27. Jiang T
28. Li X
29. Jia X
30. Yuan J
31. Li A
32. Yan J
33. Xu N
34. Xu L
35. Luo Q
36. Poo M-M
37. Sun Y
38. Li CT
39. Yao H
40. Gong H
41. Sun Y-G
42. Xu C
(2024) Whole-brain spatial organization of hippocampal single-neuron projectomes
Science 383:eadj9198.

https://doi.org/10.1126/science.adj9198
- PubMed
- Google Scholar
1. Quan T
2. Zhou H
3. Li J
4. Li S
5. Li A
6. Li Y
7. Lv X
8. Luo Q
9. Gong H
10. Zeng S
(2016) NeuroGPS-Tree: automatic reconstruction of large-scale neuronal populations with dense neurites
Nature Methods 13:51–54.

https://doi.org/10.1038/nmeth.3662
- PubMed
- Google Scholar
1. Radojevic M
2. Meijering E
(2017) Automated neuron tracing using probability hypothesis density filtering
Bioinformatics 33:1073–1080.

https://doi.org/10.1093/bioinformatics/btw751
- PubMed
- Google Scholar
Book
1. Reynolds DA
(2009) Gaussian mixture models
In: Li SZ, Jain A, editors. Encyclopedia of Biometrics. Springer. pp. 659–663.

https://doi.org/10.1007/978-0-387-73003-5_196
- Google Scholar
1. Stelzer EHK
2. Strobl F
3. Chang B-J
4. Preusser F
5. Preibisch S
6. McDole K
7. Fiolka R
(2021) Light sheet fluorescence microscopy
Nature Reviews Methods Primers 1:73.

https://doi.org/10.1038/s43586-021-00069-4
- Google Scholar
1. Sun P
2. Freund RM
(2004) Computation of minimum-volume covering ellipsoids
Operations Research 52:690–706.

https://doi.org/10.1287/opre.1040.0115
- Google Scholar
1. Sun Q
2. Li X
3. Ren M
4. Zhao M
5. Zhong Q
6. Ren Y
7. Luo P
8. Ni H
9. Zhang X
10. Zhang C
11. Yuan J
12. Li A
13. Luo M
14. Gong H
15. Luo Q
(2019) A whole-brain map of long-range inputs to GABAergic interneurons in the mouse medial prefrontal cortex
Nature Neuroscience 22:1357–1370.

https://doi.org/10.1038/s41593-019-0429-9
- PubMed
- Google Scholar
1. Türetken E
2. González G
3. Blum C
4. Fua P
(2011) Automated reconstruction of dendritic and axonal trees by global optimization with geometric priors
Neuroinformatics 9:279–302.

https://doi.org/10.1007/s12021-011-9122-1
- PubMed
- Google Scholar
1. Volgenant A
(1996) Linear and semi-assignment problems: a core oriented approach
Computers & Operations Research 23:917–932.

https://doi.org/10.1016/0305-0548(96)00010-X
- Google Scholar
1. Wang Y
2. Li Q
3. Liu L
4. Zhou Z
5. Ruan Z
6. Kong L
7. Li Y
8. Wang Y
9. Zhong N
10. Chai R
11. Luo X
12. Guo Y
13. Hawrylycz M
14. Luo Q
15. Gu Z
16. Xie W
17. Zeng H
18. Peng H
(2019) TeraVR empowers precise reconstruction of complete 3-D neuronal morphology in the whole brain
Nature Communications 10:3474.

https://doi.org/10.1038/s41467-019-11443-y
- Google Scholar
1. Wang X
2. Xiong H
3. Liu Y
4. Yang T
5. Li A
6. Huang F
7. Yin F
8. Su L
9. Liu L
10. Li N
11. Li L
12. Cheng S
13. Liu X
14. Lv X
15. Liu X
16. Chu J
17. Xu T
18. Xu F
19. Gong H
20. Luo Q
21. Yuan J
22. Zeng S
(2021) Chemical sectioning fluorescence tomography: high-throughput, high-contrast, multicolor, whole-brain imaging at subcellular resolution
Cell Reports 34:108709.

https://doi.org/10.1016/j.celrep.2021.108709
- Google Scholar
1. Wei X
2. Peng M
3. Huang H
4. Zhou Y
(2023) An overview on density peaks clustering
Neurocomputing 554:126633.

https://doi.org/10.1016/j.neucom.2023.126633
- Google Scholar
1. Winnubst J
2. Bas E
3. Ferreira TA
4. Wu Z
5. Economo MN
6. Edson P
7. Arthur BJ
8. Bruns C
9. Rokicki K
10. Schauder D
11. Olbris DJ
12. Murphy SD
13. Ackerman DG
14. Arshadi C
15. Baldwin P
16. Blake R
17. Elsayed A
18. Hasan M
19. Ramirez D
20. Dos Santos B
21. Weldon M
22. Zafar A
23. Dudman JT
24. Gerfen CR
25. Hantman AW
26. Korff W
27. Sternson SM
28. Spruston N
29. Svoboda K
30. Chandrashekar J
(2019) Reconstruction of 1,000 projection neurons reveals new cell types and organization of long-range connectivity in the mouse brain
Cell 179:268–281.

https://doi.org/10.1016/j.cell.2019.07.042
- PubMed
- Google Scholar
1. Wu J
2. He Y
3. Yang Z
4. Guo C
5. Luo Q
6. Zhou W
7. Chen S
8. Li A
9. Xiong B
10. Jiang T
11. Gong H
(2014) 3D BrainCV: simultaneous visualization and analysis of cells and capillaries in a whole mouse brain with one-micron voxel resolution
NeuroImage 87:199–208.

https://doi.org/10.1016/j.neuroimage.2013.10.036
- PubMed
- Google Scholar
1. Xiao H
2. Peng H
(2013) APP2: automatic tracing of 3D neuron morphology based on hierarchical pruning of a gray-weighted image distance-tree
Bioinformatics 29:1448–1454.

https://doi.org/10.1093/bioinformatics/btt170
- PubMed
- Google Scholar
1. Xu F
2. Shen Y
3. Ding L
4. Yang C-Y
5. Tan H
6. Wang H
7. Zhu Q
8. Xu R
9. Wu F
10. Xiao Y
11. Xu C
12. Li Q
13. Su P
14. Zhang LI
15. Dong H-W
16. Desimone R
17. Xu F
18. Hu X
19. Lau P-M
20. Bi G-Q
(2021) High-throughput mapping of a whole rhesus monkey brain at micrometer resolution
Nature Biotechnology 39:1521–1528.

https://doi.org/10.1038/s41587-021-00986-5
- PubMed
- Google Scholar
(2013) A distance-field based automatic neuron tracing method
BMC Bioinformatics 14:1–11.

https://doi.org/10.1186/1471-2105-14-93
- PubMed
- Google Scholar
1. Zeng H
2. Sanes JR
(2017) Neuronal cell-type classification: challenges, opportunities and the path forward
Nature Reviews. Neuroscience 18:530–546.

https://doi.org/10.1038/nrn.2017.85
- PubMed
- Google Scholar
1. Zeng H
(2022) What is a cell type and how to define it?
Cell 185:2739–2755.

https://doi.org/10.1016/j.cell.2022.06.031
- PubMed
- Google Scholar
1. Zhang H
2. Wang X
3. Guo W
4. Li A
5. Chen R
6. Huang F
7. Liu X
8. Chen Y
9. Li N
10. Liu X
11. Xu T
12. Xue Z
13. Zeng S
(2021) Cross-streams through the ventral posteromedial thalamic nucleus to convey vibrissal information
Frontiers in Neuroanatomy 15:724861.

https://doi.org/10.3389/fnana.2021.724861
- PubMed
- Google Scholar
1. Zhong Q
2. Li A
3. Jin R
4. Zhang D
5. Li X
6. Jia X
7. Ding Z
8. Luo P
9. Zhou C
10. Jiang C
11. Feng Z
12. Zhang Z
13. Gong H
14. Yuan J
15. Luo Q
(2021) High-definition imaging using line-illumination modulation microscopy
Nature Methods 18:309–315.

https://doi.org/10.1038/s41592-021-01074-x
- PubMed
- Google Scholar
1. Zhou Z
2. Kuo HC
3. Peng H
4. Long F
(2018) DeepNeuron: an open deep learning toolbox for neuron tracing
Brain Informatics 5:3.

https://doi.org/10.1186/s40708-018-0081-2
- PubMed
- Google Scholar
1. Zhou H
2. Li S
3. Li A
4. Huang Q
5. Xiong F
6. Li N
7. Han J
8. Kang H
9. Chen Y
10. Li Y
11. Lin H
12. Zhang Y-H
13. Lv X
14. Liu X
15. Gong H
16. Luo Q
17. Zeng S
18. Quan T
(2021) GTree: an open-source tool for dense reconstruction of brain-wide neuronal population
Neuroinformatics 19:305–317.

https://doi.org/10.1007/s12021-020-09484-6
- PubMed
- Google Scholar
1. Zingg B
2. Hintiryan H
3. Gou L
4. Song MY
5. Bay M
6. Bienkowski MS
7. Foster NN
8. Yamashita S
9. Bowman I
10. Toga AW
11. Dong H-W
(2014) Neural networks of the mouse neocortex
Cell 156:1096–1111.

https://doi.org/10.1016/j.cell.2014.02.023
- PubMed
- Google Scholar

Article and author information

Author details

Lin Cai
1. Britton Chance Center for Biomedical Photonics, Wuhan National Laboratory for Optoelectronics, Huazhong University of Science and Technology, Wuhan, China
2. MOE Key Laboratory for Biomedical Photonics, Wuhan National Laboratory for Optoelectronics, Huazhong University of Science and Technology, Wuhan, China
Contribution
Validation, Investigation, Methodology, Writing – original draft, Writing – review and editing

Competing interests
No competing interests declared

"This ORCID iD identifies the author of this article:" 0000-0002-4413-3599
Taiyu Fan
1. Britton Chance Center for Biomedical Photonics, Wuhan National Laboratory for Optoelectronics, Huazhong University of Science and Technology, Wuhan, China
2. MOE Key Laboratory for Biomedical Photonics, Wuhan National Laboratory for Optoelectronics, Huazhong University of Science and Technology, Wuhan, China
Contribution
Validation, Investigation, Writing – original draft

Competing interests
No competing interests declared
Xuzhong Qu
1. Britton Chance Center for Biomedical Photonics, Wuhan National Laboratory for Optoelectronics, Huazhong University of Science and Technology, Wuhan, China
2. MOE Key Laboratory for Biomedical Photonics, Wuhan National Laboratory for Optoelectronics, Huazhong University of Science and Technology, Wuhan, China
Contribution
Software, Investigation

Competing interests
No competing interests declared
Ying Zhang
1. Britton Chance Center for Biomedical Photonics, Wuhan National Laboratory for Optoelectronics, Huazhong University of Science and Technology, Wuhan, China
2. MOE Key Laboratory for Biomedical Photonics, Wuhan National Laboratory for Optoelectronics, Huazhong University of Science and Technology, Wuhan, China
Contribution
Validation

Competing interests
No competing interests declared
Xianyu Gou
1. Britton Chance Center for Biomedical Photonics, Wuhan National Laboratory for Optoelectronics, Huazhong University of Science and Technology, Wuhan, China
2. MOE Key Laboratory for Biomedical Photonics, Wuhan National Laboratory for Optoelectronics, Huazhong University of Science and Technology, Wuhan, China
Contribution
Investigation

Competing interests
No competing interests declared
Quanwei Ding
1. Britton Chance Center for Biomedical Photonics, Wuhan National Laboratory for Optoelectronics, Huazhong University of Science and Technology, Wuhan, China
2. MOE Key Laboratory for Biomedical Photonics, Wuhan National Laboratory for Optoelectronics, Huazhong University of Science and Technology, Wuhan, China
3. School of Computer Science and Engineering, Hubei Key Laboratory of Intelligent Robot, Wuhan Institute of Technology, Wuhan, China
Contribution
Validation

Competing interests
No competing interests declared
Weihua Feng
1. Britton Chance Center for Biomedical Photonics, Wuhan National Laboratory for Optoelectronics, Huazhong University of Science and Technology, Wuhan, China
2. MOE Key Laboratory for Biomedical Photonics, Wuhan National Laboratory for Optoelectronics, Huazhong University of Science and Technology, Wuhan, China
Contribution
Investigation

Competing interests
No competing interests declared
Tingting Cao
1. Britton Chance Center for Biomedical Photonics, Wuhan National Laboratory for Optoelectronics, Huazhong University of Science and Technology, Wuhan, China
2. MOE Key Laboratory for Biomedical Photonics, Wuhan National Laboratory for Optoelectronics, Huazhong University of Science and Technology, Wuhan, China
Contribution
Validation

Competing interests
No competing interests declared
Xiaohua Lv
1. Britton Chance Center for Biomedical Photonics, Wuhan National Laboratory for Optoelectronics, Huazhong University of Science and Technology, Wuhan, China
2. MOE Key Laboratory for Biomedical Photonics, Wuhan National Laboratory for Optoelectronics, Huazhong University of Science and Technology, Wuhan, China
Contribution
Supervision

Competing interests
No competing interests declared
Xiuli Liu
1. Britton Chance Center for Biomedical Photonics, Wuhan National Laboratory for Optoelectronics, Huazhong University of Science and Technology, Wuhan, China
2. MOE Key Laboratory for Biomedical Photonics, Wuhan National Laboratory for Optoelectronics, Huazhong University of Science and Technology, Wuhan, China
Contribution
Supervision

Competing interests
No competing interests declared
Qing Huang
1. Britton Chance Center for Biomedical Photonics, Wuhan National Laboratory for Optoelectronics, Huazhong University of Science and Technology, Wuhan, China
2. MOE Key Laboratory for Biomedical Photonics, Wuhan National Laboratory for Optoelectronics, Huazhong University of Science and Technology, Wuhan, China
3. School of Computer Science and Engineering, Hubei Key Laboratory of Intelligent Robot, Wuhan Institute of Technology, Wuhan, China
Contribution
Supervision

Competing interests
No competing interests declared
Tingwei Quan
1. Britton Chance Center for Biomedical Photonics, Wuhan National Laboratory for Optoelectronics, Huazhong University of Science and Technology, Wuhan, China
2. MOE Key Laboratory for Biomedical Photonics, Wuhan National Laboratory for Optoelectronics, Huazhong University of Science and Technology, Wuhan, China
Contribution
Supervision, Funding acquisition, Investigation, Methodology, Writing – original draft, Writing – review and editing

For correspondence
quantingwei@hust.edu.cn

Competing interests
No competing interests declared

"This ORCID iD identifies the author of this article:" 0000-0002-8393-4292
Shaoqun Zeng
1. Britton Chance Center for Biomedical Photonics, Wuhan National Laboratory for Optoelectronics, Huazhong University of Science and Technology, Wuhan, China
2. MOE Key Laboratory for Biomedical Photonics, Wuhan National Laboratory for Optoelectronics, Huazhong University of Science and Technology, Wuhan, China
Contribution
Supervision

Competing interests
No competing interests declared

Funding

National Natural Science Foundation of China (32471146)

Tingwei Quan

PLA General Hospital (N20240194)

Tingwei Quan

The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.

Acknowledgements

We thank the members of the Britton Chance Center for Biomedical Photonics for advice and help in experiments. This work was supported by the National Natural Science Foundation of China (32471146) and the project N20240194.

Ethics

All animal experiments followed procedures approved by the Institutional Animal Ethics Committee of the Huazhong University of Science and Technology.

Version history

Sent for peer review: September 20, 2024
Preprint posted: September 23, 2024
Reviewed Preprint version 1: November 11, 2024
Reviewed Preprint version 2: July 22, 2025
Version of Record published: September 1, 2025

Cite all versions

You can cite all versions using the DOI https://doi.org/10.7554/eLife.102840. This DOI represents all versions, and will always resolve to the latest one.

Copyright

This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.