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^1-4. Advances in optical imaging and molecular labeling techniques^5-10 have allowed us to observe the entire mouse brain at single-axon resolution, and provided the database for the study of neuronal projection patterns^11-18. However, the reconstruction of these long-range projected axons still requires frequent manual manipulation in tens of TBs volumetric images^7,19-22, which severely hinders high-throughput acquisition of neuronal projection patterns²³.

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 are still densely distributed due to the morphological complexity of neurons. The identification of densely distributed axons is considered an open problem in the field^23-25, 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 to other neurons or missing²⁶. 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 framework focuses on how to accurately extract skeletons and establish the connections between skeletons^2,27. The BigNeuron project²⁸ 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^29-32, the localization of the next skeleton point position requires one to compute the signal anisotropy of the image region near the current skeleton point. Localization errors typically occur when this image region includes other neurite signals. The global approach^24,33,34 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. We also note that deep learning is widely used in neuron reconstruction^35-38, mostly to segment image regions of neurons and enhance the signal intensity in the segmented regions, which helps in reducing the reconstruction error. However, in the ideal case where all neurite centers are identified and their signals are enhanced, these reconstruction methods based on skeleton extraction still suffer from many errors (Fig.S1).

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^24,39 achieves 80% reconstruction accuracy from GB-scale images under certain constrains. First, for most neurites, their terminals and junctions require to be accurately identified, which is used to calculate the angle between two neurites and to make neuronal morphological statistics available. Second, cell bodies 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. Current long-range projection reconstruction methods are semi-automatic and require substantial human intervention^19,21,22,40.

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, a single columnar region is characterized by a minimal envelope ellipsoid for constructing connections between columnar regions, which forms the neurite shape. Based on this reconstruction, 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⁴¹ to partition the foreground points into columnar regions and determine their geometrical parameters by solving the minimum-volume covering ellipsoids problem⁴². Using these geometrical parameters, we construct a 0-1 assignment problem⁴³ to establish links between these columnar regions. Finally, skeletons are extracted from these linked columnar regions to reduce data redundancy by using region growing⁴⁴. The key procedures for neuron reconstruction are presented in Figure 1A.

Figure 1.

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 smaller angle, called the node angle, corresponds to the segment that serves as part of a single neurite (Figure 1B). In this way, 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 minimal information flow tree (MIFT).

In image blocks of densely distributed neurites, we used semi-automatic software²¹ extracting 208 neuronal branches and identifying their root nodes, and calculated their information flow values, i.e., the ground-truth values. We looped through all possible structure of these branches by changing the root node in order to compute the maximum and minimum information flow scores (Figure 1C). It is evident that, for most neuronal branches, 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 Fig.S2), 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 which behaves as crossover and closely paralleled neurites. These neurites can be manually identified by visualization with different view angles (Fig.S3). We compared PointTree with several skeleton-based methods such as neuTube⁴⁵, PHDF⁴⁶, NGPST³⁹ and MOST⁴⁷ in performing this task. We manually labeled the locations where neurites are crossover or closely parallel from five 256×256×256 image blocks. For fair comparison, all methods are performed on segmented images derived from the segmentation network. Figure 2A illustrates the process of 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.

Furthermore, we present the quantitative results derived from PointTree and five widely used skeleton-based reconstruction methods. including APP2, neuTube, NGPST, PHDF, 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 dataset. 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 the good solution to separate the densely neurites. The F1-score of these reconstructions ranges from 30% to 40%, which indicates the ineffective reconstructions.

Figure 3.

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 4A 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 4A). 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 4A).

Figure 4.

We furthermore measure the enhancement in the reconstruction accuracy achieved by MIFT (Figure 4B). 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 4A), 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, F1 score is less than 0.2. We also present some reconstruction examples after two fusions in Figure 4C, which are close to the ground truth. These results suggest that MIFT model take 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 5A). We also used GTree to manually reconstruct these neurons as the ground-truth reconstruction (Figure 5B). 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 5C and Movies S1 & S2), and we compared these reconstructions with ground truth (Fig.S4). The average precision is above 85% and the average recall and F1-score are above 80% (Figure 5E). In addition, we presented the axons reconstructions from two image blocks (Figures 5C₁&C₂) 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 5.

We also applied PointTree to process another 10739×11226×3921 image blocks collected with HD-fMOST system⁴⁸. 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 (Fig.S5). The quantitative results suggest that the average F1-score is above 90% (Table S1).

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 minutes. 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 minutes (Table S2). 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 (Movie S3).

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 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 to 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 well 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 derived from the branching morphological statistics may cause inconsistencies in rare cases when correcting for reconstruction errors. 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

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⁴⁹ was used to acquire imaging data (Figures 1-5), more details can be seen in the reference ⁵⁰. For one C57BL/6J 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 (Fig.S5), more details can be seen in the reference⁴⁸.

Generation of foreground points

Our reconstruction method performs on the image foregrounds. Here, we used UNet3D⁵¹ for image stacks segmentation without network structure modification. The detailed information about UNet3D can be found in the reference⁵¹. 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 supervise the network. Initially, the semi-automatic software GTree was utilized to extract the neurites 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 um. 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 forms 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 loss function is used to supervise the network outputs to the input labels. We presented the reconstructions from two kinds of fMOST datasets. One is from the reference ⁵⁰ and the other is from the reference ⁴⁸. Therefore, we created two sets of training data, each consisting of 20 512×512×512 image blocks. 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

For the image stack, we allocated the foreground points to their respective neurites and established connections between neurites by constructing three optimization models. The constrained Gaussian mixture model divides the foreground points into a set of points, each of which has a column shape. The minimum-volume covering ellipsoids model extracts the features of the column-shaped point set. 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

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 and obtain the probability density function P(x) with respect to the foreground point x. For each foreground point x, its probabilities in K Gaussian distributions are calculated. The foreground points are assigned to the k-th clusters according to their maximum probability and form k-th columnar regions. Considering that both the number of foreground points and K are large, we have added some constrained conditions for Gaussian mixture model as follows.

N (x|μ_k, Σ_k) is the Gaussian density function with mean value μ_k and covariance matrix Σ_k. refers to the fact that the total probability distribution normalizes to 1. I (·) represents the signal intensity from segment image, ϵ₀ is the minimum signal intensity of foreground points. I (μ_i) ≥ ε₀ restrain the center of the gaussian distribution to be a foreground point. | Σ_i |≤ ε₁ restrain the determinant of the covariance matrix which control the suitable number of foreground points for each columnar region. ϵ₁ is set to cube of three times the average diameter of neurite.

Given the foreground points x₁, x₂, …, x_N, we employ maximum likelihood to estate the parameters of Gaussian mixture model.

In the solving this optimization problem, we employ peak density algorithm to compute density for each foreground points and sort them in descending order. We select signal points from the median to both sides as seed points which can decrease the situations that seed points lies 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 (μ₁, μ₂, …, μ_K). The initial setting of covariance matrix is the identity matrix. The constrained Gaussian mixture model was solved by EM algorithm⁵² (Supplementary Text).

Shape characterization of columnar regions

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 use the ellipsoid that can fully encompasses a single columnar region and meantime has the minimum volume. For, a three-dimensional ellipsoid can be defined as follow:

Here, c is the center of ellipsoid, Q represents the geometric shape, denotes the convex cone of 3×3 symmetric positive definite matrices. The volume of E_c,Q is given by the formula:

Here, Г(·) is the standard gamma function of calculus.

For a columnar region with foreground points P{x₁, x₂, … x_m }, we define the target function as follow:

Here c ∈ CHull(P_i) restrain the solved center of ellipsoid to locate within the smallest convex hull formed by the clustering points. Through solving the above minimum-volume covering ellipsoids problem (Supplementary Text), we can characterize the columnar regions more accurate.

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 proportionally enlarged, but the volume of the corresponding ellipsoid is not minimum. These two cases will reduce the accuracy of the connections between clusters, i.e., 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

The 0-1 assignment model⁴³ can robustly and accurately establish connections between particles in live-cell imaging⁵³. 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 equals to 0 or 1. According to the direction and the endpoints, we design the cost matrix for building the 0-1 assignment model.

Here D is 2n×2n auxiliary matrix all element of which are all set 100. Both i0 and j0 in EQ (11) are equal to 0 or 1, labeling the two endpoints of the longest axis of the ellipsoid. norm(t_{i,i 0}, t _{j, j 0}) represents the Euclidean distance between t_{i,i 0} and t _{j, j 0}. θ(t_i,i0, t _{j, j0}) describes the angle between two ellipsoids, i.e., two columnar regions. dir(c_i, t_{i,i 0}) represents the line from point c_i to t_{i,i 0}. ⟨dir(c_i, t_{i,i 0}), dir(c_i, t _{j, j 0})⟩ represents cosine angle between the two lines. The threshold of 100 in D in EQ (10) and EQ (11) is an experimental value designed to ensure that the terminal points of neurites do not connect to more than one other terminal points.

After set the cost matrix, the 0-1 assignment problem is defined as follow:

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. and restrain each terminal establish connection with at most one other terminal. The Lapjv algorithm⁴³ 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⁵⁴, we convert it into a binary tree structure with following steps.

Step1: select the neurite skeleton S_1. S₁ has the largest length in the neurite skeletons that connect with each other. One of its terminal nodes are recorded as the head node n₁.

Step2: generate the initial tree structure. Starting at head node n₁, search the linking nodes along the skeleton S₁, denoted by . The topology structure is .

Step3: generate new structure induced by the linking node is regarded as the head node and its corresponding neurite skeleton is denoted by S₂. Let represent the linking nodes in skeleton S₂The corresponding topology structure is .

Step4: repeat the operation in Step3 for dealing with the linking nodes . The corresponding topology structures are added into the total tree structure. After obtaining the tree structures induced by linking nodes in S₁, use the operation in Step 3 to generate the tree structures induced by linking nodes in S₂. 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 (Fig.S6).

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₁, …, n_{K -1}, n_K. The connections among these nodes are stored in a matrix W with K ’ s 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 → headnode = i, W_{i, j} = -2 indicates that j → leftnode = i, W_{i, j} = -3 indicates that j → rightnode = i.

The information flow can be computed as follow:

Here the optimization objective function in EQ (17) is called information flow. θ(·) is the angle between flow from n_i → headnode to n_i and flow from n_i to n_i → leftnode. To minimize optimization problem while ensuring that the topology matrix W does not exhibit abnormal values, we adopt the strategy of dynamic programming to update topology matrix W. Briefly, we calculate the other two possible angles θ(n_i → headnode, n_i, n_i → rightnode) and θ(n_i → leftnode, n_i, n_i → rightnode) 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₁, r₂, …r_m } are obtained (node satisfies W (r_t, i) = 0 or -1(i = 1,…n) is set root node). The pseudo-code for solving the optimization problem can be found in Supplementary Text.

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

By using 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 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 {p, p, … p },{q₁, q₂, …q } and the cost matrix is constructed as follow:

Here D_m×m, D_n×n, D_n×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 follow:

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. and restrict each node connect to one other node at most. With solved matrix A, the neurite skeletons of adjacent blocks can be merged and fused skeleton structures can be obtained.

Statistical Analysis

In this study, three commonly used metrics defined in ³⁹ were used, including precision, recall and F1-score are computed to measure the fidelity between the reconstruction results and the ground truth. They are defined as follow:

R represents the point set of reconstructed neurons, G represents the point set of the ground truth, | · | 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.

Acknowledgements

We thank the members of the Britton Chance Center for Biomedical Photonics for advice and help in experiments.

Additional information

Funding

Advanced Biomedical Imaging facilities (2105-420118-89-01-121873).

Author contributions

Methodology: T.Q., L.C.

Investigation: T.Q., L.C., T.F., X.Q., X.G., W.F.

Validation: L.C., T.F., Q.D, T.C, Y.Z.

Supervision: S.Z., T.Q., X.L., X.L., Q.H.

Writing: T.Q., L.C., T.F.

Competing interests

Authors declare that they have no competing interests.

Data and materials availability

The data for Figure 1C, Figure 2, Figure 3, Figure 4, Figure 5 and Fig.S5 is available in https://zenodo.org/records/11239828. The raw image blocks are extremely large and can be available on request from the corresponding author. The training code of the segmentation network is available on Github: https://github.com/FateUBW0227/Seg_Net. The software of PointTree and its user guideline are available on https://zenodo.org/records/11239828.