The accurate tracing of single neurons is one prerequisite for the determination of anatomical features of different neuronal cell types required for biophysical modeling of single cells and signal processing in small circuits. Several automated reconstruction approaches have been reported previously. These reconstructions, made with two-photon,1 confocal 2, 3, 4, 5, 6, 7, 8, 9, 10 or brightfield images,10, 11 focus mainly on the dendritic tree and the extraction of geometrical features such as volumes or surface areas of dendrites or spines. Hence, these approaches completely lack a reconstruction of the axonal arbor. The reason is that the axonal branching pattern is more complex and that axons spread over a much larger volume (cubic millimeters) compared to dendrites (a few hundred cubic micrometers). Furthermore, axonal staining is fainter than that of dendrites because of the smaller diameters of axons and their greater distance from the soma where a tracer is loaded into the cell. Thus, no successful automated tracing of these widely spreading neuronal projections has yet been reported.
Nevertheless, the determination of individual cells’ projection patterns,12 of connections between different areas of the brain, the number of synaptic contacts between cells and the difference in axonal growth under various developmental conditions12, 13, 14, 15 are problems that illustrate the need for detailed three-dimensional (3D) reconstruction of axonal arbors. The method presented here focuses on the accurate tracing of all neuronal projections. The dendritic tree and the exten-sively spreading axonal arbor are traced and reconstructed simultaneously.
Typically the cells are filled in vivo with a tracer, like biocytin.16 The brain is then perfused with fixative and cut in sections of about thickness. The current approach for tracing these neurons and their axons in three dimensions is a manual one, based on the Camera Lucida technique (e.g., The Neurolucida System,17 FilamentTracer.18 Here, neuronal structures in each section are traced manually. A human user interacts with a microscope that is enhanced with computer imaging hardware and software.11, 13 The user recognizes patterns and marks neuronal structures on a live camera image displayed on a computer screen. By moving the stage in the and directions and focusing through the brain slice in the direction, a progressively larger volume is inspected. The 3D tracings of neuronal branches from this volume are collected by the computer system interfaced to the camera and result in a 3D graph representation of neurons.
Manual tracings of dendritic trees are very reliable. The reliability is due to the localized branching pattern and the relatively large diameters of dendrites. However, the axonal arbor frequently extends further away from the soma. The average volume that has to be inspected by the user is for most cell types around per brain section. For accurate tracing, this volume should be inspected in a raster scan order, moving from one field of view to the next and progressively focusing through the specimen.
Using a objective and a standard charge-coupled device (CCD) camera [e.g., Q ICAM (Q Imaging, Surrey, British Columbia, Canada)], the number of fields of view is approximately 5 500, with an average sampling along the optical axis of . Because axons can have diameters less than , a objective is usually used for manual tracing. In this case, approximately 37 000 fields of view have to be inspected. Taking a typical cell (e.g., layer pyramidal neuron of rat cortex) that spreads over 10 to 20 brain sections, the number of fields of view that have to be inspected is of the order of . The manual reconstruction of axonal branching patterns is hence tedious and time-consuming. Therefore, correct manual tracing of axons requires experienced users to reach a reliable level of reconstruction quality.
We present a semiautomatic “reconstruction” pipeline (Fig. 1 ) that traces reliably both dendrites and axonal arbors by extracting their skeleton (approximate midlines). In addition, the method presented here needs significantly less time compared to manual tracing. The automated tracing is carried out on large image stacks acquired by mosaic scanning19 and optical sectioning20 using a transmitted light brightfield (TLB) microscope. The manual inspection of thousands of fields of view is replaced by the acquisition of a stack of mosaic images. A rectangular pattern of overlapping mosaic tiles (adjacent fields of view) is scanned, covering an area of a brain slice sufficient for the tracing of axonal arbors usually [Fig. 2a ]. The 3D information is obtained via optical sectioning, meaning the recording of such mosaic planes at multiple focal positions that are ideally separated by .
This large 3D image is then deconvolved21 based on measured optical properties of the microscope, such as lack of aberrations within the optical pathway. This guarantees an improvement of signal-to-noise ratio and resolution, in particular along the direction. A local threshold function that checks for connectivity extracts neuronal structures from the background. The extracted foreground objects are transformed into thinned approximate midlines (the skeleton) and yield a graph representation of dendrites and the widely spreading axonal arbor. The automatically reconstructed branches from serially sectioned brain slices are than manually edited and spliced using NEUROLUCIDA software for editing.
All cells were filled in Wistar rats (P28 to P31; Charles River Laboratory). Experiments were carried out in accordance with the animal welfare guidelines of the Max Planck Society.
Briefly, animals were anesthetized with urethane or isoflurane. Body temperature was maintained at by a servo-controlled heating blanket. A metal post for positioning the head was attached to the skull overlying the cerebellum by dental acrylic. A craniotomy was made overlying the left barrel cortex (each ) or the ventral posteriomedial nucleus of the thalamus . Thalamic craniotomies were centered at posterior to bregma and lateral of the midline. The dura was removed for some of the experiments. The craniotomy was kept covered with artificial cerebral spinal fluid [aCSF; in millimoles: 135 NaCl, 5.4 KCl, 1.8 , 1.0 , and 50 -(2-hydroxyethyl)piperazine- -2-ethanesulfonic acid (HEPEs); pH 7.2].
Cells were filled with biocytin either extracellularly using juxtasomal recording and electroporation22 or via whole-cell recording.23 In both cases, we used patch pipettes, pulled from unfilamented borosilicate glass on a Sutter P-97 puller (Sutter Instruments). The outside diameter of the shank entering the brain varied from , and the tip opening had an inside diameter less than .
Whole-cell recording pipettes were tip-filled with (in millimoles) 135 K-gluconate, 10 HEPEs, 10 phosphocreatin-Na, 4 KCl, 4 adenosine triphosphate–Mg, 0.3 guanosine triphosphate, and 0.2% biocytin (pH 7.2, osmolarity 291). Pipettes were inserted perpendicular to the pia under high pressure , and cells were searched for blindly under low pressure .23 Recordings were made in bridge mode for during which biocytin passively diffused into the cells. Seal resistance was , access resistance was , and spike height and overall were stable throughout the recording. No holding current was used.
Juxtasomal fillings were made with pipettes filled with aCSF containing 1% biocytin. No pressure was applied to the pipette at any time. Once a cell was isolated extracellularly, current pulses ( , on, off) were applied continuously for several minutes.22
At least was allowed to pass after cell filling and prior to tissue fixation to ensure sufficient diffusion of the biocytin throughout the axonal arbor. The animal was deeply anesthetized and perfused transcardially with phosphate buffer followed by 4% paraformaldehyde. Cortex and thalamus were cut tangentially and coronally, respectively, in vibratome sections. Biocytin in these sections was stained with the chromogen -diaminobenzidine tetrahydrochloride using the avidin-biotin-peroxidase method.16 Sections were sometimes counterstained for cytochrome oxidase.24 Processed sections were then mounted on slides and coverslipped with Mowiol (Hoechst, Austria). For comparison with our reconstruction approach, four well-filled cells from four different rats were manually reconstructed using a Neurolucida system with a oil immersion objective [Olympus Plan; 1.25 numeric aperture (NA).]
A standard TLB microscope (Olympus BX-51, Olympus, Japan) equipped with a motorized stage (Maerzhaeuser Wetzlar, Germany) was used for image acquisition [1 in Fig. 1a]. A bandpass illumination filter (CHROMA AF-analysentechnik, Germany), which is attached to the diaphragm of the lighthouse, provides quasi-monochromatic illumination of the specimen [2 in Fig. 1a]. This bandpass filter minimizes the chromatic aberrations of the imaging system and simplifies the theoretical description of the optical pathway from a polychromatic to a monochromatic one.
The light is transmitted by a high NA (1.4 NA) oil immersion condenser (Olympus, Japan) [3 in Fig. 1a], ensuring parallel illumination of the specimen under “Koehler”-conditions.25 The specimen is imaged by a high NA oil immersion objective (Olympus UPLSAPO; 1.4 NA) in combination with a television (TV)-mount (Olympus U-TV0.5XC-3) or a objective (Olympus UPLFLN; 1.3 NA) in combination with a nonmagnifying TV mount [4 in Fig. 1a]. The immersion oil has a refractive index of similar to glass.
The stage is navigated in three spatial directions by OASIS-4i-controller hardware and software (Objective Imaging Ltd., Cambridge, United Kingdom). It allows the acquisition of large mosaic images19 at different focal planes. Mosaic in this context refers to a two-dimensional (2D) image of overlapping tiles (i.e., adjacent fields of view) that are aligned automatically and then stitched during the image acquisition, resulting in a large composite image [Fig. 2a]. The user-defined scan area is automatically divided into a series of overlapping fields of view that we call tiles. For each tile location, a stack of images is acquired using optical sectioning20 with a typical separation of between the focal planes. For each focal plane, the corresponding tiles are then stitched together, resulting in a 3D stack of 2D mosaic images for each focal position [Fig. 2b]. This process is executed by the Surveyor Software (Objective Imaging Ltd.).
The images are recorded by a Q ICAM Fast 1394 camera (QImaging) [5 in Fig. 1a] equipped with a CCD chip, which in combination with the objective and the 0.5 TV mount yields an sampling of ( per pixel if the objective and a nonmagnifying TV mount are used). Because the illumination is limited to one wavelength, eight-bit gray value images, rather than red-green-blue color images, are acquired. Other cameras, such as the Retiga 2000R Fast 1394 Mono Cooled (QImaging), were tested as well and showed no significant improvement in resolution or signal-to-noise ratio.
To guarantee a similar dynamic range of the gray values for mosaic image stacks from different brain slices and animals, the exposure time of the CCD camera is set semiautomatically by the Surveyor Software. Therefore, the mean gray value in a typical field of view within the scan area is calculated and is set to be 190. Typical here refers to a field of view without any stained neuron somata or blood vessels. This highly important feature of the image acquisition is the basis for an optimal deconvolution and, therefore, a robust neuron tracing. However, significant exposure gradients within the scan area are a considerable constraint to our method. The consequences and limitations for neuron tracing will be discussed later.
In summary, image acquisition results in a high-resolution 3D mosaic image. The typical image stack for neuron tracing that we use is . However, the mosaic pattern can be adjusted to smaller or larger areas. Using the objective in combination with the 0.5 TV mount or the objective with a nonmagnifying TV mount, the sampling is or per voxel. Hence, the data volume for such a stack is approximately .
Hardware and Software Requirements for Automated Image Processing
In general, the data size is between 10 and per section. Therefore, a fast image processing pipeline based on multi processor computing that can handle such large data sets was developed. The algorithms are written in .26 The raster image file input/output, iteration through a raster image and the dilation as well as the closing filter, use the ITK Image Processing Library.27 The algorithms were parallelized by applying the OPENMP standard.28 They are executed on AMD dual-core Opteron servers, equipped with either 4 central processing units (CPUs) and memory (DELTA Computer Products GmbH, Reinbek, Germany) or 8 CPUs and memory (fms-computer.com, Netphen, Germany).
Once the mosaic image stack has been saved to disk, a software daemon (a custom written PERL script) detects new data on the hard drive and starts the processing pipeline. Additional image stacks are added to a queue and processed sequentially. The script initiates the subsequent image processing steps: inversion of the gray values, subdivision into bricks, deconvolution, down sampling, and finally the 3D neuron tracing [Figs. 1b, 1c, 5]. Therefore, the automatic tracing routine can run a day, because the queue of scanned image stacks is sequentially processed. The next step needing user interaction is the manual postprocessing of the final 3D reconstructions.
Diffraction as well as interference phenomena decrease the image quality and resolution of image stacks acquired through a TLB microscope. The resolution along the optical axis ( axis) is usually insufficient for an automatic tracing of faint axonal branches surrounded by more prominent neuronal structures.11 Image restoration (deconvolution) can partly compensate for this limited resolution. Deconvolution is the most important step for the performance of the automatic tracing procedure.29
The image formation is described by a convolution of the 3D intensity distribution within the object with a point spread function (PSF).30 The PSF (Fig. 3 , see Color Plate 1) describes the optical properties of the imaging system. The convolution results in the 3D intensity distribution within the image. The intensity distribution within the object is then restored by deconvolving the recorded image stack using the PSF. Deconvolution is essentially an inversion of the image formation process. Thus deconvolution yields the best results provided the PSF can be determined exactly.
Inversion of gray values
The TLB microscope is treated like a fluorescence microscope that lacks primary aberrations. Image formation can then be simplified and modeled as a monochromatic point source that is imaged by the circular entrance pupil of the objective. This can be calculated analytically for low NA objectives25 and numerically for high NA ones.31 Justification for these simplifications, as well as verification, can be found in Ref. 21. Simplifications are essentially based on considerations of propagation of mutual coherence32 and measurements of spherical aberrations, using a Shack-Hartmann wavefront sensor.33 Treating the TLB images like fluorescence data, the gray values have to be inverted. This inverted image stack will be referred to as the “original image” stack.
Subdivision into bricks and linear deconvolution
The deconvolution is carried out by the HUYGENS software.34 With the assumptions described above, the PSF is calculated in a wide-field fluorescence mode (Fig. 3), which is suggested by the HUYGENS software. An important parameter is the refractive index of the specimen. Aberration mea-surements21 yield a uniform value of for the specimen that is a -thick brain slice embedded in Mowiol . Because the simplifications of the image formation process are justified, the imaging system is well described by the modeled PSF (Fig. 3).
This PSF is applied to the recorded image stack by a linear Tikhonov-Miller restoration filter.35 It is derived from a least squares approach. This approach is based on minimizing the squared difference of the acquired image and a blurred estimate of the original object. The resultant image stack is improved significantly in terms of resolution and signal-to-noise ratio (Fig. 4 ). This deconvolved imaged stack is then of sufficient quality to apply an automatic tracing of neuronal structures (Figs. 5, 6, 7 ).
Using the Tikhonov-Miller filter, the HUYGENS software is not capable of processing images larger than at once. Therefore, the mosaic image stack is divided into smaller stacks of appropriate size (e.g., voxels), overlapping by 100 voxels, that can be processed by HUYGENS. These substacks, hereafter called bricks, are regrouped during the neuron tracing, resulting in a deconvolved and segmented mosaic image stack.
Maximum down sampling
After deconvolution, the data size is reduced to facilitate further computationally intensive steps. Maximum down sampling is applied for each deconvolved image brick in the plane. Two neighboring voxels are merged to one voxel, computing their maximum gray value. Because this procedure is done in both lateral directions ( and ), the data set is reduced in size by a factor of 4. Because neuronal structures are local intensity maxima, all branches will be conserved.
As a consequence, the new lateral sampling is per voxel if the objective in combination with the 0.5 TV mount is used or per voxel if the objective in combination with a nonmagnifying TV mount is used. Therefore the sampling in both cases is still below the physical resolution limit (approximately half the wavelength, ) according to the Rayleigh criterion.25
Raster image–based segmentation
A major feature of image acquisition is the semiautomated setting of the exposure time of the CCD camera to control the image histogram as described previously. This feature is critical to the downstream neuron tracing algorithms. Their parameters were adjusted and systematically tested to derive the best possible tracings from such images. This issue of neuron tracing will be discussed later.
The down-sampled 3D image bricks are now subjected individually to segmentation algorithms (Figs. 5, 6, 7). The purpose is to separate voxels that represent neuronal structures from voxels that are part of the background.
Local threshold filtering
A target image with the dimensions of a deconvolved image brick is created and initialized with gray value zero. The voxels of the deconvolved bricks are then separated into three groups.
One group, below a lower threshold is set to background (gray value 0 in the deconvolved brick). A second group of voxels belonging to potential neuronal structures, with values above an upper threshold, is assigned as foreground (gray value 255 in the target image). The foreground consists of disjointed voxel regions that will be referred to as foreground objects. Voxels with values between the two thresholds form the third group of intermediate value voxels. The intermediate voxels are tested for local features defined by a local property function to decide whether they belong to foreground objects or the background.
The lower threshold is determined by calculating the intensity distribution of the deconvolved image brick. The deconvolution produces a thin unimodal histogram with the background clustered near 0, structures in the high range, and some remaining intermediate gray values. The lower threshold is taken as the histogram’s mean value plus 1.5 standard deviations (SDs). This gray value is usually between 1 and 15. It was derived after systematic testing and essentially deletes the background noise from unstained tissue.
The upper threshold assigns the prominent and most intense structures to the foreground. The value of the upper threshold is determined by calculating the histogram of the maximum projection of the deconvolved image brick. The upper threshold is then taken as the mean value plus 3.0 SDs. Systematic testing yielded that this is the best value to detect dendrites, well-filled axons, and the prominent parts of fragmented filled axons [Figs. 5c, 6c, 7c]. This upper threshold usually has a gray value between 30 and 60. The remaining voxels between the background (1 to 15) and the foreground (30 to 60) margins are referred to as the intermediate voxels.
The local property function comprises two steps. First, each intermediate voxel is set as the center of an -voxel mask in the plane and the mean intensity of the voxels, regardless of group (background voxels were already set to 0), within this 2D neighborhood is calculated. If the centered intermediate voxel has a gray value that is larger than this mean intensity plus an epsilon value of 5 gray values, it is set to an intermediate gray level of 125 [Figs. 5c, 6c, 7c] in the target image. Otherwise it is set to background (gray value 0 in the target image). The optimal value of was derived by systematic testing. The group of intermediate voxels consists usually of isolated artifacts or dim bridges between bright structures. These dim bridges are often found between axonal boutons (swellings of the axon that are likely sites of synaptic contact [Figs. 6(a2) to 6(e2)] and should therefore be part of the foreground. The resultant target image comprises three gray values: 0, 125, and 255 for the background, intermediate and foreground voxels, respectively.
The second part of the local property function inspects the intermediate voxels of the target image for connectivity to a foreground object. Systematic testing suggested that an intermediate voxel is set to the foreground if 10% of the voxels from a mask in the plane centered on this intermediate voxel are part of a foreground object, otherwise it is set to background [Figs. 5d, 6d, 7d].
Hit or miss transformation
The application of the above local threshold filter results in a binary image. This is then subjected to a hit or miss transformation36 with rectangular frame masks of increasing size as structuring elements. The transformation is applied to every image plane. Isolated foreground objects that are completely surrounded by a frame are converted to background. Beginning with a radius of one voxel and increasing the frame size subsequently to three voxels, small, isolated artifacts that were introduced by the linear deconvolution are removed [Figs. 5e, 6e, 7e].
Dilation and closing
Next, dilation and closing filters36 are applied. The dilation filter bridges gaps between the axonal boutons that have not been closed by the local threshold filter. Finally, a closing filter is applied [Figs. 5f, 6f, 7f]. Its geometrical interpretation is that a “sphere” rolls along the outside boundary of a foreground object within the image. Therefore, it tends to smooth sections of contours and fuses narrow breaks as well as long thin gulfs, eliminates small holes, and fills small gaps in the contour.36 The 3D structuring element (sphere) for the closing and dilation has a radius of three voxels. Larger radii could result in a fusion of objects that should not be connected, and smaller radii would have no significant effect on the foreground objects. The radius of three, therefore, proved to be the best compromise after systematic testing.
The brickwise segmentation yields binary image stacks. Regrouping of these segmented bricks leads to a binary 3D image stack with the dimensions of a down-sampled mosaic image stack. This large 3D binary stack is subjected to an object labeling algorithm. Once a foreground voxel is detected, a region growing algorithm36 fuses all connected foreground voxels to a subregion and assigns a unique integer label. The binary image is thus transformed to an image of individually labeled subregions, each representing a foreground object. Labeled objects are sorted according to their number of voxels. The largest foreground object is labeled as number 1, and the smallest object is assigned number . The background voxels are labeled as 0.
Raster to vector image conversion
Because the neuron occupies only a small fraction of the scanned image volume, the 3D raster image representation requires inappropriately large data storage. A more sophisticated vector data folder is created that stores only foreground objects but keeps the 3D information of the image. Figure 8 illustrates schematically the architecture of the new data folder. The voxels of each foreground object are replaced by vectors, hereafter called compartments. Each compartment stores the 3D coordinates of the prior voxel and could in the future store additional morphological information, such as local properties (e.g., radii, surface distance).
Furthermore, each compartment is linked to its neighboring compartments by intercompartmental smart pointers, which allow access of neighboring compartments (e.g., top-left, bottom, left). Therefore, the implicit neighborhood representation by voxel coordinates of the 3D raster image is mapped to an explicit neighborhood construct in the compartment representation. The neighborhood pointers between the compartments preserve the 3D topology of any object. Furthermore, fast navigation through the object is possible using these intercompartmental smart pointers.
The compartments representing one object are grouped in a “double linked” list, which guarantees fast access to each compartment. Therefore, the 3D raster image of labeled foreground objects is converted to compartment lists. This transformation replaces the 3D raster image of labeled foreground objects by a double linked list, where each list item represents one of these disjointed foreground objects by a list of topology preserving compartments (vector image representation).
Vector image–based midline extraction
Representing objects by their main structure or skeleton (approximate midline) is a commonly used technique in image processing. A fast and reliable way to calculate the skeleton of an object is thinning. Generally, thinning is a layer-by-layer (boundary voxel/compartment) erosion until only a unit-width skeleton is left.37, 38, 39 Therefore, three different classes of compartment topologies are defined: 6-adjacent 18-adjacent , and 26-adjacent .38 Two neighboring compartments are if their Euclidean distance is equal to 1, if their Euclidean distance is between (or equal to) 1 and the square root of 2, and if their Euclidean distance is between (or equal to) 1 and the square root of 339 (Fig. 9 , see Color Plate 1).
A thinning algorithm has to obey the following demands that are derived from a 2D definition by Seul :40 (1) Connected objects must thin to connected line structures. (2) The thinned lines should be minimally 26-connected. Two compartments are 26-connected if they are connected by a chain of adjacent compartments. (3) Approximate end-line locations should be maintained. (4) The result of thinning should approximate the midlines of the structures. (5) Extraneous short branches introduced by thinning should be minimized.
To address these demands, we use the template matching algorithm described by Jonker.41 It is one of the most efficient thinning algorithms and can be briefly summarized as follows. First, a set of end-line locations is determined to represent the topology of the object. The end-line locations will not be removed. All other compartments, starting from the object boundary, are tested to check whether or not their removal affects the 26-connectivity of the object (Fig. 9). If compartment removal has no affect on the connectivity, this compartment is termed simple.
To determine whether a compartment is simple or not, its neighborhood is compared with all possible neighborhood templates that preserve local connectivity (Fig. 9). If the compartment neighborhood matches one of these templates, the removal of this compartment will not affect the local connectivity and the global connectivity will be preserved as well.41 Hence, this compartment is assigned to be simple and is removed. Whether or not a compartment is simple may change after deleting compartments. Therefore, the algorithm is iterative and finishes when no more compartments are removed (i.e., the thinning algorithm is an idempotent algorithm. The implementation of this thinning approach is described in detail below.
Detection of end-line locations
The end-line locations are determined by calculating the compartment, which is most distant to the object boundary (compartments with less than 26 neighbors). This is done by calculating the Euclidean distance map42 of the object. The compartment with the highest distance number is selected and called the seed point. In the case of several compartments sharing the highest distance value, one of them is selected randomly. Now the Euclidean graph distance39 from each compartment in the object to the seed compartment is computed. This is the shortest connection along the graph between two compartments. Compartments with a local maximal graph distance are assigned as line endings [Figs. 5g, 6g, 7g].
One disadvantage of thinning is the possibility of inward erosion, which could potentially create holes in the object. To prevent these artifacts, an iterative approach is used.38
Using the 6-adjacent and 18-adjacent topologies defined above, boundary layers are assigned that are peeled by an iterative algorithm:
1. Collect the set of boundary compartments (compartments where at least one of the neighbors is missing).
2. Delete all simple compartments (according to the topology) in this set, starting with those who have the lowest number of neighbors (i.e., intercompartmental pointers).
3. Collect the set of boundary compartments (compartments where at least one of the neighbors is missing).
4. Delete all simple compartments (according to the topology) in this set, starting with those who have the lowest number of neighbors.
Graph validation and pruning
Imaging and segmentation can lead to artifacts such as loops and clusters in the thinned midline representation of the objects, which have to be removed by postprocessing. Loops are removed by selecting one end-line location and calculating the shortest path from this end-line location to all other end-line locations in the object. Compartments or intercompartmental pointers that are not used by any path within this validation step [Figs. 5i, 5j, 6i, 6j, 7i, 7j] are removed from the object.
Because the approximate end-line locations were local distance maxima from an inner seed point, swellings within the neuronal branches will result in a short subbranch that we regard as artificial [Figs. 5k, 6k, 7k]. Therefore, these short branches are removed within this pruning step from the object if the distance from the according end-line location to the first intersection (Fig. 9) compartment is shorter than (derived by systematic testing). Nonartificial axonal subbranches below this length threshold will hence be pruned as well. However, for most scientific problems, this error is negligible.
Editing and Splicing
The image processing, so far, transformed the deconvolved raster image into a vector image representation of thinned foreground objects. Therefore, each object is represented as a 3D graph with end-line locations (endings). Furthermore, compartments that are connected to more than two neighbors will be set to be an intersection (Fig. 9). This set of endings, intersections, and normal compartments is converted to an ASCII (American Standard Code for Information Interchange) file that can be either visualized in AMIRA 43 or NEUROLUCIDA. Up to this point, all processing steps are automated. Final editing and splicing of the reconstructions is done manually.
Artifacts that are similar in shape or gray value to neuronal structure have to be erased. In most cases, these artifacts are astrocytes (star-shaped glial cells). This editing is done by superimposing the reconstruction with the maximum projection of the original image stack. Within the projection image, neuronal structure can be easily distinguished from artifacts. Neuronal branches appear as elongated structures within the projection image; whereas, falsely traced artifacts are either small spots (most intense parts of unstained neurons) or star-shaped if they are the remains of glial cells. These artifacts are then marked and erased in NEUROLUCIDA, which is used here as a graph editor. Axonal branches that were not well-filled and, therefore, traced as a fragmented structure can be restored by splicing their closest ending points to form a single continuous axon arbor.
Within this step, additional anatomical landmarks, such as section outlines, can be integrated. Superimposing the edited neuron reconstruction with a low magnification image, any visible anatomical structure (e.g., pia outlines) can be manually added to the automatic reconstruction of the neuron.
The results of this semiautomated reconstruction pipeline are 3D tracings of neuronal branches from adjacent, -thick sections. To obtain a complete reconstruction of a cell, the reconstructions from different sections have to be aligned and spliced (Fig. 10 ). Fortunately, large radial blood vessels running perpendicular to the pial surface retain their positions within each tangentially cut brain slice, rendering them ideal as position reference points. Therefore, the blood vessel pattern is extracted during the automatic tracing procedure. Blood vessels are detected via a region-growing algorithm in a maximum intensity projection of the original image stack, prior to the deconvolution. Voxels with values close to 0 (maximum transparent regions) are used as seed points for a region-growing algorithm. The region-growing algorithm seeks dark and almost circular shapes (shape number algorithm)36 and labels these structures as blood vessels. According to the vessel pattern, the reconstructions are aligned coarsely. More accurate alignment is then done by matching the upper branch endings from a lower slice with the lower branch endings from an upper slice. Once all branches are aligned with their counterparts from the adjacent section, the nearest ending points are spliced (Fig. 10). Starting at the slice including the soma, progressively more sections will be connected resulting in a complete 3D cell reconstruction [Fig. 11a ]. If no soma is present (e.g., the cortical part of the thalamocortical axonal arbor), the splicing starts at the deepest (most distant from the pia surface) section [Fig. 11b].
An international patent application is pending.44
Fractal Box Counting for Quantification of Relative Shape Deviation
For comparing the shape of automatically reconstructed branches with manual Camera Lucida–based reconstructions, a relative measure was used.45 The aim was to obtain the percentage of relative deviation in shape at various resolutions. The shape comparison was applied to 2D projections of the traced branches. This reduction to projections is necessary because the embedding medium (Mowiol) may expand or contract over time depending on the ambient room temperature. After several hours at room temperature, the tissue tends to shrink in the direction (along the optical axis). After cooling the slide, this shrinkage is reversed (Table 1 ). Thus a comparison of the 3D morphology could be confounded by shrinkage and expansion of the tissue.
First, a bounding box is put around the projection of an individual branch [Fig. 12a ]. Second, this box is divided into eight subboxes, and the numbers of boxes that are intersected by the branch are counted [Fig. 12b]. In subsequent steps, each box is divided into eight boxes45 until a box size between 200 and is reached. This is done for both the manually and automatically traced branches.
If two reconstructions would match perfectly in shape, the number of boxes intersected by the branch would be the same within each dividing step. If one branch differs in shape from its counterpart, the number of intersected boxes will be different from a certain box size onward. We checked the difference in the number of intersected boxes at a box size (resolution) of , , , and (Fig. 13 ), which is nearing the physical resolution limit of the imaging system (Rayleigh Criterion).25 The number of intersected boxes from the automatically reconstructed branch at the highest resolution is set to be 100% of the branch. Differences in the number of intersected boxes will be given as percentages with respect to this value. This is illustrated in Fig. 14 , where the differences in structure-intersected boxes are shown for an increasingly cropped branch [Fig. 14(a1) to 14(a1d)].
An important constraint for the image processing is based on the exposure time setting of the CCD camera before the image acquisition process. The mean gray value in a typical field of view within the scan area is calculated and set to be 190. Therefore, the dynamic gray value range for neuronal structures is maximized and comparable within image stacks acquired from different brain slices. Furthermore, we avoid clipping of gray values,34 which prevents artifacts from being introduced by deconvolution.
After acquisition of the mosaic image, the gray values are inverted and the large 3D mosaic image is subdivided into 3D bricks that can be processed by the HUYGENS software in combination with a linear Tikhonov-Miller filter. Because the deconvolution is performed on individual bricks, the gray value ranges for neuronal structure can be slightly shifted between the bricks. This is because HUYGENS estimates the background for each individual brick. The brickwise local threshold filtering compensates for this effect. Hence, the regrouped segmented mosaic image shows a uniform (not brickwise) distribution of foreground objects.
Connected subregions of foreground voxels within the segmented raster image are converted to lists of compartments with explicit intercompartmental neighborhood connections (vector image). The background is usually to times larger than the sum of the foreground voxels. Therefore, the storage size and the processing time decrease significantly after the data conversion from raster to vector image. Instead of processing gigabytes of data, the data set is reduced to a few hundred megabytes.
Template-based thinning reduces each object to its skeleton (approximate midline) [Figs. 5h, 6h, 7h]. The resulting graph is further validated and corrected until all thinning artifacts are erased. This correction usually creates smooth midline representations [Figs. 5l, 6l, 7l], neglecting morphological swellings such as dendritic spines or axonal boutons.
The tracing of poorly stained axons can result in a fragmented reconstruction of these structures (Figs. 5, 6, 7). Further small artifacts from unstained tissue, mainly astrocytes, are traced as well. Splicing of fragmented structures and erasing of artifacts during the manual postprocessing result in a 3D reconstruction of all neuronal branches within the brain section (Fig. 1).
Comparison to State-of-the-Art Technique
The state-of-the-art technique for the 3D reconstruction of neuron morphology, especially of the largely extending axonal arbors, is the computer-aided interactive Camera Lucida technique (e.g., Neurolucida).13 Therefore, we compared the performance of the semiautomated reconstruction method to the manual reconstruction approach, as implemented with the Neurolucida system.
First, we checked whether manually traced branches were detected by the automated system and vice versa. We quantified whether branches were missed by either of the reconstruction approaches. In a second comparison step, the relative deviation in shape between manually and automatically traced neuronal structures was determined. Finally, the time frame for reconstructing complete cells with their entire axonal arbor was compared for the conventional method systems. The comparison was done for 54 slices containing 4 different filled cells (from 4 different rats). They were randomly chosen from a pool of already manually reconstructed cells. Each manual reconstruction was done by a different individual. Users manually traced a pyramidal neuron in layers 2 and 3, two star pyramids from layer 4, and one ventroposteromedial (VPM) thalamocortical axon.
Number of identical branches
Figure 15 shows a comparison between an automated reconstruction of a section containing parts of a thalamocortical axon and its manually reconstructed counterpart. Example slices for the 4 neurons can be found in Figure 16 .
Each automatically reconstructed cell was inspected for missing braches that have been traced by the manual system and vice versa. The amount of manually reconstructed branches that were missed by the automatic system was similar for all four traced neurons. For the thalamocortical VPM cell, 5.4% were missing (Fig. 17 ), 3.5% and 3.4% for the layer 4 cells, and no manually reconstructed branch was missing for the cell in layers 2 and 3. In the case of branches that were missing in the manual tracing but found by the automated system, the result was less homogeneous. Approximately 4% of automatically reconstructed branches were missing in the manual version of the thalamocortical VPM cell, 29% and 47% were missed for the layer 4 cells, and almost 70% of the automatically traced branches were not found manually for the neuron in layers 2 and 3.
For the 4 cells, the automated method missed 29 of 682 manually reconstructed branches in 54 sections from 4 different rats. In other words, the new reconstruction technique found 96% of the manually reconstructed branches. The manual tracings were missing 268 of 923 automatically reconstructed branches from the same 54 slices. This means only 71% of the neuronal structures traced by the automated method could be found in the manual reconstructions. All of the missing branches, independent of whether they were missed by the manual or automatic system, were part of the axonal arbor. The dendritic branches were traced equally reliably by the two methods.
Fractal box counting
The fractal box counting method was applied to 167 randomly picked branches from the 54 slices that were reconstructed by both methods. The frequency of the same percentage in the difference of the number of intersected boxes between an automatically traced branch and its manual counterpart is shown as a histogram (Fig. 18 , see Color Plate 1). For each of the four box sizes, a histogram of intersected boxes was made and fitted with a Gaussian curve. The mean values of these Gaussian fits yielded the average deviation in shape at various resolutions. The relative deviation in shape is almost 0% at the lowest resolution, meaning a box size of . With increasing resolution (decreasing box size), the deviation increases to a value of around 2% at a box size of and further to more than 3% at a resolution of . At the highest measured resolution of , the relative shape difference was around 5%. However, it can be seen that some branches deviate a lot more in shape than the mean values suggest. These large deviations were mainly observed for branches from dendritic trees.
The average amount of time for a manual tracing of 1 cell (10 to 30 sections), including its axonal arbor, is approximately 60 to 90 working hours. This means the reconstruction of a single brain slice takes between 4 and .
Using the automatic system, an image size of proved to be sufficient for most slices. The computing time for an image stack of this size varies from , depending on the amount of structure within the volume. The human interaction time with the system is approximately per slice. This comprises the setting up of the scanning pattern , as well as the manual editing and splicing of the automatically reconstructed slices . Therefore, the average amount of time for reconstructing 1 cell with the semiautomatic system is between of computing and of human interaction. Given a working day of , the reconstruction using the manual system takes approximately , whereas the automatic system yields a comparable reconstruction within (Fig. 17).
The method of semiautomated reconstruction is based on the acquisition of stacks of rectangular mosaic images. One of the key steps is the high-resolution deconvolution, improving the image quality of stained neurons significantly. This is sufficient to extract neuronal structures from the background by applying a local threshold filter that seeks connectivity. The segmented raster image, consisting of disjointed foreground subregions is then converted to a vector image of objects. These are converted to thinned midlines, yielding a graph representation of the neuronal structure candidates. These graphs are then manually edited and spliced, resulting in a complete 3D reconstruction of the neuron.
The results showed that approximately 96% of the manually reconstructed axon branches were traced by the automated system. For each of the 4 cells, the performance of the new reconstruction technique is of a similar quality in terms of tracing branches that were found in their manual counterparts. The 4% of structure details that were not recovered fall into three groups.
The first group ( in the sample of four cells) consists of faint axonal branches. An indicator for faintness is that these branches are not local intensity maxima. Therefore, the deconvolution algorithm tends to weaken these branches instead of amplifying their intensity.
The second group (less than 1%) comprises small axonal endings that are not distinguishable from reconstructed artifacts and, therefore, are erased by the user during manual editing and splicing. To minimize the loss of small axonal endings, we superimpose the maximum projection of the original image stack with the 3D reconstruction.
The third group comprises branches that are not located within the scanning area. By default, an area of is scanned. For some sections, this area is not sufficient as the axons spread over a larger area. For these particular brain sections, the missing parts of the slice can be scanned. The branches missed initially are simply added to the rest of reconstructed branches, using NEUROLUCIDA for editing.
Summarizing these three groups, one can state that the 0 to 5% of missing automatically reconstructed branches consist either of small endings, faint branches, or not-scanned branches. By additionally scanning adjacent areas and inspecting the projection image, the number of missed branches can be minimized. Therefore, the performance of the semiautomatic reconstruction pipeline is limited by the first group, which implies the prerequisite that neuronal structures are local intensity maxima.
Comparison with Manual Tracing
Only 71% of the automatically traced branches were found by the manual tracers. The number of missing structures varied widely between 4% and almost 70%. The missed branches can be grouped into two cases.
The first group consists of branches that were filled and stained but do not belong to the cell of interest. In some cases, surrounding neurons take up the tracer molecule, and some of their dendrites or their main axon will be traced.
The second group comprises neuronal structures that were missed by the human tracer. The user usually starts tracing at the soma. If an axonal branch ending is missed in that section, the human tracer tends to ignore these areas in the subsequent sections and, therefore, sometimes misses large parts of the axonal branching patterns. In some cases, this error becomes obvious during splicing and the human tracer has to go back to the tissue and to identify the connections missed previously. The size of the group of missing branches closely correlates to the experience of the user. Only training, double checking, and anatomical knowledge will minimize the loss of branches and improve the tracing quality.
Comparing the two groups, an explanation for the large deviation from 4 to 70% in missed branches can be given. The cortical part of the axonal arbor from the thalamocortical cell is far away from the biocytin injection site (soma) in the thalamus.13 Therefore, no artificial branches from other cells were found in cortex. Furthermore, this cell was traced manually by an experienced person in the laboratory (T. K.). Both facts result in the small number (4%) of missed branches.
In contrast, the manual reconstruction of one cell in layers 2 and 3 was missing almost 70% of the automatically reconstructed axonal branches. This probably has two reasons. First, layers 2 and 3 are closer to the site of cell filling. Therefore, some artificial branches from unintentionally filled cells were traced automatically and left out intentionally by the manual tracer. The second reason is that this cell was the first reconstruction of axonal morphology this individual ever attempted.
In conclusion, the manual reconstruction of neurons and their axonal arbors using a computer-aided Camera Lucida system (e.g., Neurolucida) depends strongly on the abilities of the individual tracer, and reconstructions will be somewhat subjective.46 In contrast, the semiautomated tracing system delivers reliable and objective reconstructions even in the hands of new users (S. E.). S. E. was able to reconstruct the axonal arbor presented in Fig. 11b with no significant differences compared to a reconstruction of this neuron done by an experienced tracer (T. K.) (data not presented).
The automated system fails if the staining of the axon is insufficient, if the background is too high or in the presence of a strong gradient in thickness or contrast within the tissue. The latter will result in different optimal camera exposure times depending on the location of the field of view within the scanning pattern. Because the exposure is set to one value for the entire pattern, this may result in a failure of the subsequent deconvolution and neuron tracing. However, we are currently trying to compensate for this problem by applying of an exposure mask that compensates for inhomogeneous tissue.
Another limitation is that, due to the small diameters and tortuous paths of axons combined with physical resolution limits of the system, two close lying axons may occasionally be mistakenly connected by the automatic tracing algorithm. Because of the stereotyped nature of these errors, we may be able in the future to develop an algorithm to correct them. Nevertheless, because action potentials are propagated by saltatory conduction and the axon is surprisingly relatively electronically compact,47 modelers and experimentalists are often unconcerned with the precise locations of all branch points. Rather, the focus is usually on the length and spatial organization of the axonal arbor. Thus, despite this limitation, our current system already measures the two main parameters of interest with high reliability and accuracy.
Indeed, the box counting method yielded small deviations in shape between manually and automatically reconstructed branches. This difference mainly arises from the larger number of points used by the automatic system to represent the graph. At low resolution, this difference is almost 0. At high resolutions , close to the physical resolution limit, the relative deviation was approximately 5%. Therefore, depending on the scientific aim, these slight deviations can be irrelevant.
Nevertheless, Fig. 18 shows clearly that some branches, in particular at the highest resolution, deviate substantially from the 5% mean value given by the Gaussian fit. This group consists of dendritic branches. Because the dendrites are not as smooth as axonal branches, their automated reconstruction will also not be perfectly smooth. Reconstructed spine heads that were not removed during the graph validation and/or pruning, result in quite large relative shape deviations (Fig. 19 ). However, because the spines are part of the dendrites, the relative difference in shape can usually be neglected.
Mowiol shrinkage. Qualitative measurement of the shrinkage effect of tissue embedded in Mowiol. After 3h at room temperature, the specimen showed shrinkage of 2to4μm along the optical axis. After a few hours at 4°C this process is reverted. By keeping the room air-conditioned and using the narrow illumination bandpass filter that protects the specimen from any infrared light, the specimen shrinkage is slowed down. However, once the specimen is at thermal equilibrium with the surroundings, the shrinkage is extremely slow that it can be regarded static over the scanning period.
|Directly fromRefrigerator||3h at RoomTemperature||After 3.5h Recooling (4°C)|
An important issue is the reconstruction time needed to trace a neuron (e.g., the pyramidal neuron in layers 2 and 3 takes 10 to 30 brain slices). Because the automatic tracing is running a day and the human interaction time decreases from approximately , the reconstruction time decreases from weeks to days.
Comparison to Other Methods
Reconstructing the branching pattern of neurons is a challenging task. In the literature, several methods for different imaging and processing approaches have been described. 1, 2, 3, 4, 5, 6, 7, 8, 9, 11 Most of the methods available today are based on confocal or two-photon imaging. This is due to the good signal-to-noise ratio and the high resolution along the optical axis for both of these imaging techniques. Rodriguez9 showed that laser scanning microscopy in combination with sophisticated software algorithms is even capable of resolving the neuronal morphology of dendrites at spine level. Their method is capable of reconstructing the dendritic arbor and visualizing the surface of the spines. By scanning several adjacent image stacks, the field of view is increased, but still their method lacks the possibility of reconstructing the axonal arbor. This is due to the limited field of view of laser scanning microscopes, the bleaching of the tissue,48 and the lack of suitable mosaic systems for this imaging technique.
Another method that has been developed to analyze axonal arbors is the method described by Broser 15 Here lentiviral gene transfer is used to transfect a population of pyramidal neurons. The axonal arbors of these cells are subsequently stained. Two-dimensional mosaic microscopy in combination with extended focus imaging (EFI) is used to resolve the axons. Automatic 2D traces then quantify the axonal length in different areas of the tissue.
The appropriate reconstruction method strongly depends on the scientific task. To analyze spine shapes and distributions along a dendrite, the electron or confocal microscope methods described by Denk and Horstmann49 and Rodriguez 9 might be the most suitable. For determination of statistical connectivity patterns, the 2D EFI method proposed by Broser 15 could be used. However, if the precise 3D extension of the axonal arbor is important, one can either use the time-consuming Camera Lucida technique or one can use the alternative method presented here.
We found that the human influence on the reconstruction during the splicing of neuronal branches from individual sections may result in ambiguous tracings. To further minimize the subjective issue of pattern recognition, we currently are developing a semiautomated interactive 3D editing and splicing tool (AMIRA extension: SPATIAL GRAPH)43 in cooperation with the Konrad-Zuse Institute, which will replace the 2D editing environment of NEUROLUCIDA we used so far. The interface file format is the .hoc format (NEURON).50, 51 After editing and splicing in AMIRA, such a .hoc file can then be integrated into morphological databases or used for simulations in NEURON or NEURODUNE.52 Furthermore, an automated extraction of approximate dendritic radii will be added, because this is essential information for the simulation of neuron properties.
We conclude that the semiautomated method of 3D tracing of neuronal morphology is a reliable, objective, and less time-consuming alternative to the current manual systems. Furthermore, the quality of reconstruction depends less on the abilities of the human tracer. Our method should facilitate the creation of large databases of neuron morphologies. This will yield interesting possibilities for simulation of information processing in small neuronal circuits. Furthermore, axonal outgrowth under various conditions of development and experience can now be more easily quantified in three-dimensions and statistically analyzed.
Thanks to Marlies Kaiser for histology; Merle Kurz, Thorben Kurz (T.K.), and Maren Schmitt for manual reconstructions; Seszgin Erdogan (S.E.) for semiautomatic reconstructions; Enrico Reimer for large support in hardware and software optimization and automation of the pipeline; and Guenther Giese for his support in basic microscopy techniques.