Skip to main content

METHODS article

Front. Neuroanat., 25 November 2014
Volume 8 - 2014 |

Automated computation of arbor densities: a step toward identifying neuronal cell types

  • 1Department of Brain and Cognitive Sciences, Massachusetts Institute of Technology, Cambridge, MA, USA
  • 2Department of Ophthalmology, Harvard Medical School, Boston, MA, USA
  • 3Department of Electrical Engineering and Computer Science, Massachusetts Institute of Technology, Cambridge, MA, USA
  • 4Department of Neurobiology, Harvard Medical School, Boston, MA, USA
  • 5Princeton Neuroscience Institute and Computer Science Department, Princeton University, Princeton, NJ, USA

The shape and position of a neuron convey information regarding its molecular and functional identity. The identification of cell types from structure, a classic method, relies on the time-consuming step of arbor tracing. However, as genetic tools and imaging methods make data-driven approaches to neuronal circuit analysis feasible, the need for automated processing increases. Here, we first establish that mouse retinal ganglion cell types can be as precise about distributing their arbor volumes across the inner plexiform layer as they are about distributing the skeletons of the arbors. Then, we describe an automated approach to computing the spatial distribution of the dendritic arbors, or arbor density, with respect to a global depth coordinate based on this observation. Our method involves three-dimensional reconstruction of neuronal arbors by a supervised machine learning algorithm, post-processing of the enhanced stacks to remove somata and isolate the neuron of interest, and registration of neurons to each other using automatically detected arbors of the starburst amacrine interneurons as fiducial markers. In principle, this method could be generalizable to other structures of the CNS, provided that they allow sparse labeling of the cells and contain a reliable axis of spatial reference.

1. Introduction

The classification of neuronal types is far from complete. Advances in genetic engineering for sparse and specific labeling (Gong et al., 2003; Wickersham et al., 2006, 2007; Kim et al., 2008; Chung et al., 2013; Ke et al., 2013) offer improved data acquisition and molecular identification of neuronal classes. However, the need for structural information has not diminished because what defines a true neuronal type is not clear when only molecular information is available. One challenge facing a successful classification is to ensure that every cell type is represented in the sample set. For the structural approach, dense reconstruction of tissues imaged by electron microscopy offers a solution to this completeness problem (Denk and Horstmann, 2004; Hayworth et al., 2006; Bock et al., 2011). On the other hand, electron microscopy is not yet capable of either obtaining large enough sample sets to capture the biological variability within individual cell types, or imaging cells with very large neuronal arbors. Light microscopy offers high throughput imaging and a large field of view to complement electron microscopy. However, the time-intensive tracing step represents a bottleneck of the overall program.

Recently, it was shown that neurons in the mammalian retina can achieve submicron precision in their laminar positioning (Sümbül et al., 2014). This was done by combining an arbor density formalism (Stepanyants and Chklovskii, 2005) with a neurite based registration system for sparsely labeled neurons. The ensuing arbor density classification suggests that a robust classification of all mammalian retinal ganglion cells is within reach. However, this study and many other previous attempts (Sun et al., 2002; Badea and Nathans, 2004; Kong et al., 2005; Coombs et al., 2006; Völgyi et al., 2009) depend on manual tracing of individual neuronal arbors, which is a time-intensive task. Tracing a neuronal arbor creates a “skeleton representation” of the arbor, which consists of interconnecting line segments going through the dendrites. The thickness of dendrites along the line segments is often ignored because tissue preparation artifacts can result in unreliable estimates. In contrast, a volumetric representation includes both the skeleton and the dendrite thickness along the skeleton. Here, we propose an automated method using volumetric analysis to aid the classification of neuron types. At the heart of our approach is the simple observation that while the arbor density representation in Sümbül et al. (2014) requires a precise characterization of laminar positioning, it does not utilize detailed descriptions of arbors. In particular, we demonstrate that volumetric stratification precision of neurons can match the trace-based precision in the mammalian retina. Our method is designed for sparse imaging scenarios. It does not address the problem of separating the arbors of overlapping neurons from each other, for which tracing may still be required. Kim et al. (2014) recently used volumetric analysis and semi-manual arbor reconstruction to identify bipolar and starburst amacrine cells in an electron microscopy setting.

Volumetric reconstruction of neuroanatomy from an image stack involves obtaining a digital representation of the neuronal arbor (i.e., a voxel is “white” if it belongs to the cell, and “black” otherwise), and registering this representation to other neuronal structures to achieve a comparative description. As a first step to reconstruct a sparsely labeled neuron, we use a convolutional network (LeCun et al., 1998), which is a supervised machine-learning architecture, to enhance the image quality and suppress the acquisition noise. Although robust and accurate reconstruction of neuronal morphology is still a largely unsolved problem, it has become a bottleneck only recently as a result of the advances in high-throughput imaging. The demand for automated reconstruction prompted the Digital Reconstruction of Axonal and Dendritic Morphology Challenge (DIADEM challenge) (Brown et al., 2011). The challenge helped disseminate many novel approaches (Bas and Erdogmus, 2011; Chothani et al., 2011; Narayanaswamy et al., 2011; Turetken et al., 2011; Wang et al., 2011; Zhao et al., 2011). We anticipate that some of these approaches may be preferable to the convolutional network module of our method depending on the imaging conditions. A common problem is that when labeling is not sparse enough, cells other than the neuron of interest are also reconstructed. Our solution is to apply a post-processing routine to remove extraneous objects after the initial reconstruction step.

In the mammalian retina, the dendrites of the starburst amacrine interneuron form two parallel surfaces in the inner plexiform layer, which serve as fiducial marks (Haverkamp and Wässle, 2000). When the tissue is not flattened to preserve internal structure, it assumes a wavy form under the microscope. We solve this problem by digitally flattening (unwarping) the stack with the guidance of starburst surfaces after the imaging is done. Finally, we obtain a common depth coordinate by registering the starburst surfaces from different stacks to each other.

2. Materials and Methods

2.1. The Dataset

We use the retinal ganglion cells (RGCs) from a recent study on the classification of retinal cell types (Sümbül et al., 2014). The dataset was obtained by confocal microscopy at a voxel size of 0.4 μm × 0.4 μm × 0.5 μm. This dataset also includes the relative positions of On and Off starburst amacrine interneurons for each RGC, by staining for choline acetyltransferase (Haverkamp and Wässle, 2000), thereby allowing a stratification analysis of RGCs based on starburst amacrine arbors. The methodological bottleneck of that study was the semi-automated tracing of RGC arbors, which required an average time of 40 min per trace with experienced tracers. The full dataset includes five strongly defined cell types, which have consistent and specific functional, molecular, and structural identifiers. We focus here on this subset, and omit the stacks where labeling is too dense (i.e., existence of many neurites in close proximity from more than one neuron) or too dim for fully automated analysis. In a few cases, the starburst surfaces were weakly stained; these were also omitted. After this culling, two neuron types did not have enough representatives for statistical analysis and were omitted altogether. The final dataset comprises 50 neurons that form three molecularly, physiologically, and structurally homogeneous cell types.

The JAM-B neurons express the junction adhesion molecule JAM-B, respond to offset of upward moving stimuli, and their arbors are asymmetric in the dorsal-ventral axis (in the central retina) (Kim et al., 2008). The W3 neurons express the TYW3 gene, are sensitive to local edges, and have one of the smallest arbor sizes in the mammalian retina (Kim et al., 2010). The BDa neurons express the FSTL4 gene, are On-Off direction sensitive, and arborize twice (Kim et al., 2010). Finally, these cell types are known to stratify at characteristic depths in the inner pexiform layer with submicron precision [distance from the On starburst surface: 15.6 μm (JAM-B), 5.5 μm (W3), 0.3 μm (BDa)—BDa neurons stratify again 0.3 μm distal to the Off starburst surface] (Sümbül et al., 2014).

2.2. Volumetric Reconstruction of Sparsely Labeled Neurons from Manual Traces

We use the concept of simple pixel from digital topology (Bertrand and Malandain, 1994) to probe whether neuronal mass attains the stratification precision achieved by the arbor traces (skeletons). A simple pixel is defined as a pixel that does not change the topology of the digital image when its value is flipped. (i.e., does not create/remove objects, holes, splits, mergers) Similar approaches were previously used in the reconstruction of dense electron microscopy images of neuronal tissue (Jain et al., 2010; Helmstaedter et al., 2013). Specifically, we inflate the individual traces by respecting the topology of the traces (via simple pixel characterization), and the geometry of the neurons (via thresholding the brightness values in the raw image). We use 60% of the maximum brightness value in an image stack as the threshold. We iterate the inflation process 62 times, potentially inflating by a single layer of voxels at each step so that somata as large as (62 × 2 + 1) × 0.4 μm= 50 μm in diameter are properly characterized. Algorithm 1 presents a pseudocode of the steps. The resulting three dimensional binary stacks are seemingly perfect characterizations of neuronal morphology based on the raw image stacks and the arbor traces (Figure 1) because they respect both the tree structure (through tracing, Figure 1B), and the dendritic widths (through inflation, Figure 1D). The caveat is that the resulting volumetric representations depend on the laborious task of (semi-) manual tracing.


Algorithm 1. Pseudocode for topologically constrained inflation of a trace. Binary operations on same-size arrays are to be interpreted elementwise. ¬ (⊕) denotes negation (exclusive-or). imdilate dilates its first argument using its second argument as the kernel. nnz returns the number of nonzero (true) entries in an array. Matlab notation is used in the array on line 12.


Figure 1. Volumetric reconstruction of a BDa neuron starting from a trace. Maximum intensity projections of the raw image stack (A), the manually traced arbor (B), the inflated trace after one round of topologically constrained inflation (C), and the inflated trace after 62 rounds of topologically constrained inflation (D). In each panel, large image: xy projection, bottom: zy projection, right: xz projection. Scale bar, 40 μm; bottom-right, xy projection image in panel D.

2.3. Automated Enhancement and Post-Processing of RGC Arbors

Various approaches have been developed recently for automated reconstruction of neuronal morphology from sparsely labeled image stacks (Al-Kofahi et al., 2002, 2008; Schmitt et al., 2004; Zhang et al., 2007; Losavio et al., 2008; Peng et al., 2010, 2011; Srinivasan et al., 2010; Bas and Erdogmus, 2011; Turetken et al., 2011; Wang et al., 2011; Xie et al., 2011; Choromanska et al., 2012; Turetken et al., 2012; Gala et al., 2014). While these methods can capture the geometrical layout of neuronal arbors, imperfections in tissue handling and imaging (e.g., non-uniform labeling of neurites, high density labeling, low signal-to-noise ratio images) often result in topological errors such as missing branches and extraneous structures. On the other hand, blurring and projection operations are robust against local mistakes. Therefore, topological imperfections in the reconstruction may be acceptable for cell type identification purposes so long as the general morphology of a neuron is captured properly. As a first step, we use the convolutional network based enhancement of RGC arbors reported in (Sümbül et al., 2014). A convolutional network is a feed-forward network of convolutional filters whose outputs are transformed by a non-linearity (e.g., sigmoid). An advantage of such a supervised machine learning approach is that it does not have free parameters to adjust. Rather, the paradigm depends on the existence of a labeled training set through which the various parameters are automatically optimized. The network is trained to transform noisy gray-scale images of sparsely labeled neurons into cleaner binary images. Here, we improve the architecture and filter sizes, and provide an efficient implementation that does not need specialized hardware ( The resulting network has 8 layers with 8 perceptrons in each hidden layer except for the last hidden layer, which is a fully connected layer of 100 perceptrons. The filter sizes within each layer are identical and are as follows: 5 × 5 × 1, 5 × 5 × 1, 3 × 3 × 3, 5 × 5 × 1, 3 × 3 × 3, 3 × 3 × 3, 1 × 1 × 1, 1 × 1 × 1. Therefore, the overall patch size to decide whether the central voxel of the patch belongs to a neurite or not is 19 × 19 × 7 voxels (7.6 μm × 7.6 μm × 3.5 μm). The network has all-to-all connectivity between subsequent layers, and is trained by backpropagation learning LeCun et al. (1998).

When the density of labeling is not low enough, somata and neurites of other neurons appear in the image stacks. On the other hand, the reconstructed arbors may have breaks due to dim/inhomogeneous labeling. Therefore, we devise a simple post-processing routine to isolate the neuron of interest. The algorithm uses connected component analyses and basic morphological image operations (i.e., opening and dilation) to remove extraneous structures and somata. In particular, the algorithm detects the largest object in the image stack and removes the objects that are smaller than a given size and farther from the largest object than a given distance. Somata are removed by locating and removing the white regions that are large enough to fully enclose a given ellipsoid (Algorithm 2). While soma size is known to carry information on neuronal identity, it is a weak classifier (Sun et al., 2002; Coombs et al., 2006; Völgyi et al., 2009). The final image stacks may include axonal projections from other neurons, imperfectly suppressed noise, missing small branches, extraneous branches from other neurons, and splits/mergers of the neuronal arbor depending on the image quality and the sparsity of labeling in the tissue. Nevertheless, the next few subsections demonstrate that the reconstruction quality is high enough to study stratification patterns and probe neuronal identity.


Algorithm 2. Pseudocode for post-processing a binary volume. Various binary operations are as defined in Algorithm 1. Matlab notation is used for brevity. bwlabeln returns an array the same size as its argument, where voxels are assigned different values iff they belong to different objects (26-connectivity). regionVolumes returns a list of object sizes. bwareaopen removes from its first argument all objects whose volumes are smaller than the second argument. imopen performs a morphological opening operation on its first argument using a cubic kernel whose edge length is given by the second argument.

2.4. Quasi-Conformal Unwarping of Volumetric Data and Laminar Registration

We use the automatically detected starburst surfaces in individual stacks as fiducial marks (Figure 2). We find quasi-conformal mappings that independently transform the detected starburst surfaces into flat surfaces as described in Sümbül et al. (2014) to maximally preserve local angles within the surfaces (Levy et al., 2002). The two flattened surfaces are registered to each other in-plane by matching the xy coordinates of the patch in which both starburst layers are the flattest. We extend the resulting transformation to other points in the image stack by using local polynomial approximations (quadratic in xy, linear in z). In particular we apply the transformation to individual voxels of the binary three-dimensional representation of a neuron, rather than its trace points. The transformed voxels are scaled and shifted in z so as to place the flattened On starburst surface at z = 0 μm and the flattened Off starburst surface at z = 12 μm. Figure 3 depicts the dramatic effect of unwarping on a BDa neuron. Finally, the histogram of depth positions of the voxels (depth profile) is obtained by gridding onto a Cartesian grid with a resolution of 0.5 μm (Figure 3D). The gridding step uses a Kaiser-Bessel kernel (Jackson et al., 1991) to maintain accuracy in laminar registration, and applies weights to individual voxels to compensate for the distance fluctuations between warped voxels. Note that if the arbor density function is obtained by blurring in xy only (and not in z), then the depth profile is the projection of the three-dimensional arbor density function.


Figure 2. Fully automated enhancement and post-processing of RGC arbors (green), and detection of starburst surfaces (red). Left: xy (A), xz (B), and zy (C) projections of the raw image of an RGC. Right: xy (D), xz (E), and zy (F) projections of the processed arbors and detected starburst surfaces of the same RGC. Starburst surfaces within a slab are shown and starburst somata are removed for better visualization. Scale bar, 40 μm; lower-right, panel D.


Figure 3. Warping neuronal mass reveals volumetric depth profiles. Left: xy (A), xz (C), and yz (E) projections of a BDa cell reconstruction obtained by inflating its trace. Right: xy (B), xz (D), and yz (F) projections of the reconstructions after soma removal and quasi-conformal unwarping of the white voxels. Note the depth alignment of neurites after warping. (G) The normalized depth profiles of the RGC based on its trace (skeleton) and trace-based volumetric reconstruction. Scale bar, 40 μm; bottom-right, panel B.

2.5. Statistical Measures and Other Metrics

The peak position of a depth profile is the signed distance from the On starburst layer at which the profile achieves its maximum value. [The On (Off) layer is located at z = 0 μm (z = 12 μm.)] For the bistratified BDa cells, a second peak position is also reported. This second peak position is defined as the depth value at least 6 μm away (half the distance between starburst layers) from the first peak position, at which the remaining profile achieves its maximum value.

We assume that the peak positions of the depth profiles of individual neurons of a given type are independent and identically distributed (i.i.d.) with N(μ, σ2). The distribution of the sample variance of n i.i.d. N(μ, σ2) observations is given by χ2n − 1(t(n − 1)/σ2)(n − 1)/σ2, where χ2n − 1(t) denotes the chi-squared distribution with n − 1 degrees of freedom. The symmetrical 95% confidence interval for σ, given the sample standard deviation s, is

[(n1)s2Xn12(0.975),(n1)s2Xn12(0.025)],    (1)

where X−2n − 1 denotes the inverse cumulative distribution function of χ2n − 1.

We use the Brown-Forsythe test to infer whether the different reconstruction methods return significantly different variance values for the peak stratification position of cells of the same type.

We define the “signal” in each normalized depth profile for cells of a given type as the average normalized profile over cells of that type. Then, the “noise” in each profile is the difference of the profile from the “signal” component. We define the signal-to-noise ratio (SNR) for a cell type as the average, of the Euclidean norm of the signal divided by the Euclidean norm of the noise, over all cells of that type.

The Crest factor is defined as the peak amplitude of a profile divided by the root-mean-square value of the profile. That is, it is the ratio of the peak value to the average value. It indicates how extreme the peak is in a given depth profile. Since narrow and sharp peaks in a profile where the “background” regions are small makes it easy to detect a cell type in the presence of many cell types, we use the Crest factor as a figure of merit for the different approaches analyzed in this paper.

3. Results

3.1. Projections of Volumetric Data Preserve the Stratification Precision of RGCs

We obtain the volumetric reconstructions of all 50 neurons in the dataset by inflating their manually reconstructed traces as described in Algorithm 1. These volumetric reconstructions were unwarped and registered, to obtain depth profiles of all neurons in the dataset. Figure 4 shows the average profiles for each neuron type generated from the volumetric reconstructions together with the average profiles of the traces (skeletons). Two qualitative observations emerge: (i) The peak positions of the average profiles are preserved across the two methods. The distribution of mass along the skeletons of neurons preserve the peak stratification depth of the skeletons. (ii) The peaks of the normalized profile averages are lower in the volume-based approach because the branches close to the soma are typically thicker than the distal branches. Table 1 tabulates the Crest factors for both methods, and quantifies the observation that the trace profiles have slightly sharper peaks. Profiles with sharper peaks are preferable when identifying cell types in the presence of a heterogeneous dataset, similar to spectroscopy.


Figure 4. Depth profiles of trace-based volumetric reconstructions maintain the stereotypy attained by the depth profiles of arbor traces while having lower peaks. (A) Depth profiles of the arbor traces, (B), depth profiles of the topology preserving inflations of the traces.


Table 1. Mean and standard deviation values for the peak positions and norms of the depth profiles.

The specificity of stratification peaks -not just their average- is important to be able to identify cell types. The sample standard deviations of the peak position for each cell type do not change significantly between the skeleton-based and volume-based depth profiles (Brown-Forsythe test—See Tables 1, 2 for individual n and p-values). This suggests that neurons of a given type are as precise about distributing their neuronal volume in depth as they are about distributing their skeleton-based presence.


Table 2. Statistical measures of the variability in peak positions.

Table 1 also tabulates the mean SNR values over the three cell types using both the trace profiles and the trace-based volumetric profiles (Methods). Higher SNR values indicate stereotypical distributions. While n = 3 pairs are too few to probe statistical significance with rigor, the trace-based volumetric profiles do not seem to have lower SNR values than the trace profiles. (Right-sided Wilcoxon signed rank test, n = 3 pairs, p = 0.75)

3.2. An Automated Method to Probe RGC Identity

The stereotypy of the profiles of volume-based reconstructions obtained by inflating manual traces suggests that it may be possible to avoid the laborious task of manual tracing altogether for cell type identification purposes. For comparison, we begin by implementing simple thresholding: Each image stack is thresholded at 60% of its maximum brightness value. Then, somata are removed and the resulting binary stack is unwarped and registered as described in the Methods. The results are not impressive: The extraneous structures and imaging artifacts contaminate many stacks significantly. As a simple proxy, we observe that the mistakes perturb the depth profiles enough to create spurious peaks far away from the original stratification peak in 11 out of 50 stacks (Figure 5A). Therefore, this simple approach is not suitable for automation. This is especially clear for cells that stratify close to the ganglion cell layer. After removing the part of the profiles where z < −6 μm, we find that the stratification stereotypy of the depth profiles is essentially preserved (Tables 1, 2). Nevertheless, even in this restricted region, the presence of objects that do not belong to the neuron makes type identification harder, as reflected by the low SNR values and the Crest factors (Table 1 and Figure 5C).


Figure 5. Enhancement and post-processing enable automated peak detection and produce sharper depth profiles with higher SNR values. (A) Normalized depth profiles obtained by thresholding the raw stack display large build-ups close to the ganglion cell layer that create spurious profile peaks and low SNR values. (B) Enhancement and post-processing remove objects not belonging to the neuron of interest while retaining the neuron to produce depth profiles with consistent peak positions and higher SNR values. C, Peak values vs. peak positions of all four methods. Symbols indicate the mean values, and the lower/upper bounds indicate the 10th and 90th percentiles. Each color indicates a different method as defined in the legend, and each panel depicts a single cell type as indicated at the top. Both peaks are shown for the bistratified BDa neurons. Peaks at z < −6 μm are not considered for the threshold method.

While the threshold based approach may allow for cell type identification when the image stacks are sparsely labeled and have very low noise, insufficient suppression of the background noise and failure to isolate the neuron of interest from other structures prevents it from working reliably on our dataset. Therefore, we apply the convolutional network described in the Methods on the image stacks to suppress the background noise, retain the neuronal structures, and connect the occasionally disconnected neurite pieces. Subsequently, we apply the post-processing routine (Methods) to remove the extraneous structures from the image stack that are not critically close to the neuron of interest. Notably, no manual labor is used in this scheme.

A drawback is that the automated approach occasionally causes splits and mergers in the reconstruction and includes extraneous structures. On the other hand, the depth profiles—one-dimensional arbor densities that serve as proxies for the three-dimensional arbor densities—identify the stratification peaks correctly (Figure 5B). Moreover, the sample standard deviation of the peak position did not change significantly in any of the three neuron types (Brown-Forsythe test—See Tables 1, 2 for individual n and p-values). The Crest factors for this automated method are lower than those of the trace profiles, but they are roughly the same as those of the trace-based volumetric profiles. Lastly, the mean SNR value for the automated method is lower than that for the trace-based approaches, but it is higher than the threshold method's mean SNR value (Table 1).

4. Discussion

Identifying and providing experimental access to homogeneous cell types of nervous systems is a prerequisite to understanding the fundamental principles of brain function in health and disease. Recently, it was shown that a method using a neurite based registration system and an arbor density representation of neurons is capable of robustly identifying the mammalian RGC types in a highly heterogeneous sample set (Sümbül et al., 2014). Notably, that study relied on traces of neuronal arbors, which are time consuming to obtain. Here, we show that the spatial distribution of the arbor volume attains a stratification precision similar to that of the arbor trace. Based on this observation, we describe an automated method that can remove the time intensive tracing step in identifying cell types. We anticipate our approach to be useful in integrating structural information to studies that investigate the molecular or functional dynamics of neurons, as well as purely anatomical pursuits.

We quantify the stratification precision as the standard deviation of the peak position of the depth profiles. We do not observe significant differences between the stratification precisions of the depth profiles of the traces and the volumes obtained by inflating the traces or by our automated method (Table 2). This suggests that the depth distribution of the overall mass can be as stereotyped as that of of the skeletal mass. Another observation suggesting volumetric stereotypy is the lack of a significant difference between the mean SNR values of the normalized depth profiles of the traces and the volumes obtained by inflating the traces.

We have argued that the presented method can be useful in identifying cell types using three-dimensional arbor densities. However, we have not attempted a formal classification of the cells used in this study. While Figure 5C, Tables 1, 2 clearly suggest that such an attempt would be successful, classification becomes a hard task only in the presence of a highly heterogeneous dataset. On the other hand, considering that the automated approach can maintain the stratification precision attained by the trace based analysis and the arbor density representation in Sümbül et al. (2014) used substantial in-plane blurring (and no axial blurring), it is plausible that arbor densities generated from the output of our automated method can be classified successfully not only in the presently studied dataset of three cell types, but also in a more heterogeneous sample set.

We observe that the peak values of the normalized volumetric profiles are smaller than those of the normalized trace profiles. This can be explained by the fact that branches closer to the soma are typically thicker than the distal branches, presumably to minimize signal propagation delays while keeping arbor volume to a minimum (Chklovskii and Stepanyants, 2003).

Dim or inhomogeneous labeling of neurites, denser (not sparse enough) labeling of neurons, and high noise levels often result in imperfect reconstructions with the current state-of-the-art automated approaches. Our convolutional network implementation is not immune to such imperfections, either. Removal of failing image stacks decreases the throughput of the overall method. On the other hand, standard approaches in machine learning, such as boosting and training deeper networks with larger training data, suggest ways of increasing the throughput by providing better noise suppression and better reconstruction of arbor topology. Moreover, while other automated reconstruction methods often require manual tuning of free parameters, they can be inserted instead of our convolutional network implementation as well (Al-Kofahi et al., 2002, 2008; Schmitt et al., 2004; Zhang et al., 2007; Losavio et al., 2008; Peng et al., 2010, 2011; Srinivasan et al., 2010; Bas and Erdogmus, 2011; Turetken et al., 2011, 2012; Wang et al., 2011; Xie et al., 2011; Choromanska et al., 2012; Gala et al., 2014).

While we investigate retinal ganglion neurons in this study, our approach only assumes (i) the existence of an arbor marker specific to a cell type and (ii) a method of labeling cells sparsely in a laminar structure. Therefore, it is readily extendible to other neuron classes of the retina. In particular, the same fiducial marks (starburst amacrine cells) and very similar sparse labeling methods can be used to study the classification and co-stratification of bipolar and amacrine cell classes. The effort required to trace a neuron increases as the complexity of its arbor increases. Hence, the potential impact of our method is higher for neurons whose total dendritic lengths are larger. Cortical neurons are typically much larger than retinal neurons, and classifying them is an impending problem (Ascoli et al., 2008). Traditionally, obtaining datasets of cortical neurons that capture their diversity has been a practical challenge. However, recent advances in tissue clarification and a multiplicity of genetic or viral methods (Gong et al., 2003; Wickersham et al., 2006, 2007; Kim et al., 2008; Chung et al., 2013; Ke et al., 2013) enable high-throughput structural imaging of such neurons. Therefore, we speculate that our approach can be useful in automating the discovery and identification of cortical cell types if the two requirements mentioned above are met.

Data Sharing Statement

The ZNN library is available at The subset of trace files, the associated automatically detected starburst surfaces, the trace based reconstructions, and the software used in this study are available at The original dataset of trace files is available at Raw image stacks are available upon request.

Author Contributions

All authors contributed to the conception of the study and the editing of the manuscript. Uygar Sümbül designed and performed the analysis. Aleksandar Zlateski conceived and implemented the ZNN package. Uygar Sümbül and Ashwin Vishwanathan wrote an initial draft of the manuscript.


We are grateful for financial support from the Harvard Neuro-Discovery Center, the U.S. Army Research Office (W911NF-12-1-0594), DARPA (HR0011-14-2-0004), NIH/NINDS, the Howard Hughes Medical Institute, the Gatsby Charitable Foundation, and the Human Frontier Science Program.

Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.


Al-Kofahi, K. A., Lasek, S., Szarowski, D. H., Pace, C. J., Nagy, G., Turner, J. N., et al. (2002). Rapid automated three-dimensional tracing of neurons from confocal image stacks. IEEE Trans. Inf. Technol. Biomed. 6, 171–187. doi: 10.1109/TITB.2002.1006304

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Al-Kofahi, Y., Dowell-Mesfin, N., Pace, C., Shain, W., Turner, J. N., and Roysam, B. (2008). Improved detection of branching points in algorithms for automated neuron tracing from 3D confocal images. Cytometry A 73, 36–43. doi: 10.1002/cyto.a.20499

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Ascoli, G. A., Alonso-Nanclares, L., Anderson, S. A., Barrionuevo, G., Benavides-Piccione, R., Burkhalter, A., et al. (2008). Petilla terminology: nomenclature of features of gabaergic interneurons of the cerebral cortex. Nat. Rev. Neurosci. 9, 557–568. doi: 10.1038/nrn2402

CrossRef Full Text | Google Scholar

Badea, T., and Nathans, J. (2004). Quantitative analysis of neuronal morphologies in the mouse retina visualized by using a genetically directed reporter. J. Comp. Neurol. 480, 331–351. doi: 10.1002/cne.20304

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Bas, E., and Erdogmus, D. (2011). Principal curves as skeletons of tubular objects: locally characterizing the structures of axons. Neuroinformatics 9, 181–191. doi: 10.1007/s12021-011-9105-2

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Bertrand, G., and Malandain, G. (1994). A new characterization of 3-dimensional simple points. Pattern Recogn. Lett. 15, 169–175. doi: 10.1016/0167-8655(94)90046-9

CrossRef Full Text | Google Scholar

Bock, D. D., Lee, W. C. A., Kerlin, A. M., Andermann, M. L., Hood, G., Wetzel, A. W., et al. (2011). Network anatomy and in vivo physiology of visual cortical neurons. Nature 471, 177–182. doi: 10.1038/nature09802

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Brown, K., Barrionuevo, G., Canty, A., De Paola, V., Hirsch, J., Jefferis, G., et al. (2011). The diadem data sets: representative light microscopy images of neuronal morphology to advance automation of digital reconstructions. Neuroinformatics 9, 143–157. doi: 10.1007/s12021-010-9095-5

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Chklovskii, D. B., and Stepanyants, A. (2003). Power-law for axon diameters at branch point. BMC Neurosci. 4:18. doi: 10.1186/1471-2202-4-18

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Choromanska, A., Chang, S. F., and Yuste, R. (2012). Automatic reconstruction of neural morphologies with multi-scale tracking. Front. Neural Circuits 6:25. doi: 10.3389/fncir.2012.00025

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Chothani, P., Mehta, V., and Stepanyants, A. (2011). Automated tracing of neurites from light microscopy stacks of images. Neuroinformatics 9, 263–278. doi: 10.1007/s12021-011-9121-2

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Chung, K., Wallace, J., Kim, S., Kalyanasundaram, S., Andalman, A., Davidson, T., et al. (2013). Structural and molecular interrogation of intact biological systems. Nature 497, 332–337. doi: 10.1038/nature12107

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Coombs, J., van der List, D., Wang, G.-Y., and Chalupa, L. (2006). Morphological properties of mouse retinal ganglion cells. Neuroscience 140, 123–136. doi: 10.1016/j.neuroscience.2006.02.079

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Denk, W., and Horstmann, H. (2004). Serial block-face scanning electron microscopy to reconstruct three-dimensional tissue nanostructure. PLoS Biol. 2:e329. doi: 10.1371/journal.pbio.0020329

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Gala, R., Chapeton, J., Jitesh, J., Bhavsar, C., and Stepanyants, A. (2014). Active learning of neuron morphology for accurate automated tracing of neurites. Front. Neuroanat. 8:37. doi: 10.3389/fnana.2014.00037

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Gong, S., Zheng, C., Doughty, M., Losos, K., Didkovsky, N., Schambra, U., et al. (2003). A gene expression atlas of the central nervous system based on bacterial artificial chromosomes. Nature 425, 917–925. doi: 10.1038/nature02033

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Haverkamp, S., and Wässle, H. (2000). Immunocytochemical analysis of the mouse retina. J. Comp. Neurol. 424, 1–23. doi: 10.1002/1096-9861(20000814)424:1<1::AID-CNE1>3.0.CO;2-V

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Hayworth, K. J., Kasthuri, N., Schalek, R., and Lichtman, J. W. (2006). Automating the collection of ultrathin serial sections for large volume tem reconstructions. Microsc. Microanal. 12(Suppl. 2), 86–87. doi: 10.1017/S1431927606066268

CrossRef Full Text | Google Scholar

Helmstaedter, M., Briggman, K. L., Turaga, S. C., Jain, V., Seung, H. S., and Denk, W. (2013). Connectomic reconstruction of the inner plexiform layer in the mouse retina. Nature 500, 168–174. doi: 10.1038/nature12346

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Jackson, J. I., Meyer, C. H., Nishimura, D. G., and Macovski, A. (1991). Selection of a convolution function for fourier inversion using gridding. IEEE Trans. Med. Imaging 10, 473–478. doi: 10.1109/42.97598

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Jain, V., Bollmann, B., Richardson, M., Berger, D., Helmstaedter, M., Briggman, K., et al. (2010). “Boundary learning by optimization with topological constraints,” in IEEE Conference on Computer Vision and Pattern Recognition (CVPR) (San Francisco, CA), 2488–2495.

Google Scholar

Ke, M. T., Fujimoto, S., and Imai, T. (2013). Seedb: a simple and morphology-preserving optical clearing agent for neuronal circuit reconstruction. Nat. Neurosci. 16, 1154–1161. doi: 10.1038/nn.3447

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Kim, I., Zhang, Y., Meister, M., and Sanes, J. (2010). Laminar restriction of retinal ganglion cell dendrites and axons: subtype-specific developmental patterns revealed with transgenic markers. J. Neurosci. 30, 1452–1462. doi: 10.1523/JNEUROSCI.4779-09.2010

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Kim, I., Zhang, Y., Yamagata, M., Meister, M., and Sanes, J. (2008). Molecular identification of a retinal cell type that responds to upward motion. Nature 452, 478–482. doi: 10.1038/nature06739

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Kim, J. S., Greene, M. J., Zlateski, A., Lee, K., Richardson, M., Turaga, S. C., et al. (2014). Space-time wiring specificity supports direction selectivity in the retina. Nature 509, 331–336. doi: 10.1038/nature13240

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Kong, J. H., Fish, D. R., Rockhill, R. L., and Masland, R. H. (2005). Diversity of ganglion cells in the mouse retina: unsupervised morphological classification and its limits. J. Comp. Neurol. 489, 293–310. doi: 10.1002/cne.20631

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

LeCun, Y., Bottou, L., Bengio, Y., and Haffner, P. (1998). Gradient-based learning applied to document recognition. Proc. IEEE 86, 2278–2324. doi: 10.1109/5.726791

CrossRef Full Text | Google Scholar

Levy, B., Petitjean, S., Ray, N., and Maillot, J. (2002). Least squares conformal maps for automatic texture atlas generation. ACM Trans. Graphic. 21, 362–371. doi: 10.1145/566654.566590

CrossRef Full Text | Google Scholar

Losavio, B. E., Liang, Y., Santamaria-Pang, A., Kakadiaris, I. A., Colbert, C. M., and Saggau, P. (2008). Live neuron morphology automatically reconstructed from multiphoton and confocal imaging data. J. Neurophysiol. 100, 2422–2429. doi: 10.1152/jn.90627.2008

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Narayanaswamy, A., Wang, Y., and Roysam, B. (2011). 3-D image pre-processing algorithms for improved automated tracing of neuronal arbors. Neuroinformatics 9, 219–231. doi: 10.1007/s12021-011-9116-z

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Peng, H., Long, F., and Myers, G. (2011). Automatic 3D neuron tracing using all-path pruning. Bioinformatics 27, i239–i247. doi: 10.1093/bioinformatics/btr237

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Peng, H., Ruan, Z., Atasoy, D., and Sternson, S. (2010). Automatic reconstruction of 3d neuron structures using a graph-augmented deformable model. Bioinformatics 26, i38–i46. doi: 10.1093/bioinformatics/btq212

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Schmitt, S., Evers, J. F., Duch, C., Scholz, M., and Obermayer, K. (2004). New methods for the computer-assisted 3-D reconstruction of neurons from confocal image stacks. Neuroimage 23, 1283–1298. doi: 10.1016/j.neuroimage.2004.06.047

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Srinivasan, R., Li, Q., Zhou, X., Lu, J., Lichtman, J., and Wong, S. T. (2010). Reconstruction of the neuromuscular junction connectome. Bioinformatics 26, i64–i70. doi: 10.1093/bioinformatics/btq179

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Stepanyants, A., and Chklovskii, D. B. (2005). Neurogeometry and potential synaptic connectivity. Trends Neurosci. 28, 387–394. doi: 10.1016/j.tins.2005.05.006

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Sümbül, U., Song, S., McCulloch, K., Becker, M., Lin, B., Sanes, J. R., et al. (2014). A genetic and computational approach to structurally classify neuronal types. Nat. Commun. 5:3512. doi: 10.1038/ncomms4512

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Sun, W., Li, X., and He, S. (2002). Large-scale morphological survey of rat retinal ganglion cells. Visual Neurosci. 19, 483–493. doi: 10.1017/S0952523802194107

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Turetken, E., Benmansour, F., and Fua, P. (2012). “Automated reconstruction of tree structures using path classifiers and mixed integer programming,” in IEEE Conference on Computer Vision and Pattern Recognition (CVPR) (Providence, RI), 566–573.

Google Scholar

Turetken, E., Gonzlez, G., Blum, C., and Fua, P. (2011). Automated reconstruction of dendritic and axonal trees by global optimization with geometric priors. Neuroinformatics 9, 279–302. doi: 10.1007/s12021-011-9122-1

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Völgyi, B., Chheda, S., and Bloomfield, S. (2009). Tracer coupling patterns of the ganglion cell subtypes in the mouse retina. J. Comp. Neurol. 512, 664–687. doi: 10.1002/cne.21912

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Wang, Y., Narayanaswamy, A., Tsai, C., and Roysam, B. (2011). A broadly applicable 3-D neuron tracing method based on open-curve snake. Neuroinformatics 9, 193–217. doi: 10.1007/s12021-011-9110-5

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Wickersham, I. R., Finke, S., Conzelmann, K. K., and Callaway, E. M. (2006). Retrograde neuronal tracing with a deletion-mutant rabies virus. Nat. Methods 4, 47–49. doi: 10.1038/nmeth999

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Wickersham, I. R., Lyon, D. C., Barnard, R. J., Mori, T., Finke, S., Conzelmann, K. K., et al. (2007). Monosynaptic restriction of transsynaptic tracing from single, genetically targeted neurons. Neuron 53, 639–647. doi: 10.1016/j.neuron.2007.01.033

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Xie, J., Zhao, T., Lee, T., Myers, E., and Peng, H. (2011). Anisotropic path searching for automatic neuron reconstruction. Med. Image Anal. 15, 680–689. doi: 10.1016/

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Zhang, Y., Zhou, X., Degterev, A., Lipinski, M., Adjeroh, D., Yuan, J., et al. (2007). Automated neurite extraction using dynamic programming for high-throughput screening of neuron-based assays. Neuroimage 35, 1502–1515. doi: 10.1016/j.neuroimage.2007.01.014

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Zhao, T., Xie, J., Amat, F., Clack, N., Ahammad, P., Peng, H., et al. (2011). Automated reconstruction of neuronal morphology based on local geometrical and global structural models. Neuroinformatics 9, 247–261. doi: 10.1007/s12021-011-9120-3

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Keywords: cell types, classification, retinal ganglion cells, reconstruction, stratification, laminar structures

Citation: Sümbül U, Zlateski A, Vishwanathan A, Masland RH and Seung HS (2014) Automated computation of arbor densities: a step toward identifying neuronal cell types. Front. Neuroanat. 8:139. doi: 10.3389/fnana.2014.00139

Received: 17 July 2014; Accepted: 06 November 2014;
Published online: 25 November 2014.

Edited by:

Hermann Cuntz, Ernst Strüngmann Institute in Cooperation with Max Planck Society, Germany

Reviewed by:

Karl Farrow, Neuroelectronics Research Flanders, Belgium
Marta Costa, University of Cambridge, UK

Copyright © 2014 Sümbül, Zlateski, Vishwanathan, Masland and Seung. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

*Correspondence: Uygar Sümbül, Grossman Center for the Statistics of Mind and Department of Statistics, Columbia University, New York, NY 10027, USA e-mail:

Present address: Uygar Sümbül, Grossman Center for the Statistics of Mind and Department of Statistics, Columbia University, New York, NY, USA

These authors have contributed equally to this work.