Three-dimensional immersive virtual reality for studying cellular compartments in 3D models from EM preparations of neural tissues.

Advances in the application of electron microscopy (EM) to serial imaging are opening doors to new ways of analyzing cellular structure. New and improved algorithms and workflows for manual and semiautomated segmentation allow us to observe the spatial arrangement of the smallest cellular features with unprecedented detail in full three-dimensions. From larger samples, higher complexity models can be generated; however, they pose new challenges to data management and analysis. Here we review some currently available solutions and present our approach in detail. We use the fully immersive virtual reality (VR) environment CAVE (cave automatic virtual environment), a room in which we are able to project a cellular reconstruction and visualize in 3D, to step into a world created with Blender, a free, fully customizable 3D modeling software with NeuroMorph plug-ins for visualization and analysis of EM preparations of brain tissue. Our workflow allows for full and fast reconstructions of volumes of brain neuropil using ilastik, a software tool for semiautomated segmentation of EM stacks. With this visualization environment, we can walk into the model containing neuronal and astrocytic processes to study the spatial distribution of glycogen granules, a major energy source that is selectively stored in astrocytes. The use of CAVE was key to the observation of a nonrandom distribution of glycogen, and led us to develop tools to quantitatively analyze glycogen clustering and proximity to other subcellular features.

Over a century ago, neurohistologists such as Camillo Golgi and Santiago Ramón y Cajal formulated hypotheses about brain physiology based on detailed structural observations. Today, direct observations powered by high‐resolution imaging techniques remain an important tool to identify and analyze cellular and subcellular elements; however, recent advances in automated serial imaging are facilitating analyses of entire volumes of tissue samples with unprecedented detail (Hell and Wichmann, 1994; Denk and Horstmann, 2004; Hess et al., 2006; Rust et al., 2006; Knott et al., 2008; Briggman and Bock, 2012). These 3D imaging methods can produce precisely aligned image stacks that can be segmented to show subcellular features with high resolution (Jorstad et al., 2015).

Electron microscopy (EM) is the best technique to observe subcellular details of the nervous system: synaptic contacts, synaptic vesicles, endoplasmic reticulum, and mitochondria. Single‐plane transmission electron microscopy (TEM) images provides the best spatial resolution (to about 1 nm per pixel size), making it the technology of choice for studying brain ultrastructure. It has been a challenge to produce accurate 3D images that represent the complexity and level of detail at this resolution without it being a tedious and time‐consuming process.

Denk and Horstmann (2004) created the first commercial serial block‐face electron microscopy (SBEM) based on a concept proposed by Leighton (1981) that automatically cuts thin (<15‐nm) serial sections and creates images using backscatter detectors. The quality of the images matched that of TEM images, but they also aligned with each other, were nearly isotropic, and covered a broad field of view. Together with SBEM, focused ion beam scanning electron microscopy (FIB‐SEM) (Knott et al., 2008) was used as an alternative sectioning approach allowing thinner sectioning (down to 5 nm) in a more reliable manner, with a narrower field of view (around 15 µm). The two approaches are complementary, the first being more suitable for studying large field of view in detail (i.e., tracing long‐distance connections; Lichtman et al., 2014), and the second more suitable for looking at regional cellular features.

With the challenge of image production simplified, the focus moved to improving the management of large datasets and the segmentation of complex structures. 3D reconstruction and visualization of all elements included in the imaged volumes of tissues is pursued by all large‐scale brain projects such as the Human Brain Project (https://www.humanbrainproject.eu/), the BRAIN initiative (http://braininitiative.nih.gov/), and the Human Connectome project (http://www.humanconnectomeproject.org) (Kandel et al., 2013). In the past few years, this had led to an increased effort to develop engineering software suited for this task.

Small image stacks containing a limited number of objects can be reliably reconstructed manually, as long as the user correctly identifies the objects. However, whereas increasing image size requires only moderate additional access to computer resources, manually segmenting structures is considerably more difficult. The most common manual segmentation tools TrakEM2 (http://fiji.sc/TrakEM2, RRID:nlx_151924; Cardona et al., 2012) and reconstruct (http://synapses.clm.utexas.edu/tools/reconstruct/reconstruct.stm, RRID:nif‐0000‐23420; Fiala, 2005) can be very time‐consuming. In an effort to engage more manpower to share the load of manual segmentation work, scientists have developed collaborative annotation tools, mostly in the context of connectomics projects. Catmaid (www.catmaid.org; RRID:nlx_151922; Saalfeld et al., 2009) allows expert users to create a skeleton segmentation, but is not well suited for proofreading. Another tool, RhoANA, allows fully automated segmentation, scalable for very large datasets, that can then be proofread using Mojo, which is a stand‐alone software, or Dojo, which is a web‐based application allowing a collaborative proofreading framework (Kaynig et al., 2013; Knowles‐Barley et al., 2014, http://www.rhoana.org/dojo/). The Seung lab realized an interesting pipeline by turning the task into an online game (www.eyewire.org, RRID:nlx_144421), in which players are challenged to recognize and extend branches of neurons. Although the task may be somewhat repetitive, gamification encourages participation through rewards and competition between users (Helmstaedter et al., 2013).

Here, we took advantage of the semiautomated tool ilastik 0.5 (www.ilastik.org; Sommmer et al., 2011). The user initially identifies a structure's profile, which the program can use to follow the structure across an image stack to create a 3D reconstruction for proofreading. Here we describe our pipeline involving both TrakEM2 and the semiautomatic tool.

The recent technical solutions to 3D reconstructions that are making extremely detailed brain‐cell imaging available pave the way to quantitative analysis in 3D. The choice of tools that allow 3D analysis is very limited; for instance, Avizo or Imaris can produce measurements and annotations in 3D models. Here we have used the package NeuroMorph (www.cvlab.ch/NeuroMorph, RRID:SciRes_000156; Jorstad et al., 2015), which is freely available for measuring meshes in Blender (www.blender.org, RRID:nif‐0000‐31943), a 3D modeling software that offers the advantage of being customizable using a Python programming environment.

In particular, we used these tools to analyze the spatial relationship between glycogen granules contained in astrocytic processes and synaptic elements such as axonal boutons and dendritic spines. Blood‐derived glucose is stored in astrocytes as glycogen; the breakdown of glycogen (glycogenolysis) results in the production of lactate, a molecule that can be shuttled to neurons through a process known as the astrocyte–neuron lactate shuttle (Magistretti and Allaman, 2015). Recent studies show that glycogen‐derived lactate is necessary for synaptic plasticity (Yang et al., 2014) and memory formation in the hippocampus (Newman et al., 2011; Suzuki et al., 2011), but the predominant subcellular target, pre‐ and postsynaptic profiles, of lactate remains unknown. Understanding the polarization of glycogen toward different synaptic components could help identify which synaptic element is likely to be preferentially targeted by lactate.

The considerable detail involved in such a complex 3D model can make it difficult to turn observations into coherent hypotheses. For this reason, we took advantage of a cave automatic virtual environment (CAVE): a fully immersive virtual reality (VR) environment whose unique image projection allows the viewer to literally step inside and navigate within the reconstructed 3D model. The possibility of visualizing 3D cellular and subcellular data from such a perspective makes it a unique tool for explorative analysis of an apparently random and irregularly shaped environment like the brain. We propose that the use of VR is a crucial step in formulating appropriate questions and provides the means to address them.

Here we show how to conveniently observe and interactively analyze a highly complex dataset in a 3D virtual reality environment. In addition, we show that the approach is generally applicable for detailed studies of subcellular features (Fig. 10).

Figure 10

Visual schematics of the analysis pipeline. After the 3D model is extracted from serial electron micrographs, models can be loaded and visualized in CAVE at all times for multiple purposes (i.e., reconstruction proofreading, explorative analysis of the dataset). Quantifications can be done by means of custom‐made tools embedded in blender software. See text for a detailed description of how we observed a nonrandom distribution of the glycogen granules, and used a DBSCAN algorithm to infer the clustering of such cloud of points.

MATERIALS AND METHODS

Sample preparation and imaging

The EM stack was kindly provided by Graham Knott (BioEM imaging facility, EPFL, Lausanne) (Fig. 1), using the same method as described previously (Maco et al., 2014). Briefly, one adult (P90) Sprague–Dawley rat was deeply anesthetized by isoflurane inhalation and transcardially perfused using a fixative (2% paraformaldehyde [PFA], 2.5% glutaraldehyde [GA], EMS in PB 0.1 M 200 ml). The brain was removed 2 hours after perfusion, and 100‐μm coronal slices were cut by using a Leica VT1000 vibratome. Sections were collected in phosphate buffer (PB) 0.1 M, and those including the dorsal hippocampus were selected for staining and embedding in Durcupan (Fluka, Buchs, Switzerland) This procedure was in accordance with the Swiss Federal Laws on Animal Experimentation administered by the Véterinaire Cantonale Lausanne.

Figure 1

High‐resolution image stack from adult rat hippocampus. Note that orthogonal projections along the z‐axis are indistinguishable from the imaging plane (xy).

Tissue preparation

Sections were washed in cacodylate buffer (0.1 M, pH 7.4) and then postfixed and stained first in a reduced osmium solution (1.5% potassium ferrocyanide with 1% osmium tetroxide in cacodylate buffer) and then with 1% osmium tetroxide in cacodylate buffer alone for 40 minutes. Finally, sections were left in 1% aqueous uranyl acetate for 40 minutes and then briefly washed prior to dehydration.

Embedding in Durcupan

After membrane staining, tissues were embedded in Durcupan (Fluka/Sigma‐Aldrich, St. Louis, MO). Sections were dehydrated in aqueous solutions containing increasing concentrations of ethanol (50%, 70%, 96%, and 100%) prior to placing them into a 50% Durcupan–ethanol mix. The mix was then replaced with increasing concentrations of Durcupan until they reached a pure resin form. The tissues were then left overnight. The next morning, the sections were embedded in a thin layer of fresh resin in an aluminum weigh boat and placed in a 60°C oven for about 24 hours.

FIB‐SEM imaging

We used the following procedure to mount specimens with minimal specimen charging. Regions of the hippocampus were first dissected under a stereoscopic microscope using a razor blade. Next, samples were mounted onto a flat, blank resin slab and trimmed with a glass knife using an ultramicrotome in preparation for the electron microscope. The block containing hippocampal CA1 stratum radiatum was imaged inside an NVision 40 FIB‐SEM (Carl Zeiss) using an acceleration voltage of 1.5 kV, a current of 350 pA, and a dwell time of 10 μs/pixel. With the repeated use of a FIB, directed parallel to a block face to remove (or mill) a thin layer of embedded tissue, serial images were collected at a magnification of 6 nm/pixel. The milling depth between images was 5 nm, resulting in perfectly isotropic voxels. Each exposed surface was imaged using backscattered electrons.

3D Reconstruction, rendering, and analysis

Serial micrographs were registered using MultiStackReg (Thévenaz et al., 1998), a freely available plugin for Fiji software. To annotate and measure features of interest we took advantage of a combination of programs: TrakEM2 (http://fiji.sc/wiki/index.php/Fiji), a free plugin available from Fiji image processing platform (http://fiji.sc/TrakEM2, RRID:SciRes_000137; Cardona et al., 2012); ilastik (www.ilastik.org; Sommer et al., 2011) for segmentation and generation of the 3D models; and Blender (www.blender.org) for visualization and rendering with the NeuroMorph package (www.cvlab.ch/NeuroMorph; Jorstad et al., 2015).

We then expanded NeuroMorph by coding and embedding new add‐ons in Blender (see Fig. 6). The Blender/Python API enabled us to incorporate individual scripts into a single graphical interface within Blender (a so‐called add‐on), which is a Python module that provides a stable workflow of data analysis operations. The Glycogen Analysis plugin starts by calculating the centroids of the objects of interest (in our case glycogen granules or their clusters, spines, and boutons), calculates their distances, and identifies the closest spine or bouton from each glycogen granule or glycogen cluster. The results can be exported on a tab‐limited text file that contains the name of each measured couple of object and their distance, for further analysis. The Glycogen Analysis add‐on is freely available at https://github.com/daniJb/glyco‐analysis. It uses the following algorithms:

Figure 6

Visualization of glycogen clusters. Single glycogen granules (rendered as spheres) appeared as nonrandomly clustered when observed inside CAVE. Clusters have been generated using custom‐made tools embedded with Blender (left). Each cluster has a different color. Right, top; example of one dendrite (blue) and one axon (white) associated with two clusters; right, bottom; only granules and clusters are visible.

Nearest neighbor calculations were based on the KD Tree algorithm. It is a fast tree‐based data structure algorithm that indexes a k‐dimensional array of points that can be queried for nearest neighbors of any point (Bentley, 1975). Both algorithms are implemented in the scikit‐learn Python library (http://scikit‐learn.org/stable/modules/neighbors.html#k‐d‐tree).

Cluster analyses were based on the DBSCAN algorithm (Ester et al., 1996), an automated method to determine clusters from a set of points. This algorithm is well suited for noisy datasets because it assumes nothing about number or shape of clusters, but instead uses minimum distance between minimum number of members to determine clusters. These two parameters (minimum distance between members and minimum number of members) were optimized by using the silhouette coefficient (Rousseeuw, 1987), which gives a metric on the fraction of points included in all clusters. For a minimum number of members, the add‐on can look for the maximum silhouette coefficient within a range of minimum distances. For visual inspection, we colored the granules belonging to each cluster and enclosed them with similarly colored ellipsoids.

ilastik scaling

We tested the ilastik 0.5 carving workflow on a 1.7‐GB dataset and four progressively scaled‐down samples to quantify ilastik's performance on a 40‐core Intel Xeon CPU E5–2680 v2 @ 2.80 GHz with 64 GB RAM, using Scientific Linux 6.5. The results were extrapolated to a hypothetical stack of 26 GB with 4,096 × 4,096 × 1,500 pixels, ∼1.75 pixels/cell boundary, and covering ∼1003 μm to represent a minimum resolution for capturing all details of an entire astrocyte.

CAVE and zSpace

CAVE is a fully immersive virtual reality system that allows for the convenient and interactive visualization of 3D models. Our CAVE space is a cubic room with 3‐m sides that are each projected with two 8‐megapixel streams. With this setup, several people are able to simultaneously stand inside our 3D model, where details are 1 million‐fold larger (μ to m).

We also took advantage of zSpace (http://zspace.com/), a monitor with enhanced 3D capabilities, using a head‐tracking system to project a holographic‐like stereo image. Although the field of view of the image was limited, it was small and portable, able to be conveniently placed on any workspace desk for interactive segmentation proofreading.

The 3D model was projected directly from the 3D modeling Blender window interface by using TechViz software (www.TechViz.net). The user is then able to manipulate the model inside the room, or within the zSpace, with a controller (see Supplementary Material). A head‐tracking system can modify the prospective of the object based on the position of the user, who is able to move freely in the room and observe the model from different points by keeping the model still and moving around it.

RESULTS

Image stack

Isotropic EM stack (Fig. 1) was acquired through hippocampal CA1 stratum radiatum. This allowed us to image a 220 μm3 volume (7.07 μm × 6.75 μm × 4.75 μm) at a voxel resolution of 6 nm.

The imaged volume contained portions of axons, dendrites, a large blood vessel surrounded by endothelium and pericytes, astrocytic processes with intracellular glycogen granules, and synapses (Figs. 2, 3). To measure the distribution of granules compared with the other components in the EM stack, we segmented and reconstructed each element of interest in 3D (Supplementary video 1). Axons and dendrites have a well‐known morphology. Axons appear as relatively straight, tubular structures, with a narrow diameter, swelling at sites where spheroidal vesicles accumulate (axonal boutons) (Shepherd and Harris, 1998). Dendrites appear as tubular structures, with a larger and regular diameter compared with axons, often showing protrusions where typically synaptic contacts occur, called dendritic spines (Stuart and Häusser, 2007). Glial processes were identified by their strategic localization among neurites, the presence of glycogen granules, and the relatively clear cytosol (Witcher et al., 2007, 2010). Portions of endothelial cells were covering the lumen of the blood vessel, and processes of pericytes were at the interface between the brain parenchyma and the endothelial cells (Hamilton et al., 2010) (Fig. 4). Serial sections were useful to clarify doubts about the nature of the different structures, which were annotated in the resulting 3D model (Fig. 5) using terms from public ontologies (Gene Ontology Consortium et al., 2013; Meehan et al., 2011) (see Supplementary Material).

Figure 2

TrakEM2 interface. One synapse (red) and glycogen granules (yellow spheres, arrowhead) have been segmented in one plane. Scale bar = 500 nm.

Figure 3

Identification of synapses and glycogen granules. Examples of synapses (A,D; red arrows), glycogen granules (D, yellow arrows), axons (A), boutons (B), dendrites (D), spines (s), and astrocytic processes (Ast) on EM micrographs. E: The lumen of the blood vessel is surrounded by one endothelial cell interfacing with a pericytic process. Excitatory asymmetric synapses (B) or inhibitory symmetric synapses (C) can be discriminated by the thickness of the synaptic density, as well as the shape and size of synaptic vesicles. Scale bar = 500 nm in A,D,E; 300 nm in B,C.

Figure 4

ilastik 0.5 Interface. The down‐sampled stack has been loaded and one dendrite (green) and one axon (yellow) have been segmented over all the serial images and then rendered in 3D.

Figure 5

3D model from a fully reconstructed volume of rat brain created from the image stack. All elements (axons, white; dendrites, blue; astrocyte, green; capillary, red; synaptic densities, red) were combined and rendered using Blender.

NeuroMorph add‐ons. The import add‐on sets the right scale for models (top image, bottom inset). The superimposition add‐on shows the selected EM source image together with the model (bottom image). The measure add‐on does spatial measurements (lengths, surfaces, and volumes) on user‐selected objects (bottom image).

Glycogen and synapse segmentation

We used TrakEM2 to segment synaptic densities and glycogen granules (Fig. 2).

Synaptic densities appear as a dark (i.e., electrodense) membrane thickening between axonal boutons and postsynaptic spines or shafts (Figs. 2, red, 3A,D, red arrows). Synapses were classified based on their appearance and the shape and size of the vesicles in the presynaptic terminals (excitatory ones have a thick asymmetric postsynaptic density and presynaptic terminals with round clear vesicles, whereas inhibitory ones have thin symmetric postsynaptic density with flattened vesicles at the presynaptic site; Fig. 3B,C). Out of 93 synapses in this sample, 58 were on spines and 35 were on shafts. The majority (79) were presumed excitatory, 35 of which were on shafts, whereas 12 were inhibitory, 9 of which were on shafts. Only two synapses were unidentifiable based on these morphological features, with pleiomorphic vesicles at a presynaptic site, and we classified them as “unknown.” Each synaptic density was manually marked in each serial image.

Glycogen granules appear inside astrocytic processes as small, dark spheres with diameters ranging from 10 to 60 nm (Figs. 2, yellow sphere, 3D, yellow arrows). Even though the resolution of the microscope was 6 nm/voxel, a feature that makes this technique less appropriate for analysis of intracellular organelles, the size distribution was appropriate for imaging glycogen granules using FIB‐SEM. We used the ball tool in TrakEM2 to render individual granules as spheres. We identified a total of 1,007 granules in this stack.

Cellular segmentation

All other elements (axons, dendrites, blood vessels, and astrocytic endothelial and pericytic processes) were reconstructed by using the carving tool from ilastik version 0.5 (Fig. 4). The underlying algorithm is described in detail in Straehle et al. (2011). The stack contained 11 dendrites with a total of 56 spines, 161 axons whose classification is based on previous assumption (68 excitatory, 9 inhibitory, 67 unknown, and 17 showing no varicosities), and 94 boutons (73 excitatory, 10 inhibitory, and 12 unknown). Out of the 12 unknown axonal profiles with identifiable varicosities, only two were associated with a synaptic density. The remaining ones were en passant terminals, with small, pleiomorphic flattened vesicles, with or without dense‐core granules, not associated with any synaptic density, whose morphology is consistent with that of monoaminergic terminals (Gaugler et al., 2012). Non‐neuronal processes included one astrocytic process, an epithelial cell enclosing a blood vessel, and two pericytic processes interfacing the epithelium and the neuropil.

Glycogen distribution

Analysis of glycogen granule distribution showed that glycogen is more frequently found around axonal boutons (total count 523) than dendritic spines (158). The vast majority of granules were nearest to excitatory boutons (347), followed by unknown boutons (97), and inhibitory boutons (79) (Fig. 8B). The remaining granules were closest to the perivascular cells (115 to pericytes and 211 to endothelial cells; Fig. 8A). We then counted the number of granules closest to each individual spine and bouton: spines had the lowest number of granules, with an average of 5.9 ± 1.4, followed by excitatory boutons (7.9 ± 1.1) and inhibitory boutons (8.8 ± 2.3). The highest number of granules colocalized was with unknown varicosities (19.4 ± 6.8).

Figure 8

Distribution of glycogen granules around synaptic features. A,B: Analysis on single granules. C,D: Analysis on clusters.

The parameters of the clustering algorithm DBSCAN (Ester et al., 1996) (0.31 μm maximum distance for adjacent granules in a cluster and a minimum cluster size of 7; Fig. 6) were optimized with the silhouette coefficient (Rousseeuw, 1987) for the maximum fraction of granules within clusters. Of the 24 resulting clusters, 14 (52% of the total) were associated with axonal boutons, for a total of 571 granules; the remaining 10 (containing 319 granules) were associated with dendritic spines (Fig. 8C,D). The preference toward presynaptic terminals confirms previous observations.

DISCUSSION

We have described a workflow that is suitable for analyzing the distribution of neural profiles in 3D space. The high resolution of the imaging procedures allowed us to identify and precisely locate axons, dendrites, astrocytic processes, and synaptic connections, as well as minute subcellular features such as glycogen granules. We performed an initial study on the distribution of this energy source in relation to specific cellular compartments.

Analysis of glycogen in a 3D volume

Although the creation of 3D models is becoming increasingly common, limited tools are available for detailed morphometric 3D analyses. The NeuroMorph package (Jorstad et al., 2015) for Blender analyzes lengths, surfaces, and volumes of meshes representing 3D objects. We expanded this analytic capability with our own plugin: Glycogen Analysis. Our plugin calculates distances either between points in space and meshes or between points and centroids of several objects, while identifying the nearest neighbor. With this approach, it is possible to analyze glycogen distribution by calculating 1) the distance between single granules and their closest individual element; or 2) the distance between the center of the cluster and the centroid of each synaptic element at granule clusters.

Imaging a volume of tissue using FIB‐SEM

EM is the gold standard microscopy technology to identify cellular and subcellular features with nanometer resolution. Single‐planar images originating from tissue slices that have been cut in series can be registered and aligned for further processing to reconstruct 3D volumes. The classic approach to the successful acquisition of large sets of serial section electron micrographs (SSEMs) consists of manually cutting tissues from a resin bloc and imaging them under a TEM microscope. However, SSEM is challenging and time consuming even for experienced microscopists. Here, we took advantage of FIB‐SEM microscopy, an imaging technique that automatically produces serial electron micrographs at near‐isotropic resolution. Compared with classic TEM, FIB‐SEM takes advantage of backscattered electrons to generate the image. Resin‐embedded sample block face is sequentially scanned using a low‐voltage electron beam. Although sample preparation and image acquisition and registration in FIB‐SEM requires considerably less time than in SEM, the number of images and hence volume are considerably larger, making image processing the most demanding step in the process.

Caveats of FIB‐SEM for intracellular organelles

Before taking into consideration the use of this technique for imaging subcellular organelles, a few caveats must be considered. First, the spatial resolution of block‐face EM does not reach TEM resolution; therefore imaging of nanometer‐scale organelles is a problem.

More importantly, block‐face SEM fails to image several nonmembranous subcellular compartments or organelles. One example is provided by ribosomes, which are simply not visible using block‐face SEM. Intracellular filaments, such as actin bundles, have also been reported to be hard to distinguish using block‐face SEM. The reason might be the higher detected electron density of the cytosolic compartments; to overcome this issue, special techniques are required to enhance the contrast of these compartments (Tanaka and Mitsushima, 1984; Tanaka et al., 1986). Nevertheless, the convenience of this imaging technique is driving researchers, mostly in the field of cell biology, to find solutions that allow the use of this approach for intracellular studies; one solution for instance is provided by comparative SEM–TEM studies (Puhka et al., 2012; Watanabe et al., 2012; Ichimura et al., 2015).

In our specific case, the size of the glycogen granules (about 10 nm minimum) and the pixel magnification we used (6 nm) allowed us to reliably detect even the smallest granules. TEM has been successfully used on glycogen detection by others (Marchand et al., 2002; Obel et al., 2012).

Because the stack is not stained for any particular structure, cellular features need to be recognized by their morphology, meaning that reconstruction work requires supervision by an expert who must proofread the results. Although pre‐embedding immunohistochemistry using diaminobenzidine (DAB) or gold particles to reveal primary antibodies has been reported (Sonomura et al., 2013), permeabilization of tissue might impair the fine structural preservation of the specimen. In our sample, the portion of the astrocytic process surrounding the capillary is filled with granules, a site where glycogen is known to accumulate (Broadwell et al., 1983; Cataldo and Broadwell, 1986). Our analysis confirmed that a significant portion of the total number of granules is associated with the blood vessel wall, suggesting that the granules formed following uptake of glucose from the bloodstream (Fig. 8). A recent article (Takado et al., 2015) reports the accumulation of 13C‐labeled glucose into glycogen granules that could be imaged either by NanoSIMS imaging and EM, revealing that 13C‐labeled glucose accumulates in particles that have the morphology of glycogen granules. Taken together, these data indicate that the particles characterized here are to be considered as glycogen granules, without any staining, and can be identified solely by the morphological features discussed above.

Manual versus semiautomated segmentation

Our primary challenge was to reduce the time‐consuming process (Kreshuk et al., 2011) of generating a large‐scale 3D model that includes all elements potentially interacting with glycogen granules. Fully automated approaches have not delivered reliable results, even for common structures like axons and dendrites (Chklovskii et al., 2010; Mishchenko et al., 2010). Semiautomated segmentation tools are a reasonable compromise between speed and reliability. The user is only required to identify the inside and outside of the structures of interest in a few sections for the software to automatically follow membrane boundaries along the z‐axis to create a draft of the segmentation.

TrakEM2 is currently the most widely used tool for the segmentation of EM stacks. Because it is entirely manually operated, each segmentation step can immediately be validated. Provided that the datasets are small, the carving module of ilastik can be an excellent time‐saving tool for image segmentation.

Compared with traditional manual approaches that require manual drawing of thousands of serial images to achieve segmentation, ilastik only requires a few labeling clicks on the structure of interest to segment it. The work involving segmentation of a few hundred objects over hundreds of serial images—a task that might take many months using TrakEM2—can instead be done by a single user in less than a week using ilastik.

The ilastik scaling issue

Two primary factors determine the minimum EM‐image stack size: the resolution of the boundary features (which must be >1 pixel, but more is better) and the total volume in pixels.

Boundary feature resolution is important because it controls the accuracy of the seeded watershed algorithm and determines the number of training labels needed. More distinct, and, up to a point, larger boundaries are better. Artifacts from the original staining and imaging process also affect the results. For example, if the cell boundaries are not completely continuous, the segmentation bleeds through the gaps; this can only be repaired by manual intervention.

Computing time increases linearly with the volume of the image stack (i.e., cubically with stack dimensions; Fig. 9). A good practice is to down‐sample the stack, toreach the optimal balance among image resolution, preprocessing time, and segmentation accuracy. After down‐sampling our high‐resolution stack (1,178 × 1,125 pixels, 789 serial images) by 2.5 times, the subcellular details were still identifiable, and the stack preprocessing time was reduced to a reasonable 10 minutes using a laptop. Segmentation work created 160 axons, 10 dendrites, 1 astrocytic process, 2 endothelial processes, and 1 pericytic process. The semiautomated approach with ilastik enabled one person to finish the work in 4 days.

Figure 9

ilastik performance as function of size, memory usage (A), and processing time (B). Run on a 40‐core (Intel Xeon E5–2680 v2 @2.80GHz) workstation 64 GB RAM under Scientific Linux 6.5.

Applications of ilastik version 0.5 are limited by its memory and computational requirements. The seeded watershed algorithm uses a complete graph of all connected supervoxels in an image stack, and its current implementation requires that the entire image and supervoxel data be kept in memory. Because image stack size increases cubically when the sides of the captured stack grow, the need for memory and CPU power quickly exceeds any existing workstation‐class hardware (see Materials and Methods). By extrapolating performance from a downscaled stack, it appears that to capture all the details of an entire astrocytic cell would require a 26‐GB stack, which would call for ∼1,000 GB to preprocess, ∼700 GB to segment, another ∼24 hours to preprocess, and ∼15 minutes for interactive label segmentation. When analyzing our sample, we bypassed this problem by down‐sampling the stack because the resolution of the stack was more than adequate to resolve all the features.

Running ilastik on a multicore computer had little effect because most of the processing of the carving module runs single threaded. This limitation is currently being addressed such that future versions of ilastik will parallelize its workload.

Why is it difficult to describe a complex 3D model?

The resulting 3D model (Fig. 5) appears as a very complex forest of connected axons (Fig. 5, white) and dendrites (Fig. 5, blue), packed within and around a complex astrocytic perivascular endfoot (Fig. 5, green) that closely envelops them. Because of their complexity, qualitative observations are difficult, even for an experienced neuroscientist. Nevertheless, there are obvious benefits from any system that allows a clear and complete view of the entire model.

Benefits of CAVE and zSpace systems

In our work, CAVE was pivotal to assess the spatial arrangement of glycogen granules (Supplementary Video 2). The use of an interactive VR room and the more portable system zSpace (Supplementary Video 3) allows multiple users to share and discuss minute details, potentially setting up a new standard in the observation of densely populated neurobiology datasets. First, the granules appear as a collection of random points in space, but upon further observation when fully immersed inside the virtual 3D space, the cloud of granules appears to be grouped in various sized clusters. Granules and clusters of granules also appear to have a certain spatial relationship with spines, boutons, or the blood vessel, when visualized together. These initial observations led us to design an analytical strategy to unveil the spatial association between glycogen granules or clusters, and the different cellular elements (Fig. 6).

Analysis

Our results indicate that significant portions of granules are within astrocytic profiles associated with the vasculature, a fact that may reflect the immediate conversion of glucose entering the brain parenchyma from the bloodstream. Furthermore, energy stores, represented by glycogen granules, accumulate preferentially around presynaptic boutons rather than dendritic spines. Bouton distribution analysis revealed that boutons of uncertain nature, most likely monoaminergic varicosities, are associated with energy stores, suggesting a functional interaction between the two elements. This is not surprising, as monoamines such as noradrenaline and serotonin modulate glycogen metabolism in the brain (Magistretti and Allaman, 2015).

The clustering reveals a nonrandom distribution of glycogen granules, suggesting the existence of cellular mechanisms that regulate the process. We can speculate that two separate processes, each with its own role, are involved—one triggered by high glucose concentration near the blood vessel and resulting in glycogen storage, and the other resulting in the accumulation of glycogen in the proximity of neuronal elements that contain neuromodulators that can mobilize glycogen upon increased energy demands (Phelps, 1972; Roach et al., 2012; Magistretti and Allaman, 2013).

CONCLUSIONS

We have described a series of procedures and software applications that allow 3D reconstruction of an entire brain volume from serial EM sections, and how to take full advantage of such a model using a fully immersive VR facility (i.e., CAVE) as a powerful explorative tool. Also, because ilastik or TrakEM2 could process any kind of image stack and generate meshes that Blender could process, our pipeline could be applied to any imaging technology producing an image stack. The narrow field of view of a 2D monitor cannot compare with the advantage of being immersed in a complex dataset, similarly to an astronomer looking at the stars on the monitor of the screen rather than in the open air through a powerful telescope. In addition to providing an unprecedented opportunity to examine the relationships of cellular and subcellular elements of the neuropil, with particular emphasis on relationships at the neuron–glia–vascular unit level, this approach has allowed us to identify a particular association between glycogen clusters contained in astrocytic profiles and presynaptic boutons of likely monoaminergic varicosities. Our approach represents a useful tool to analyze glycogen localization and dynamics, particularly in behavioral paradigms in which its role has already been demonstrated, such as for synaptic plasticity, memory consolidation, and sleep–wake regulation (Suzuki et al., 2011; Magistretti and Allaman, 2015; Petit et al., 2015).

CONFLICT OF INTEREST STATEMENT

The authors declare no conflicts of interest.

ROLE OF AUTHORS

All authors had full access to all the data in the study and take responsibility for the integrity of the data and the accuracy of the data analysis. CC wrote the manuscript; CC and JB worked on the reconstruction; CC, JB, HL, and DB performed the analysis; DB, JB, and HL wrote the analysis scripts; MS and GH supervised the virtual reality environments; GH tested the performance of ilastik; AK and FAH wrote ilastik 0.5 software; CC, HL, and PJM initiated and supervised the project and contributed to the manuscript.

Supplementary Information Movie 1

Click here for additional data file.

Supplementary Information Movie 2

Click here for additional data file.

Supplementary Information Movie 3

Click here for additional data file.