|Home | About | Journals | Submit | Contact Us | Français|
Motivation: Unraveling the structure and behavior of the brain and central nervous system (CNS) has always been a major goal of neuroscience. Understanding the wiring diagrams of the neuromuscular junction connectomes (full connectivity of nervous system neuronal components) is a starting point for this, as it helps in the study of the organizational and developmental properties of the mammalian CNS. The phenomenon of synapse elimination during developmental stages of the neuronal circuitry is such an example. Due to the organizational specificity of the axons in the connectomes, it becomes important to label and extract individual axons for morphological analysis. Features such as axonal trajectories, their branching patterns, geometric information, the spatial relations of groups of axons, etc. are of great interests for neurobiologists in the study of wiring diagrams. However, due to the complexity of spatial structure of the axons, automatically tracking and reconstructing them from microscopy images in 3D is an unresolved problem. In this article, AxonTracker-3D, an interactive 3D axon tracking and labeling tool is built to obtain quantitative information by reconstruction of the axonal structures in the entire innervation field. The ease of use along with accuracy of results makes AxonTracker-3D an attractive tool to obtain valuable quantitative information from axon datasets.
Availability: The software is freely available for download at http://www.cbi-tmhs.org/AxonTracker/
The modern microscopic image acquisition techniques provide the neuroscientists the visual perception of axons in 3D space. Of particular interest is the knowledge of connectivity of the neuronal components of the nervous system, the connectome. Analyzing neuronal structures is instrumental in completely understanding the functional properties of the mammalian central nervous system (White et al., 1986). One such phenomenon that can be studied is the developmental reorganization of connectivity in the mammalian nervous system, also known as synapse elimination (Purves and Lichtman, 1980). Pruning of synaptic connections through competitive interactions during neurodevelopment and its effect on the neuronal connectivity is one such example.
A study of connectomes over various time points during brain development can reveal the mysteries about how the neural circuits are established and the way it changes over the lifespan of the animal under study. Moreover, mapping wiring diagram mouse models can help understand how an abnormality in the neural circuits can lead to psychiatric disorders, such as schizophrenia and autism (Lichtman et al., 2008). Another potential application is in the study of the effect of new drugs on the mammalian connectomes. Analysis of the connectomes at different time points during the drug-testing could lead to valuable insights into the effectiveness of the new drugs designed to treat disorders related to brain and central nervous system.
The wiring diagram of the neuromuscular circuit is believed to be a good starting point for this, as it might lead to discoveries valuable to this study. The neuromuscular circuits are well-isolated and accessible. Its connectome is well defined. Due to the large size of the adult motor neurons (2–10 μm in diameter), which are separated from neighboring axons by Schwann cells, and the large size of neuromuscular junctions (~20 μm in length), the neuromuscular connectome is easier to resolve, analyze and can be easily visualized. Moreover, since the morphological features of the neuromuscular connectome can be directly linked with its functional properties; for example the morphology of individual neuromuscular junction has been found to be directly related to its synaptic strength (Herrera et al., 1985), it renders itself as an attractive platform to conduct connectome studies. Many questions related to synapse elimination can be answered from the study of mammalian neuromuscular connectome. As the spatial organization of the axons is very specific in nature, down to the level of individual axons, it becomes important to isolate each axon from the connectomes for extraction of morphological features for these studies.
Studies in neuronal connections are believed to help unravel the functional properties of the nervous system. The manual effort needed, however, is excessive, as the neurobiologist has to manually track individual axons from a large set of data (~30–50 GB of volumetric data). This might amount to a few months of painstaking labor. Of note are similar researches in the mapping of neural circuits for the study of neuronal connectivity. A study of potential connectivity among neurons was done based on neuron reconstructions of the cortical circuits in cat primary visual cortex (Stepanyants et al., 2008). Probabilistic synaptic density maps were derived from the neural circuitry of the olfactory centers of drosophila to predict connectivity (Jefferis et al., 2007). Drosophila olfactory circuits were mapped and rendered to study the mechanisms of olfactory coding in higher brain centers (Lin et al., 2007).
Significant research efforts have been attempted to extract the geometric structure of axons from microscopic image stacks. Some of the algorithms proposed in literature that are directly related to the work in this article are the Kalman filter-based axon tracking (Jurrus et al., 2006), 3D template matching approach to detect neuronal structures (Al-Kofahi et al., 2002), 3D rayburst sampling for detection of neuron structures (Rodriguez et al., 2006), voxel coding for detection of centerlines from segmented neurons (Vasilkoski and Stepanyants, 2009) and use of generalized cylinders for semi-automatic 3D reconstruction of neurons (Schmitt et al., 2004). As the problems exhibited by our datasets are completely different from the datasets for which these methods were proposed and due to the sheer complexity of the datasets presented in the articl, none of these methods can be directly applied to label individual axons.
The authors' labs have also published their research work in this long-standing problem of extraction and visualization of the geometrical features and topological characteristics of axons. Repulsive force based snake algorithm (Cai et al., 2006) was proposed to solve the close lying axons in a sequence of cross-sectional images. Although this method works for axons that run approximately perpendicular to the cross-sectional plane, it is unsuitable to complicated structures. Probabilistic region growing algorithm (Srinivasan et al., 2007) was proposed to optimize the time taken to track the axons by utilizing a pseudo-3D approach. This method tracks the centerlines of the axons in the maximum intensity projection (MIP) image and switches to the cross-sectional domain in case of ambiguity. This method is suitable for datasets where the axons are somewhat distinguishable in the MIP image and is also unsuitable for use in the datasets presented in this article. Dynamic programming algorithm (Zhang et al., 2008) introduces an approach to track the axons using constraints for smoothness and continuity in the cross-sectional images along with marker-controlled Watershed algorithm. This approach also requires the axons to run in bundles and roughly maintain their orientation throughout the dataset, and hence, cannot be directly applied to our case. Semi-automated reconstruction in serial images (Lu et al., 2009) was introduced as an extension to the already existing software; Reconstruct (Fiala, 2005), to interactively segment the axons in a sequence of cross-sectional images. This approach requires a re-sampling of the datasets so that majority of the axons appear in the cross-sections. This is not always achievable, especially in the datasets presented in this article, as the axons do not exhibit a specific orientation pattern.
Freely available image processing software, such as ImageJ (Abramoff et al., 2004) and Reconstruct (Fiala, 2005), either deal with 2D images/sequence of 2D images in a semi-automatic fashion or fail in case of sharp change in axon orientation, which makes it extremely time consuming and infeasible to track complex structures. On the other hand, commercially available software, such as Imaris (http://www.bitplane.com/go/products/imaris), assumes the dataset to have an excellent contrast and cannot scale up for tracking large neuronal circuitry. Due to the unpredictable orientation change, connection complexity and low-contrast signal of axons, none of current methods or approaches offer an efficient pipeline for tracking and reconstructing wiring diagram of axons in a large-scale, quantitative neurobiology studies.
In this work, we introduce a new tool, AxonTracker-3D, which does not require any parameter adjustment by the user and can drastically reduce the effort and time required to obtain valuable quantitative information from volumetric axon imagery. From our study, the AxonTracker-3D, which incorporates a 3D, real-time tracking approach based on diffusion, greatly reduces the time required in processing the large volumes of axon images from months via manual analysis to minutes on personal computers.
The interscutularis muscles were removed from thy-1-YFP-16 transgenic mice (Feng et al., 2000) for diffraction-limited confocal imaging of entire motor axonal connectomes. For identification of muscle fiber type, muscles were removed and post fixed with 1% PFA for 7 min, frozen and sectioned at 20 μm using a Leica Cryostat (Leica, Germany). Following that, sections were incubated with blocking solution (2% BSA + 1% goat serum + 0.3% triton) at 25○C for 3 h, and incubated with monoclonal antibodies against myosin type I and 2A (mouse anti-myosin I IgG1, 1 : 20, Novocastra, UK; mouse anti-myosin 2A IgG1, 1:10, Iowa Hybridoma Bank, USA) at 4○C for 6–8 h. After several washes in 0.1% PBS-triton, samples were incubated with secondary antibody (Alexa-488 anti-mouse IgG1, 1:1000; Molecular Probes, USA) for 3 h. Finally, muscle sections were rinsed in PBS (25○C, 20 min ×2) and mounted on slides with Vectashield mounting medium.
The motor neurons expressed yellow fluorescent proteins (YFP) uniformly inside each cell. The YFP was excited with a 488 nm line of argon laser using a 60× (1.4NA) oil objective and detected with a 520–550 nm band-pass emission filter. The images were sampled at Nyquist limit in the x − y direction and oversampled by a factor of 1.5 in the z direction (voxel size = 0.1 μm ×0.1 μm × 0.2 μm) with a 12 bit dynamic range. In order to obtain the complete connectome (consisting of ~200 stacks), the motorized stage controlled by the MutliTimeZ macro, developed by Zeiss, was used.
Tracking and segmentation are the two major steps in solving the connectome and extracting quantitative information from these datasets. In this article, AxonTracker-3D, a tool incorporating the tracking and visualization of axons, has been developed to facilitate the connectome function analysis in large-scale, quantitative neurobiology studies. The datasets are processed one-by-one and are then stitched together to form a collage for further analysis. A segmentation approach is also proposed, which uses the tracking results from AxonTracker-3D. The workflow is described in Figure 1.
The user first loads the dataset that needs to be tracked, which is followed by preprocessing for noise removal. The user then selects starting points for all the axons present in the dataset followed by the initiation of the tracking algorithm. The initialization of the algorithm consists of computation of the candidate center points for all the axons present in the dataset, which is done automatically. As the tracking proceeds, the results can be viewed as an overlay on the rendered data in 3D in real time. In case of a problem in tracking due to low contrast, dramatic change in the orientation of the axons or due to very low voxel intensity, the algorithm stops and prompts the user for correction and guidance. The interactive tracking has the provision of deletion of wrongly tracked points and addition of new points to the centerline. After the user guides the algorithm by clicking a few new points along the centerline, the algorithm resumes the tracking process.
Once the tracking is complete, the datasets are stitched together to obtain the big picture (the connectome). The centerlines can be used to segment the axons and visualize them in 3D space. The user interface of AxonTracker-3D, which has been developed using C++, ITK, VTK and QT on a PC with 1.86 GHz Intel® Core™2 CPU and 2 GB of RAM configuration, is shown in Figure 2.
In a few cases when the algorithm fails to find the correct center due to bad contrast between adjacent axons or low intensity, a provision for manual intervention is provided in the software, using which, the user can add or delete the center points of the axons to guide the algorithm. AxonTracker-3D can display the three orthogonal views of the current location of the centerline along with a visualization of the tracking results (Fig. 3). The user can then add or remove voxels from the centerline, hence guiding the software to accurately track the centerlines of the axons. After all the axons in the dataset are tracked, segmentation can be performed.
The raw image datasets obtained are prone to noise and need to be preprocessed. We use anisotropic diffusion image filtering (Perona and Malik, 1990) to reduce the noise while preserving contrast near the edges of the axons. In order to extract the structures from the datasets for quantitative analysis, centerlines or the backbones of the objects have to be extracted first, which can then be used for segmentation.
The process of centerline extraction involves the estimation of most likely voxels in the dataset that belong to the medial axis of the axons. From the initial point (defined manually), we employ a tracking method to extract the centerlines of each axon present in the dataset.
In order to compute the candidate center points of the axons, we define a scoring method for the voxels in the dataset using a combination of the gradient vector flow (GVF) (Xu and Prince, 1997) and object orientation vectors. In their original work, these forces were generated to guide the initial contour to the boundaries of the objects. We use the inverse of the approach where the forces help the boundaries to collapse to the center. To generate these candidate points, we first compute the GVF in all the cross-sectional images along the three orthogonal directions, along the x-, y- and z-axes, respectively. Each cross-sectional image along a particular axis is then binarized and the boundaries of the structures present in the cross-section are shrunk based on the forces of the GVF. An example is shown in Figure 4. As seen in Figure 4c, the resulting scores contain noise around the actual center and hence, non-maxima suppression is used.
Once the computation of possible centers along all three directions is complete, we merge them together using a directional vector-based approach. The entire dataset is binarized and is first divided into smaller subsets. In each such portion, the geometrical orientation vector of the binarized objects contained in the small dataset is computed. The orientation of an object in 3D space can be computed using moment of inertia matrix. The elements of the matrix are the moment of inertia the mass exhibits in 3D space about different combinations of the three orthogonal axes (x, y and z). By computing the eigensystem of the matrix and by looking at the eigenvector with the least eigenvalue, we can determine the direction in which the mass exhibits the least amount of inertia. In other words, this direction is nothing but the orientation of the mass present inside the sub-volume. We first compute the moment of inertia matrix as:
where Ixx = ∑rR V (r) (d(r)y2 + d (r)z2) and similarly Iyy and Izz are defined; Ixy = Iyx = −∑rR V (r)d (r)xd(r)y and similarly, Iyz=Izy and Izx=Ixz are defined in a symmetrical way; R is the set of all the voxels in the binarized volume under consideration; d(r)=[d(r)x, d(r)y, d(r)z] is the distance of the voxel, r, from the center of mass of the binarized volume; and V(r) is the value of the voxel at r in the binarized volume.
We compute the eigenvalues and the corresponding eigenvectors of the system and select the eigenvector with the lowest absolute eigenvalue as the orientation of the object in the sub-volume. Intuitively, based on the projection of this vector along the three axes (x, y and z), the three scores along orthogonal directions are combined using:
where is the orientation vector of the sub-volume. This score is an indication of likelihood of a voxel belonging to the axon centerline. Figure 6b–d are examples of candidate center points calculated using GVF along the x, y and z-axes for a dataset. We can see these candidate center points are very close to actual center points even in a complicated dataset shown in the figure. Figure 6e is obtained after combining the scores in (b–d).
Once the candidate centerlines are computed, the tracking process starts with the manually chosen initial point for each axon. The centerlines of the axon are traced using a moving sphere whose center always lies on current center of the centerline that has been traced so far. To compute the next center on the centerline, the directional vector for the sphere (computed using the current center and nine previous centers) is computed. The radius of the sphere in our study has been set to five voxels. The search for the next center is performed in a preset angle range (set to 45○ in our case) about the directional vector inside the sphere so that the centerline does not trace back itself. In each step, inside the search region, the voxel with the highest GVF score is chosen as the next center. If many voxels are found to have the same value, we resolve them using:
where Ci+1 is the (i+1)-th center (or the new center on the centerline to be determined); R is the set of all the candidate voxels that need to be resolved; Gr is the normalized value of the maximum value of gradient along the path connecting Ci (current center on the centerline) and r; which is a smoothness term based on the dot product of the directional vectors (the directional vector of the sphere and the directional vector formed between the current center and the voxel under consideration respectively); a and b are the weights that denote accuracy and smoothness of the centerline. In our study, these parameters have been set as 0.65 and 0.35, respectively. An illustration of the tracking algorithm is shown in Figure 5. The shaded conical region inside the spheres shown in Figure 5 are the search region at different points along the centerline.
In many cases, axons twist and turn, and lie close to each other in the image stack. The challenges in segmentation include weak boundaries between close-lying axons and large intensity variations in some axons, which cause many segmentation methods to fail. In order to deal with this problem, the segmentation scheme consists of two steps.
In the first step, local threshold 3D region growing is applied to get the rough segmentation result which will be used as the initial segmentation for the subsequent step. The local region is defined as spheres with certain radius (set to two to five voxels in our study), and the centers of the spheres are located on these center points. This radius was determined based on the radius of the thinnest axon present in the datasets.
Once the rough segmentation of axons is obtained, an interactive level set model (Yan et al., 2008) was adopted to detect boundaries in 3D space. The interaction involves two types of mechanisms: repulsion and competition. The repulsion term is for separating the touching axons and the competition is for defining the axon boundaries. These mechanisms were formulated as an energy functional, which is then minimized by using the multiphase level set method. Figure 6g shows the segmentation result obtained using the centerlines of the axons from Figure 6f.
After the datasets have been segmented, they are stitched together to form a collage. As the datasets are contiguous with a variable amount of overlap between them, we process the dataset stack by stack to form a montage by using the normalized cross-correlation measure to determine contiguous stacks and the extent of overlap between them. The cross-correlation measure between an image, I, and a template, t, is defined as:
where (u, v) is the coordinate where the template is centered at; I(x, y) is the image region under the moving template; is mean of the template; and Īu,v is the mean of the image region under the template. The cross-correlation measure has been applied in xy plane and then in xz plane to complete the alignment in 3D. This reduces the computational time involved in alignment as opposed to a completely 3D alignment. As the datasets were down-sampled by a factor of 2, an alignment error of one voxel is unavoidable during alignment.
To record morphological properties of neuronal processes, images were acquired from neonatal mice using laser scanning immunofluorescent confocal microscopes. The connectome of interscutularis muscle of the mouse was used for this study. The complete connectome consisted of ~200 data stacks. Adjacent stacks had an overlap of around 10% for alignment to form the montage. The data acquisition lasted around 24–48 h and accumulated ~50 GB of data. The datasets that did not contain any axons were discarded. The resulting set of datasets contained 154 stacks which amounted to around 35 GB of data.
The volumetric datasets are a large number of 3D axon image stacks, each typically being 1024 × 1024 × (100–200) voxels in size (voxel size of 0.1 μm × 0.1 μm × 0.2 μm), which when put together form a collage of the connectome. The datasets were down-sampled along the xy direction, yielding datasets that were 512 × 512 ×(100–200) voxels in size (voxel size of 0.2 μm × 0.2 μm × 0.2 μm), which improved the tracking and segmentation speed without losing accuracy and lowered the computer memory requirement.
Figure 7 shows the tracking and segmentation results of 30 contiguous stacks from the collage. The centerlines were tracked and segmented automatically using AxonTracker-3D. The datasets were then stitched together using the cross-correlation measure and are visualized together as shown in Figure 7b. The results show that touching axons and axons with short turns can be accurately extracted by the software.
On an average, each axon takes <30 s to track if no manual intervention/correction is required. Owing to the complexity of axon, each dataset had ~25–40% of the axons that needed manual correction at three to four places for tracking. The average reconstruction time for one image stack is ~10–15 min. Once tracked, each dataset in the collage took 5 min for segmentation.
In this article, we presented a novel interactive tool to extract and reconstruct axons obtained from 3D confocal microscopy image stacks. Compared to the previously proposed methods, our algorithm has two major contributions: firstly, the software can successfully extract centerlines of axons with minimal user intervention. Rendering of the data in 3D along with overlay of tracking results in real-time is a useful feature offered by the software. By introducing user interaction, we can see that our algorithm can successfully extract and reconstruct touching axons and axons with dramatic orientation changes. Secondly, by using local region growing and interactive level set segmentation methods, we can get the actual boundary of each axon which can provide us valuable quantitative information such as the length, radius for axon function and connectivity analysis. As a continuation of the work presented in this article, we will discuss algorithms to detect and track branches in order to minimize user-intervention needed for AxonTracker-3D. The following are the planned improvements to the current version of the software.
Axons in the connectome of neuro-muscular junctions tend to branch into complex patterns. The current version of AxonTracker-3D does not have a provision to automatically detect branches. At branch points, the centerlines were joined together manually to complete the wiring circuitry of the axons. As an extension to our current work, we plan to design and implement new algorithms to detect branching points along the axons.
Methods such as AdaBoost classifier (Freund and Schapire, 1999) could be a potential solution to this problem. This consists of three steps: (i) re-slice the axon tubes along its orientation; (ii) extract 2D and 3D features from the slices and spheres rounding the center points; (iii) select samples to train AdaBoost classifier (manual selection of ~100 samples for training). The AdaBoost method can be used to classify positive training examples from negative examples by selecting a small number of critical features from a huge feature set previously designed and creating a weighted combination of them to use as a strong classifier to detect the branch points along axons.
We will include a provision for manual intervention in branch detection as well, but the need for it will be rather minimal. AxonTracker-3D has been designed such that the user interface and interaction is kept as simple, intuitive and user-friendly as possible. Any addition or modification to AxonTracker-3D modules will not affect the basic work-flow as described in Figure 1.
The current version of the software does not make use of multi-threading, and hence there is a good scope of increasing the execution speed manifolds by making use of multiple CPU's if available on the PC. As the size and number of the datasets in the collages are huge, a major improvement to the existing software would be to include a provision for batch processing. The datasets in the collage could be processed automatically as much as possible to reduce the time and effort needed for tracking the connectome. Along with this, the existing algorithms will also be improved to increase automation without compromising accuracy of results.
Most importantly, using the quantitative data available from tracking and segmentations of complete connectomes, many biological questions can be generated, such as: (i) the relationship between the number of branches and the length of each segment in a population of axons; and (ii) the spatial distribution of motor neurons or the optimal layout of the motor neurons in a population of axons. Questions such as whether the spatial distribution of the axons are random in nature or follow a certain pattern can be answered. We intend to employ data mining and biostatistics methodologies on the datasets and wiring diagrams extracted from AxonTracker-3D to address these questions.
The authors would like to thank various members of the Wong lab at The Methodist hospital Research Institute and of the Lichtman lab at Harvard University for their technical comments over the years on the AxonTracker-3D project.
Funding: This research is funded by grants from Ting Tsung and Wei Fong Chao Center for Bioinformatics Research and Imaging in Neurosciences (BRAIN), The Center for Bioengineering and Informatics at The Methodist Hospital Research Institute, and Harvard Neurodiscovery Center (previously Harvard Center for Neurodegeneration and Repair) (Wong).
Conflict of Interest: none declared.