Research Article

Dense connectomic reconstruction in layer 4 of the somatosensory cortex

See allHide authors and affiliations

Science  29 Nov 2019:
Vol. 366, Issue 6469, eaay3134
DOI: 10.1126/science.aay3134

Brain anatomy revealed in startling detail

The mammalian cerebral cortex is an enormously complex network of neuronal processes that are long and thin, branching, and extremely densely packed. This high packing density has made the reconstruction of cortical neuronal networks challenging. Motta et al. used advanced automated imaging and analysis tools to reconstruct with high spatial resolution the morphological features of 89 neurons and their connections in the mouse barrel cortex. The reconstruction covered an area more than two orders of magnitude larger than earlier neuroanatomical mapping attempts. This approach revealed information about the connectivity of inhibitory and excitatory synapses of corticocortical as well as excitatory thalamocortical connections.

Science, this issue p. eaay3134

Structured Abstract

INTRODUCTION

The brain of mammals consists of an enormously dense network of neuronal wires: the axons and dendrites of nerve cells. Their packing density is so high that light-based imaging methods have so far only been able to resolve a very small fraction of nerve cells and their interaction sites, the synapses, in mammalian cortex. Recent advances in three-dimensional (3D) electron microscopy allow researchers to image every nerve cell and all chemical synapses in a given piece of brain tissue, opening up the possibility of mapping neuronal networks densely, not just sparsely. Although there have been substantial advances in imaging speed, the analysis of such 3D image data is still the limiting step. Therefore, dense reconstructions of cortical tissue have thus far been limited to femtoliter-scale volumes, keeping the systematic analysis of axons, neuronal cell bodies and their dendrites of different types, and the dense connectome between them out of reach.

RATIONALE

Image analysis has made decisive progress using artificial intelligence–based methods, but the resulting reconstructions of dense nerve tissue are still too error-prone to be scientifically meaningful as is. To address this, human data analysis has been integrated into the generation of connectomes and it is the efficiency of this human–machine data analysis that now determines progress in connectomics. We therefore focused on efficiency gains by: (i) improving the automated segmentation quality, (ii) analyzing the automated segmentation for locations of likely errors and directing the human work to these locations only, and (iii) optimizing human data interaction by helping annotators to immediately understand the problem to be solved, allowing fast, in-browser parallel data flight, and by minimizing latency between annotator queries. With this, close to 100 student annotators solved hundreds of thousands of reconstruction problems within just 29 s each, including all preparation and transition time.

RESULTS

We reconstructed 2.7 m of neuronal wires densely in layer 4 of mouse somatosensory cortex within only ~4000 invested human work hours, yielding a reconstruction ~300 times larger than previous dense cortical reconstructions at ~20-fold increased efficiency, a leap for the dense reconstruction of connectomes. The resulting connectome between 6979 presynaptic and 3719 postsynaptic neurites with at least 10 synapses each, comprising 153,171 synapses total, was then analyzed for the dense circuit structure in the cerebral cortex. We found that connectomic data alone allowed the definition of inhibitory axon types that showed established principles of synaptic specificity for subcellular postsynaptic compartments, but that at scales beyond ~5 μm, geometric predictability of the circuit structure was low and coarser models of random wiring needed to be rejected for dense cortical neuropil. A gradient of thalamocortical synapse density along the cortical axis yielded an enhanced variability of synaptic input composition at the level of single L4 cell dendrites. Finally, we quantified connectomic imprints consistent with Hebbian synaptic weight adaptation, obtaining upper bounds for the fraction of the circuit that could have undergone long-term potentiation.

CONCLUSION

By leveraging human–machine interaction for connectomic analysis of neuronal tissue, we acquired the largest connectome from the cerebral cortex to date. Using these data for connectomic cell-type definition and the mapping of upper bounds for the learned circuit fraction, we establish an approach for connectomic phenotyping of local dense neuronal circuitry in the mammalian cortex, opening the possibility for the connectomic screening of nervous tissue from various cortices, layers, species, developmental stages, sensory experience, and disease conditions.

Dense reconstruction of ~500,000 cubic micrometers of cortical tissue yielding 2.7 m of neuronal cables (~3% shown, front) implementing a connectome of ~400,000 synapses between 34,221 axons and 11,400 postsynaptic processes (fraction shown, back).

These data were used for connectomic cell-type definition, geometrical circuit analysis, and measurement of the possible plastic fraction (the “learnedness”) of the circuit.

Abstract

The dense circuit structure of mammalian cerebral cortex is still unknown. With developments in three-dimensional electron microscopy, the imaging of sizable volumes of neuropil has become possible, but dense reconstruction of connectomes is the limiting step. We reconstructed a volume of ~500,000 cubic micrometers from layer 4 of mouse barrel cortex, ~300 times larger than previous dense reconstructions from the mammalian cerebral cortex. The connectomic data allowed the extraction of inhibitory and excitatory neuron subtypes that were not predictable from geometric information. We quantified connectomic imprints consistent with Hebbian synaptic weight adaptation, which yielded upper bounds for the fraction of the circuit consistent with saturated long-term potentiation. These data establish an approach for the locally dense connectomic phenotyping of neuronal circuitry in the mammalian cortex.

The cerebral cortex of mammals houses an enormously complex intercellular interaction network implemented with neuronal processes that are long and thin, branching, and extremely densely packed. Early estimates indicated that 4 km of axons and 400 m of dendrites are compressed into a cubic millimeter of cortical tissue (1). This high packing density of cellular processes has made the locally dense mapping of neuronal networks in the cerebral cortex challenging.

So far, reconstructions of cortical tissue have been either sparse (27) or restricted to small volumes of up to 1500 μm3 (810). Consequently, the detailed network architecture of the cerebral cortex is unknown. Particular open questions are to what degree local neuronal circuits are explainable by geometric rules alone (1, 2, 1113) and on which spatial scales cortical connectivity is only explainable by innervation preferences beyond such geometric models (5, 8, 9, 14, 15). Similarly, although numerous cortical neuronal cell types have been described based on protein expression, morphology, and electrophysiological characteristics (16), and these have been shown to have particular synaptic target patterns (17), the inverse question—whether, at the level of the dense cortical circuit, axons represent a continuum of synaptic preference or a set of distinct innervation paradigms that would allow for a purely connectomic cell type definition [as has been successful in the retina (18, 19)]—is still open. Next, at the level of synaptic input to the primary dendrites of cortical excitatory cells, it is not known whether the typically three to 10 primary dendrites of a cortical neuron that leave the cell body homogeneously sample the available excitatory and inhibitory synaptic inputs or if there is an enhanced heterogeneity of synaptic input composition, making it possible to exploit the numerous mechanisms that have been discussed for the nonlinear integration of local synaptic inputs (2023). Finally, whereas the change of synaptic weights in response to electrical and sensory stimulation has been widely studied (2428) and connectomic data consistent with LTP have been described (29, 30), the fraction of a given cortical circuit that is plausibly shaped by processes related to Hebbian learning under undisturbed conditions is still unknown.

We used dense connectomic reconstruction to quantitatively address these questions about the formational principles of a dense cortical circuit.

Results

We acquired a three-dimensional (3D) EM dataset from upper layer 4 of primary somatosensory cortex of a 28-day-old mouse (Fig. 1, A to D, likely located within a barrel, see supplementary materials) using serial block-face electron microscopy [SBEM (31); dataset size: 61.8 × 94.8 × 92.6 μm3; voxel size: 11.24 × 11.24 × 28 nm3]. For dense reconstruction (Fig. 1, E to H), we 3D aligned the images and applied a sequence of automated analyses [SegEM (32), SynEM (33), ConnectEM, and TypeEM; Fig. 2, supplementary materials and methods, and table S2], followed by focused manual annotation (FocusEM). We reconstructed 89 neurons that had their cell body in the dataset (Fig. 1, E and F). These neurons constituted only 2.6% of the total path length (69 mm; Fig. 1G). To reconstruct axons, which constitute most of the wiring in the dense circuit (1.79 m, 66.6%, Fig. 1H), we applied a scalable distributed annotation strategy that identified locations of uncertainty in the automated reconstruction, which were then resolved by targeted manual annotation. To reduce the required manual annotation time, it was critical to obtain an automated reconstruction with low error rates, to use efficient algorithms for identifying locations for focused manual inspection (queries), and to minimize the time spent per user query. For this (Fig. 2A), we developed artificial intelligence–based algorithms that evaluated the EM image data and convolutional neural network (CNN)–filtered versions of the image data in the surrounding of interjunctions between segmented pieces of neurites (Fig. 2B). Together with classifiers that computed the probability of volume segments belonging to an axon, a dendrite, a spine head, or a glial process (using, among others, shape features; Fig. 2C), this allowed us to automatically connect parts of dendrites, attach spine heads to dendritic shafts (by a greedy stepwise agglomeration initiated at the spine head, Fig. 2D; 58.9% of spine heads unaffected by the dataset boundary were automatically attached), and reconstruct parts of axons. Similarly, synapses were automatically detected by evaluating pre- and postsynaptic volumes at neurite interfaces [Fig. 2E and figs. S2 to S5 (33); for shaft synapses, additional CNN-based classifiers for vesicle clouds and mitochondria were used]. To manually correct remaining errors in axons (Fig. 2, F to H), we detected ending locations of automatically reconstructed axon pieces (Fig. 2F) and directed user queries to these locations. For this, we used an egocentric 3D image display mode [“flight mode,” Fig. 2G (34)] and oriented the user annotation along the axis of the neurite for which a local annotation (“query”) was requested (movie S2). Together with data preloading, this yielded a low-latency, targeted neurite annotation in which individual user queries took 29.4 s to resolve (traveled path length per query: 5.49 μm). These queries could be easily distributed among 87 annotators. Similarly, we detected locations of likely mergers between axons (Fig. 2H, “chiasmata”) and directed user queries to reconnect the chiasma exits along actual axons. Using this scalable annotation architecture, we obtained a dense reconstruction of 2.69 m of neuronal processes (Fig. 1, G and H) with a total investment of 3981 human work hours, ~10 times faster than a recent dense reconstruction in the fly larval brain (35) (Fig. 2, I and J), ~20 times faster than the previous dense reconstruction in the mammalian retina (18), and ~25 times faster than the previous dense reconstruction in mammalian cortex (9) (Fig. 2, I and J). To quantify remaining reconstruction error rates in this dense neuropil reconstruction, we measured the remaining errors in a set of 10 randomly chosen axons and found 12.8 errors per millimeter of path length (of these, there were 8.7 continuity errors per millimeter; see materials and methods). This is indistinguishable from the error rates previously found in fast human annotations (18, 34, 36).

Fig. 1 Dense connectomic reconstruction of cortical neuropil from layer 4 of mouse primary somatosensory cortex.

(A to D) Location [(A), red] of the 3D EM dataset (B). WM, white matter. High-resolution example images are shown in (C) and (D). Asterisks indicate examples of dendritic spines. Direct links to data browser webKnossos are as follows: https://wklink.org/9276 (B), https://wklink.org/7101 (C), and https://wklink.org/8906 (D). (E) Reconstruction of n = 89 neurons with a cell body and dendrites in the dataset. (F) Three spiny neurons (SpNs) and two INs (see movie S1). (G) Quantification of circuit components in the dense reconstruction. Most of the circuit path length (total: 2.69 m) is contributed by nonproximal axons (1.79 m, 66.6%), spine necks (0.55 m, 20.5%), and dendritic shafts (0.28 m, 10.3%) not connected to any cell body in the volume. (H) Display of all 34,221 reconstructed axons contained in the dataset. Scale bars in (D) are as in (C); scale bar in (F) is 10 μm.

Fig. 2 Methods for the efficient dense connectomic reconstruction.

(A) Simplified diagram of reconstruction steps [fig. S1, detailed in (B) to (H)]. wh, annotation work hours. (B) ConnectEM classifier for combining neurite pieces from the CNN-based volume segmentation (32): at junctions of volume segments (bottom right), raw data, CNN, and shape features were evaluated. (C) TypeEM classifier for assigning cellular identity to volume segments: the probability of axons, dendrites, spine heads, and glial processes was evaluated. Shown is an illustration of spine head (purple) and astrocyte (cyan) classification; one of the 985 features is illustrated (segment thickness). Numbers indicate the probability of the segment being a spine head. Precision and recall of spine head detection were 92.6 and 94.4%, respectively. (D) Process for automatically attaching spine heads to the dendritic shaft by stepwise agglomeration of volume segments along the highest-probability transition between neighboring segments [according to the ConnectEM score (B)]. An example of six neighboring spine heads that were all automatically attached is shown. In total, 58.9% of spine heads were automatically attached (A). (E) Automated detection of spine and shaft synapses [here, vesicle clouds (green) and mitochondria (blue) were detected and used as additional features for the SynEM (33) classifier]. (F to H) Focused annotation strategy for directing human annotation queries (Q) to ending locations of the automatically reconstructed axon pieces [(F), blue], oriented along the axon’s main axis [traced in webKnossos using flight mode (G {34}), yielding flight paths of 5.5 ± 8.8 μm length (21.3 ± 36.1 s per ending annotation, n = 242,271, movie S2)]. Neurite mergers (H) were detected as “chiasmatic” configurations, and queries (Q) directed from the exits of the chiasma toward its center were used to determine correct neurite continuities (fig. S1). (I and J) Quantification of circuit size and invested work hours for dense circuit reconstructions in connectomics and resulting order-of-magnitude improvement provided by FocusEM compared with previous dense reconstructions (m). Fish o.b., zebrafish olfactory bulb (59); M. retina, mouse retina IPL (18); Fly larva, mushroom body in larval stage of D. melanogaster (35); M. cortex, mouse somatosensory cortex [(9) and this study (magenta)]. Only completed dense reconstructions were included in the comparison.

We obtained a connectome (Fig. 3) between 34,221 presynaptic axonal processes and 11,400 postsynaptic processes [6979 × 3719 connectivity matrix (Fig. 3E) when restricted to those pre- and postsynaptic neurites that established at least 10 synapses]. Among the postsynaptic processes, we classified n = 169 apical dendrites (ADs) that traversed the dataset along the cortical axis without connection to one of the neuronal cell bodies in the dataset (Fig. 2A), 246 smooth dendrites (SDs, Fig. 2B), 80 somata, 116 axon initial segments (AISs; Fig. 2C), and 89 proximal dendrite (PD) trees connected to a soma in the dataset (movie S1; note that some of these neurons also had ADs that were classified as PDs and not included in the AD definition above; see materials and methods and tables S1 and S2).

Fig. 3 Postsynaptic target classes and dense cortical connectome.

(A to D) Display of all ADs [(A), magnified one-AD bundle (left) and top view in tangential plane illustrating AD bundles], SDs [(B), magnification inset illustrating low rate of spines], AISs (C), and their respective path length and spine density distributions (D). Note that spine density is underestimated by ~20% (table S1). (E) Display of connectome between all axons (n = 6979) and postsynaptic targets (n = 3719) in the volume with at least 10 synapses each, establishing a total of 153,171 synapses (of 388,554 synapses detected in the volume). For the definition of postsynaptic target classes, see (A) to (D); for the definition of presynaptic axon classes, see Fig. 4 and fig. S6. AISs with fewer than 10 input synapses are also shown. SOM, neuronal somata; Note that some of these PD dendrites are L4 ADs not included in the AD definition above. Asterisks indicate remaining unassigned axons.

Connectomic definition of axon types

We investigated whether, based solely on connectomic information (Fig. 3), we could extract the rules of subcellular innervation preference described for inhibitory axons in the mammalian cortex (17) and if such synaptic target preference could also be found for excitatory axons. We first measured the preference of each axon for innervating dendritic spine heads versus dendritic shafts and other targets (Fig. 4, A and B) because, in the mammalian cortex, most axons of inhibitory interneurons (INs) preferentially innervate the dendrites’ shafts or neuronal somata (17) and most excitatory glutamatergic axons preferentially innervate the spine heads of dendrites (1). The fraction of primary spine synapses per axon (out of all synapses of that axon) accordingly allowed the identification of spine-preferring, likely excitatory axons with at least 50% primary spine innervations (n = 5894 axons) and shaft-preferring, likely inhibitory axons with <20% primary spine innervations (n = 893 axons, or 13.2% of all axons; for exceptions to this rule and control measurements, see the supplementary materials and tables S1 and S2).

Fig. 4 Connectomic definition of axon classes.

(A) Example axons with high (top) and low (bottom) fraction of output synapses made onto dendritic spines. (B) Distribution of spine-targeting fraction over all n = 6979 axons; dashed lines indicate thresholds applied to distinguish non-spine-preferring, likely inhibitory axons (<20% spine innervation, n = 893, 12.8% of all axons) from spine-preferring, mostly excitatory axons (>50% spine innervation, n = 5,894, 84.5%). Diagram shows the definition of primary spine innervations. (C to I) Connectomic definition of axon classes by preferential synaptic innervation of subcellular targets. (C) Two example axons innervating three somata [left, n = 6 synapses onto somata (S) of 14 total, arrows] and an AD (right, n = 2 synapses onto AD of 13 total), respectively. All other cell bodies and ADs are shown in gray. (D) Fraction of synapses onto somata, PDs, ADs, SDs, and AISs for all axons. Binomial probabilities are shown over axons to establish at least one synapse onto the respective target (arrows: magenta, excitatory; black, inhibitory). Black lines indicate the average over axons. (E) Comparison of predicted synapse fraction onto target classes per inhibitory axon on the basis of the binomial probability to innervate the target at least once [gray shading; see arrows in (D)] and measured distribution of synapse fractions onto targets (black lines). (F) Same as (E) but for excitatory axons. (G) Fraction of target-preferring excitatory (Exc.) and inhibitory (Inh.) axons identified using the false detection rate criterion [q = 5 to 30% (39)]. Colored bars indicate the distribution for q = 5% (left) and q = 30% (right). Mixed colors indicate axons specific for both somata and PDs. [(H) and (I)] Second-order innervation preference by target-specific axons; numbers indicate fractional innervation by remaining synapses per axon; colors indicate underfrequent (black) or overfrequent (blue) innervation. Diagonal entries are the fraction of synapses onto the same target (black boxes).

We then determined for each of the subcellular synaptic target classes defined above (Figs. 3 and 4C) the per-synapse innervation probability that would best explain whether an inhibitory axon establishes at least one synapse onto each of these targets. These inhibitory “single-hit” binomial innervation probabilities were 4.2% (somata), 17.8% (PD), 4.9% (SD), 3.3% (AD), and 0.5% (AIS) (Fig. 4D). We then computed the expected distribution of synapses per axon made onto each target class assuming the double-hit, triple-hit, etc., innervation probabilities are the same as the probability to establish at least one synapse onto that target. When comparing these target distributions with the measured distributions of synapses per axon onto each target class (Fig. 4E), we found that inhibitory axons established enhanced preference for somata (p = 2.4 × 10−34, n = 893, one-sided Kolmogorov–Smirnov test), PDs (p = 6.0 × 10−77), ADs (p = 2.5 × 10−4), and to a lesser degree for SDs (p = 1.7 × 10−3, table S1), but no enhanced preference for AISs in L4 (p = 0.648). AISs were synaptically innervated by 0.172 input synapses per micrometer of AIS length, but these innervations were not made by axons with an enhanced preference for AISs, unlike in supragranular and infragranular layers (37).

When performing the same analysis for excitatory axons (Fig. 4F), we found clear target preference for ADs (p = 2.5 × 10−34, Fig. 4F), SDs (p = 7.6 × 10−25), and PDs (p = 1.3 × 10−169). By contrast, thalamocortical (TC) axons [detected using the criteria reported in (38); see fig. S6 and materials and methods] indicated a target preference for PDs (p = 2.5 × 10−31), but not for ADs (p = 0.019) or SDs (p = 0.723). To determine the fraction of inhibitory and excitatory axons that had an unexpectedly high synaptic preference for one (or multiple) of the subcellular target classes, we applied the false detection rate criterion used for the determination of significantly expressed genes [q value (39); see materials and methods] and obtained lower bounds on the fractions of axons in the tissue that preferentially innervate the various subcellular target classes (Fig. 4G; at least 58.0% of inhibitory and 24.4% of excitatory axons). Inhibitory axons (Fig. 4H), but not excitatory axons (Fig. 4I), showed higher-order innervation preferences, indicating that at the level of the dense cortical circuit, synaptic target preferences established by axons were not a continuum but allowed cell-type classification without the need for measurements of neuronal morphology, electrical activity, protein expression, or transcription levels.

Geometric sources of synaptic innervations

Could these local connectivity rules have been derived solely from the geometry of axons and dendrites? We first quantified the overall relation between the spatial distribution of axons and dendrites and the establishment of synapses between them (Fig. 5). One paradigm, originally proposed by Peters (11), states that TC axons entering a certain cortical tissue volume would sample the available cortical dendrites for synaptic innervation according to their relative prevalence in the tissue (1). This model (Fig. 5A) predicted the TC innervation of most cortical dendrites rather well, with the exception of smooth dendrites [an exception reported by White et al. (14)] and the enhanced TC innervation of PDs of layer 4 cells (21). When applied to corticortical excitatory and inhibitory axons (Fig. 5A), we found that this model predicted excitatory innervation of most spiny dendrites rather well, but again failed to predict innervation of SDs and the proximal bias of inhibitory synapses. Because this model [which has been most widely used for circuit inference (2, 12)] implicitly accounts for the density of synapses along the presynaptic axons, it was capable of capturing the increased synapse density of TC axons (Fig. 5A). A simpler variant of the Peters model (1, 15) (Fig. 5B), which uses the density of pre- and postsynaptic path length as basis for the synaptic innervation prediction, failed at predicting the TC innervation but captured the corticocortical innervation of spiny dendrites (Fig. 5B). We then analyzed whether a Peters model normalized for postsynaptic synapse density (Fig. 5C) would better capture synaptic innervation and found that, in fact, the dendritic model was a far better predictor of synaptic innervation (compare Fig. 5, C and B). This indicated that SDs and ADs sampled synaptic input according to the relative path length of the presynaptic axons (Fig. 5C). We then investigated whether a Peters model accounting for pre- and postsynaptic synapse densities would improve the innervation prediction (Fig. 5D). In this model, both the output and the input of cortical excitatory neurites were properly predicted, but the suppressed innervation of SDs and ADs by TC axons and the proximal bias of inhibitory axons was not. Notably, none of the Peters models could account for this proximal bias of inhibitory synapses [Fig. 5D; for other failures of Peters predictions, see, e.g., (3, 6, 8, 9)].

Fig. 5 Contribution of neurite geometry and membrane availability to cortical wiring.

(A to D) Quantitative test of various formulations of Peters’ rule: comparison of actual synaptic innervation to the prediction of synaptic innervation on the basis of the availability of postsynaptic path length in the dataset (A), the product of pre- and postsynaptic path length (B), the sampling of presynaptic partners by their relative prevalence (C), and the product of pre- and postsynaptic synapse density (D). Log likelihood ratios were as follows: –1.1×103 (A), –11×103 (C), and –12×103 (D), all compared with the simple model in (B); p < 10−14 (corrected for degrees of freedom). (E to H) Prediction of single-axon synaptic target preference by distance-dependent postsynaptic surface sampling. (E) Diagram of the surface area of the various subcellular postsynaptic target classes (colors) within a distance rpred from a given axon (black) and example surfaces around two axons within a prediction radius rpred = 5 μm. (F) Surface fraction of target classes around all n = 6979 axons in dependence of rpred around axons. Colors indicate the fraction of synapses of a given axon actually innervating the respective target. (G) Relationship between the surface fraction around all axons and synaptic innervation by these axons for each target (rpred = 10 μm). Black lines indicate linear regression for geometrical innervation prediction. (H) R2 reporting the fraction of synaptic innervation variance [over all axons; see (G)] explained by a multivariate linear innervation model using the available postsynaptic surface area around axons [shaded areas: red, excitatory axons (Exc.); blue, inhibitory axons (Inh.)]; lower end of shades indicates prediction; upper ends indicate correction by the variance contributed by the multinomial sampling of targets along axons; solid lines represent direct prediction of innervation from surface fraction. Dashed lines indicate modeled prediction for a purely geometric forward model at rpred = 10 μm. Insets (right) show sampling-corrected predictive power of excitatory (top) and inhibitory (bottom) axons for the innervation of target classes.

More recently, the Peters model has been investigated for the close proximity between axons and dendrites on the scale of few micrometers (8, 9) and concluded poor (8) or absent (9) geometric predictability of synaptic innervation. We used our larger dense reconstruction to investigate the geometric prediction over a substantially broader spatial scale from 1 to ~30 μm and accounted for inhibitory axons, excitatory axons, and postsynaptic target types (Fig. 5, E to H). We measured whether the postsynaptic membrane surface available within a certain radius rpred around a given axon (Fig. 5E) would be a predictor of synaptic innervation for that given axon. We measured the available membrane surface belonging to the five subcellular target classes around all 6979 axons (Fig. 5F) and used a linear multinomial regression model to predict synaptic innervation from these data (Fig. 5G). Then, we computed the coefficient of determination (R2) reporting the fraction of axonal synaptic innervation variance that could be explained purely based on the geometrical information (Fig. 5H; for details, see the materials and methods). In fact, for small spatial scales of 1 to 5 μm, the membrane surface available around an axon was a rather good predictor of synaptic innervation from excitatory axons (range, 16 to 90%, Fig. 5H; less so for inhibitory axons: range, 23 to 79%).

Would this imply that axonal and dendritic proximity at the single-axon level can be used to infer synaptic connectivity in the cortex (13)? We found that for the spatial alignment scales that can be achieved in light-microscopy–based neuron reconstructions from multiple animals (10 to 20 μm), predictability dropped substantially (Fig. 5H), making circuit inference by an emulation of growth processes based on light-microscopically aligned data (13, 40) implausible.

Subcellular synapse placement

We used our dense reconstruction to study the spatial distribution of synapses along somata and dendrites in the cortical neuropil. The density of TC synapses had a substantial dependence on cortex depth (Fig. 6, A to D): the absolute density of TC synapses in the volume increased by ~93% over 50 μm cortex depth (Fig. 6, A and B); the TC excitatory synapse fraction TC/(TC+CC) (where CC is cortico-cortical) increased by 82.6%, corresponding to an absolute increase in the TC synapse fraction of 5.8% per 50 μm cortex depth (Fig. 6D). This gradient was consistent with light-microscopic analyses of TC synapses indicating a decrease of TC synapse density from lower to upper L4 (41). Neither the inhibitory nor the corticocortical excitatory synapse densities showed a comparable spatial profile (Fig. 6C).

Fig. 6 Gradient of TC synapse density in L4 and ensuing variability of synaptic input composition in L4 neurons.

(A to D) Distribution of TC synapses within the L4 dataset (A): gradient along the cortical axis (B), which is absent for inhibitory (yellow) or CC (blue) synapses (C). (D) Resulting gradient in TC synapse fraction [increase by 83% from 7.0 to 12.8% (+5.8%) within 50 μm along the cortical axis; line fit, p < 1.1 × 10−12, n = 134,537 synapses]. (E to G) Analysis of the variability of TC input onto the primary dendrites of neurons possibly resulting from the TC synapse gradient (D): example reconstructions (E) aligned to the somata; (F) fraction of excitatory input synapses originating from TC axons evaluated for each primary dendrite, plotted according to the direction of the dendrite relative to cortical axis (–1, aligned toward pia; +1, aligned toward WM). TC input fraction [TC/(TC+CC)] of each dendrite compared with the TC input fraction of its entire parent neuron (ratios shown). (G) Summary analysis of relation between dendrite direction and relative TC input fraction showing that the TC input fraction is determined by the dendrites’ orientation relative to the cortex axis (1.28-fold higher relative TC fraction for downward- than upward-pointing dendrites, n = 183, p = 0.026, two-sided t test for dendrites with a normalized absolute projection >0.5; bars correspond to ranges –1 to –0.5; –0.5 to 0.5; and 0.5 to 1). (H to K) Enhanced TC synaptic input (red spheres) is correlated to reduced inhibitory input from AD-preferring inhibitory axons (purple spheres and arrows in H) at the level of single dendrites (r = –0.24, p = 0.0095, n = 183, Pearson’s correlation after Bonferroni’s correction) and for neurons [(J), r = –0.27, p = 0.01, n = 84], but not soma-preferring inhibitory axons [green in (H) and (K), r = 0.08, p = 0.49, n = 84].

How is the synaptic TC gradient mapped onto the dendrites of L4 neurons along the cortex axis (Fig. 6, E to G)? One possibility is that the TC synapse gradient is used to enhance the variability of synaptic input composition between different primary dendrites of the L4 neurons such that a neuron’s dendrites pointing upward toward the pia would sample relatively less TC input than dendrites pointing toward the white matter. Alternatively, mechanisms to establish synaptic target preference (such as those reported in Fig. 4) could be used to counterbalance this synaptic gradient and equilibrate the synaptic input fractions on the differently oriented dendrites. Our analysis showed that, in fact, even at the level of single primary dendrites, TC input fractions were 1.28-fold higher for dendrites pointing upward toward the cortical surface versus downward toward the white matter (Fig. 6, F and G; TC input fractions of each dendrite were corrected for the entire neuron’s TC input fraction; for this analysis, see the materials and methods). We then investigated whether this differential composition of the excitatory inputs is accompanied by different compositions of the inhibitory input synapses (Fig. 6, H to L). We found that the fraction of TC input to a neuron’s dendrites was anticorrelated to the fraction of inhibitory synapses that originated from AD-preferring inhibitory axons (Figs. 6I and 4), both at the level of the input to L4 neurons and at the level of single primary dendrites of L4 neurons (Fig. 6, I and J). The effect was absent for all other synapse classes, most notably the soma-preferring inhibitory axons (Fig. 6K; see discussion).

Connectomic mapping of the plasticity-consistent circuit fraction

The concept of Hebbian plasticity, thought to be at the core of experience-dependent changes of synaptic weights in the brain, makes predictions about the temporal evolution of synaptic weights in multiple synaptic contacts between the same pre- and postsynaptic neurons (AADD joint synapses; Fig. 7, A and B): Because Hebbian synaptic plasticity is dependent on the electrical activity of the pre- and postsynaptic neurons, which in a first approximation can be assumed to be similar at joint synapses, long-term potentiation (LTP) predicts joint synapses to become stronger and relatively more similar in weight (especially if synaptic weight saturates) and long-term depression (LTD) predicts joint synapses to become weaker and relatively more dissimilar in weight (but more similar if synaptic weights saturate; Fig. 7A). Models of LTP and LTD thus make particular predictions about the temporal evolution of joint synaptic weights, and the mapping of synaptic weights and synaptic weight similarity in the connectome allows the quantification of upper bounds on the fraction of the circuit that can have undergone such particular patterns of weight change before the connectomic experiment (we denote those synapse pairs for which such patterns of weight change occurred to a sufficient degree as “having undergone LTP/LTD”; see discussion).

Fig. 7 Connectomic mapping of the plasticity-consistent circuit fraction.

(A) Hebbian LTP makes predictions about the temporal evolution of synaptic size and size similarity in AADD synapse pairs (green; insets show example model trajectories of synapse pairs exposed to LTP with and without weight saturation), yielding a region in the size-similarity plane (right) where synaptic pairs that have undergone LTP are predicted to be found (colors in right panel as in temporal plots on the left). For Hebbian LTD, pairs of synapses behave accordingly only if synaptic size saturates at low values (red). Arrows indicate trajectories of synapse pairs with randomly drawn initial size that undergo LTP with (dark green) or without (light green) weight saturation; LTD with (red) and without (pink and yellow indicate linear and exponential decay, respectively) weight saturation. (B) Example AADD synapse pair (arrows) onto dendritic spines between the same axon (blue) and same dendrite (red). Direct links to datasets are as follows: https://wklink.org/3356 (synapse 1) and https://wklink.org/6145 (synapse 2). (C) Frequency of joint synapse pairs in the dataset (n = 5290 spine–synapse pairs, shaded, analyzed here). (D) ASI as a representative measure of synapse weight (42, 43), dataset link https://wklink.org/5780 (E) Distribution of mean synaptic size and synaptic size similarity for all pairs of AADD synapses from excitatory axons; each dot corresponds to one synapse pair. Isolines indicate statistical regions defined in (F). (F) Map of the relation between synaptic size and synaptic size similarity in AADD pairs, reported as the difference of (E) to random synapse pairs (fig. S7, C and D). Isolines indicate significance levels (p = 0.05 and 0.005 for outer and inner isolines, respectively) outlining overfrequency of synapse pairs that are similar in size and large (upper area) and similar in size and small (lower area). (G and H) Analysis of AADd and AaDD synapse pairs that would indicate a contribution of cell-type-dependent connection size differences. No oversimilarity can be found in these cases (H). (I) Analysis as in (E) and (F) but for TC connections. Note upper bound of 16% of connections consistent with stabilized LTP. (J) Summary of fraction of synapse pairs that resided in the regions identified in (F) and (I) as upper bounds (for the interaction between the two upper bounds, see the supplementary materials). Numbers indicate the ranges for different significance thresholds [see (F) and (I) and fig. S7, C and D]. (K) Analysis as in (G) and (H) but for CC-to-L4 neuron connections only, refuting subtypes of CC connections as the source of the observed oversimilarity (see fig. S7, A and B). Image width is 2 μm in (B) and (D).

We set out to leverage our large connectomic dataset (n = 5290 excitatory joint synaptic pairs onto spines; Fig. 7C) to map the relation between synaptic size and synaptic size similarity in joint synapse pairs [Fig. 7E; for visualization, the figure reports relative synaptic size dissimilarity on the x-axis; for the utilization of the axon–spine interface area (Fig. 7D) as an indicator of synaptic weight (42, 43)]. These data would allow us to determine upper bounds on the plasticity-consistent fraction of the circuit beyond the previous finding that in joint synapse pairs, synaptic size is more similar than for randomly shuffled synapse pairs (9, 29, 30, 44).

Synaptic size similarity in joint synapse pairs showed a broad distribution (Fig. 7E). When comparing this distribution with the synaptic size and synaptic size similarity distribution obtained from a random assignment of the same synapses into “random pairs” (Fig. 7F and fig. S7, C and D), we observed that the population of oversimilar synapse pairs (Fig. 7F) was split into a region of oversimilar and large synapses (mean synaptic size 0.23 to 1.19 μm2; 16 to 20% of all joint synapse pairs are found in this region; the above-random synapse pairs constitute 3.6 to 3.9% of all joint synapse pairs; see fig. S7, C and D, and the materials and methods for details of the region definition and statistics), and oversimilar and small synapses (mean synaptic size 0.06 to 0.2 μm2; 15 to 19% of all joint synapse pairs were found in this region; 3.0 to 3.4% of all joint synapse pairs were above random in this region). With this information, we obtained upper bounds on the fraction of the circuit that can have undergone LTP and LTD with weight saturation (compare Fig. 7, F and A).

To what degree was the observed synaptic weight similarity a result of subtypes of neurons establishing differently sized synapses? Although the quantification of the upper bounds of the plasticity-consistent circuit fraction would remain unaffected, we could use this more detailed analysis to understand whether the plasticity-consistent circuit fractions were specific to types of neuronal connections.

First, we considered the possibility that certain presynaptic cell types made consistently larger or consistently smaller synapses (Fig. 7G). In this case, the distribution of synaptic weight similarity for same-axon different-dendrite (AADd) synapse pairs would also show a bias toward more similarly sized synapses. However, we found no such evidence (Fig. 7H), excluding cell-type-specific synapse size of either presynaptic (axonal) or postsynaptic (dendritic) origin as the cause of the observed oversimilar synapse pairs.

Next, we separated those connections established by TC axons from those made by the remaining excitatory (i.e., CC) axons (Fig. 7I and fig. S7, A and B). We found an excess of oversimilar synapse pairs in the TC connections as well, with 7 to 16% of pairs found in a region of overly similar and large synapses (i.e., an upper bound of 16% on LTP). The region of overly similar and small synapse pairs, however, only comprised 2 to 7% of joint synapse pairs. This remaining number of overly similar small synapse pairs could in fact be induced by the overly similar large synapse pairs (see the supplementary materials). At 28 days of age, ~3 weeks after the proposed critical period during which LTP can be induced in TC connections (45, 46), a fraction of up to 16% of joint synapse pairs was still consistent with previous episodes of LTP that led to stabilized potentiated synapse pairs at dendritic spines (47, 48), but 84% were not.

Repeating these analyses for other combinations of pre- and postsynaptic neurite types (Fig. 7J), we found upper bounds for LTP and LTD of ~10 to 20%. For each of these subtype-specific connections, we could then again analyze whether any purely presynaptic or purely postsynaptic subtype within the already type-selected connections (corresponding to squares in the table of Fig. 7J) could be the cause of the observed synapse similarity. For example, the connections from corticortical axons onto spiny L4 neurons (49) showed no evidence for presynaptic axonal subtypes yielding oversimilar synapses (Fig. 7K; for additional controls of these findings, see the supplementary materials; fig. S7, A and B; and table S1).

Together, these results provided a first quantitative upper bound on the fraction of the circuit consistent with previous episodes of saturated Hebbian synaptic plasticity leading to strengthening or weakening of synapses (a “connectomic fingerprint” of the maximum possible plasticity fraction of the circuit) and excluded obvious cell-type-based connection strength differences as the origin of these observations. Because these results were obtained from brains of untrained animals and were not the result of electrical or other stimulation (“plasticity induction”), these data may represent an unbiased screening of upper bounds of plasticity traces in local cortical circuits, for which the dense connectomic mapping was essential.

Discussion

Using FocusEM, we obtained the first dense circuit reconstruction from the mammalian cerebral cortex at a scale that allowed the analysis of axonal patterns of subcellular innervation, ~300 times larger than previous dense reconstructions from cortex (9). Inhibitory axon types preferentially innervating certain postsynaptic subcellular compartments could be defined solely on the basis of connectomic information (Figs. 3 and 4). In addition to inhibitory axons, a fraction of excitatory axons also exhibited such subcellular innervation preferences (Fig. 4). The geometrical arrangement of axons and dendrites explained only a moderate fraction of synaptic innervation, revoking coarse random models of cortical wiring (Fig. 5). A substantial TC synapse gradient in L4 gave rise to an enhanced heterogeneity of synaptic input composition at the level of single cortical dendrites (Fig. 6), which was accompanied by a reduced innervation from AD-preferring inhibitory inputs. The consistency of synapse size between pairs of axons and dendrites signified fractions of the circuit consistent with saturated synaptic plasticity, placing an upper bound on the “learned” fraction of the circuit (Fig. 7). FocusEM allowed the dense mapping of circuits in the cerebral cortex at a throughput that enables connectomic screening.

Synaptic input composition along L4 dendrites

Our finding of a covariation of enhanced TC inputs to L4 excitatory cells with reduced direct inhibitory input from AD-preferring INs (Fig. 6, H to K) could be interpreted in the context of a disinhibitory circuit described previously (50, 51). Taking into account the preferential targeting of ADs and of soma-preferring parvalbumin (PV)–positive INs by somatostatin (SST)–positive INs, this could imply that SST-IN–based disinhibition can enhance TC input by silencing perisomatic PV inputs recruited by feedforward inhibition (52) and concomitantly reducing the direct inhibitory component from SST INs. In any case, this finding of per-dendrite input variation points to a circuit configuration in which TC input variability is enhanced between neurons of the same excitatory type in cortical layer 4, and furthermore provides evidence for a per-dendrite synaptic input composition of enhanced heterogeneity.

Connectomic traces of plasticity

We interpreted the joint synapse data (Fig. 7) in terms of upper bounds of synapse pairs that could have undergone certain models of plasticity. Although this analysis detects those synapse pairs that were exposed to saturating plasticity (i.e., the possible plasticity event led to a final weight state of both synapses), an alternative interpretation is a dynamic circuit in which at any given point in time, only a fraction of synapses has expressed saturated plasticity, whereas other (or all) synapses are in the process of undergoing plastic changes. We expect that more elaborate plasticity models of entire circuits will also make testable predictions that are accessible by connectomic snapshot experiments as shown here.

Outlook

The presented methods and results open the path to the connectomic screening of nervous tissue from various cortices, layers, species, developmental stages, sensory experiences, and disease conditions. The fact that even a small piece of mammalian cortical neuropil contains a high density of relevant information so rich as to allow the extraction of possible connectomic signatures of the “learnedness” of the circuit makes this approach a promising endeavor for the study of the structural setup of mammalian nervous systems.

Materials and Methods

Animal experiments

A wild-type (C57BL/6) male mouse was transcardially perfused at postnatal day 28 under isoflurane anesthesia using a solution of 2.5% paraformaldehyde and 1.25% glutaraldehyde (pH 7.4) following the protocol in (53). All procedures followed the animal experiment regulations of the Max Planck Society and were approved by the local animal welfare authorities (Regierungspräsidien Oberbayern and Darmstadt).

Tissue sampling and staining

The brain was removed from the skull after 48 hours of fixation and sliced coronally using a vibratome. Two samples were extracted using a 1-mm biopsy punch (Integra Miltex, Plainsboro, NJ) from a 1-mm-thick slice 5 mm from the front of the brain targeted to layer 4 in the somatosensory cortex of the right hemisphere. The corresponding tissue from the left hemisphere was further sliced into 70-μm-thick slices followed by cytochrome oxidase staining, indicating the location of the coronal slice to be in barrel cortex.

Next, the extracted tissue was stained as described previously (53). Briefly, the tissue was immersed in a reduced osmium tetroxide solution (2% OsO4, 0.15 M CB, 2.5 M KFeCN), followed by a 1% thiocarbohydrazide step and a 2% OsO4 step for amplification. After an overnight wash, the sample was further incubated with 1.5% uranyl acetate solution and a 0.02 M lead(II) nitrate solution. The sample was dehydrated with propylene oxide and EtOH, embedded in Epon Hard (Serva Electrophoresis GmbH, Germany), and hardened for 48 hours at 60°C.

3D electron microscopy experiment

The embedded sample was placed on an aluminum stub and trimmed such that the tissue was directly exposed on all four sides of the sample. The sides of the sample were covered with gold in a sputter coater (Leica Microsystems, Wetzlar, Germany). Then, the sample was placed into an SBEM setup [(31), Magellan scanning electron microscope, FEI Company, Hillsboro, OR, equipped with a custom-built microtome courtesy of W. Denk]. The sample was oriented so that the radial cortex axis was in the cutting plane. The transition between L4 and L5A was identified in overview electron microscopy (EM) images by the sudden drop in soma density between the two layers (Fig. 1C). A region of size 96 × 64 μm2 within L4 was selected for imaging using a 3 × 3 image mosaic, a pixel size of 11.24 × 11.24 nm2, an image acquisition rate of 10 MHz, a nominal beam current of 3.2 nA (thus a nominal electron dose of 15.8 e/nm2), an acceleration voltage of 2.5 kV, and a nominal cutting thickness of 28 nm. The effective data rate, including overhead time spent during motor movements for cutting and tiling, was 0.9 MB/s. A total of 3420 image planes were acquired, yielding 194 GB of data.

Image alignment

After 3D EM dataset acquisition, all images were inspected manually and marked for imaging artifacts caused by debris present on the sample surface during imaging. Images with debris artifacts were replaced by the images at the same mosaic position from the previous or subsequent plane. First, rigid translation-only alignment was performed based on the procedures in (53). The following modifications were applied. When shift vectors were obtained that yielded offsets of >100 pixels, these errors were iteratively corrected by manually reducing the weight of the corresponding entry in the least-squares relaxation by a factor of 1000 until the highest remaining residual error was <10 pixels. Shift calculation of subsequent images in cutting direction was found to be the most reliable measurement and was therefore weighted 3-fold in the weighted least-squares relaxation. The resulting shift vectors were applied (shift by integer voxel numbers) and the 3D image data were written in KNOSSOS format (34, 36). For further improvement, subimage alignment was applied (see the supplementary materials and methods).

Methods description for software code

All routines described in the following are available as software at https://gitlab.mpcdf.mpg.de/connectomics/L4dense, which is the relevant reference for the exact sequence of processing steps applied. The following descriptions and the more detailed ones in the supplementary materials and methods are aimed at pointing to the key algorithmic steps rather than enumerating all detailed computations.

Workflow for dense circuit reconstruction

The workflow for volume reconstruction of the acquired 3D EM volume (Fig. 2 and fig. S1) was as follows. We first detected blood vessels and cell bodies using automated heuristics, and then processed the remaining image volume using machine-learning-based image segmentation [CNN and watershed as described in SegEM (32)]. The result of this processing was 15 million volume segments corresponding to pieces of axons, dendrites, and somata (volume: 0.0295 ± 0.3846 μm3, mean ± SD). We then constructed the neighborhood graph between all these volume segments and computed the properties of interfaces between directly adjacent volume segments. On the basis of these features, we trained a connectivity classifier (ConnectEM; Fig. 2, A and B) to determine whether two segments should be connected (along an axon or a dendrite or a glial cell) or if they should be disconnected. Using the SynEM classifier (33), we determined whether an interface between two disconnected processes corresponded to a chemical synapse and, if so, which was the presynaptic and which was the postsynaptic neurite segment (see below for more details). We furthermore trained a set of classifiers (TypeEM; Fig. 2C) to compute for each volume segment the probability of being part of an axon, a dendrite, a spine head, or a glia cell (precision and recall were 91.8 and 92.9% for axons, 95.3 and 90.7% for dendrites, 97.2 and 85.9% for astrocytes, and 92.6 and 94.4% for spine heads, respectively; see table S2).

Cell body–based neuron reconstruction

We next reconstructed those neurons that had their cell bodies in the tissue volume (Fig. 1, E and F, cell gallery in movie S1; n = 125 cell bodies; of these, 97 were neuronal and of these 97, 89 were reconstructed with dendrites in the dataset). For this, we used a set of simple growth rules for automatically connecting neurite pieces on the basis of the segment-to-segment neighborhood graph and the connectivity and neurite type classifiers (fig. S1, “automated agglomeration”; see the supplementary materials and methods). As a result, we obtained fully automated reconstructions of the neuron’s soma and dendritic processes. With a minimal additional manual correction investment of 9.7 hours for 89 cells (54.5 mm dendritic and 2.1 mm axonal path length), the dendritic shafts of these neurons could be reconstructed without merge errors, but there were 37 remaining split errors, at 87.3% dendritic length recall (table S2). This reconstruction efficiency compares favorably to recent reports of automated segmentation of neurons in 3D EM data from the bird brain obtained at ~2-fold higher imaging resolution (54), which reports soma-based neuron reconstruction at an error rate of beyond 100 errors per 66 mm dendritic shafts at lower (68%) dendritic length recall with a similar resource investment (see the supplementary materials and methods).

In addition to the dendritic shafts, the dendritic spines constitute a major fraction of the dendritic path length in cortical neuropil (Fig. 1G). Using our spine head classifier (part of the TypeEM classifiers; Fig. 2C), we found 415,797 spine heads in the tissue volume, which is a density of 0.784 per μm3 (0.98 per μm3 of neuropil, when excluding somata and blood vessels). To connect these to the corresponding dendritic shafts, we trained a spine neck continuity algorithm that was able to automatically attach 58.9% of these spines (evaluated in the center of the dataset at least 10 μm from the dataset border), yielding a dendritic spine density of 0.672 per μm dendritic shaft length [comparable to spine densities in the bird brain (55)]. However, in mammalian cerebral cortex, the density of spines along dendrites is even higher (at least 1 per μm dendritic shaft length). The remaining spine heads were then attached to their dendritic shafts by seeding manual reconstructions at the spine heads and asking annotators to continue along the spine necks to the dendritic shafts. This annotation was performed in the “orthogonal mode” configuration of webKnossos (34), in which the annotator viewed three orthogonal image planes to decide where to continue the respective spine neck [as in KNOSSOS (36)]. The annotation of all remaining spine necks consumed an additional 900 hours of human work for the attachment of 98,221 spines, resulting in a final overall spine density of 0.959 per μm dendritic shaft length.

Dense tissue reconstruction

The reconstruction of neurons starting from their cell bodies, however, was not the main challenge. Rather, the remaining processes, axons and dendrites not connected to a cell body within the dataset and densely packed in the tissue constitute ~97% of the total neuronal path length in this volume of cortex (Fig. 1G). To reconstruct this vast majority of neurites (Fig. 1H), we first used our connectivity and neurite type classifiers (ConnectEM and TypeEM, respectively; Fig. 2) to combine neurite pieces into larger dendritic and axonal agglomerates (“automated agglomeration,” fig. S1 and supplementary materials and methods). Then, we took those agglomerates that had a length of at least 5 μm (n = 74,074 axon agglomerates), detected their endings that were not at the dataset border, and directed focused human annotation to these endings (“queries,” Fig. 2, F and G).

For human annotation, we used an egocentric directed 3D image data view (“flight mode” in webKnossos), which we had previously found to provide maximized human reconstruction speed along axons and dendrites in cortex (34). Here, however, instead of asking human annotators to reconstruct entire dendrites or axons, we only queried their judgment at the endings of automatically reconstructed neurite parts. To make these queries efficient, we made three additions to webKnossos: (i) we oriented the user along the estimated direction of the neurite at its ending, reducing the time the user needs to orient within the 3D brain tissue; (ii) we dynamically stopped the user’s flight along the axon or dendrite whenever another of the already reconstructed neurite agglomerates had been reached; and (iii) we preloaded the next query while the user was annotating (Fig. 2, F and G). Movie S2 illustrates this annotation process for cases of splits and mergers, respectively. Note that the user was able to switch quickly to the next query and, based on its 3D orientation, spent little time orienting in the tissue at the new location. With this, the average user interaction time was 21.3 ± 36.1 s per query, corresponding to an average of 5.5 ± 8.8 μm traveled per query. In total, 242,271 axon-ending queries consumed 1978 paid-out work hours (i.e., including all overheads, 29.4 s per query).

However, we had to account for a second kind of reconstruction error, so-called mergers, which can originate from the original segmentation, the agglomeration procedure, or erroneous flight paths from human queries (Fig. 2H). To detect such mergers, we started with the notion that most of these merger locations will yield a peculiar geometrical arrangement of a 4-fold neurite intersection once all neurite breaks have been corrected (Fig. 2H, “chiasma”). Because such chiasmatic configurations occur rarely in branching neurites, we directed human focused annotation to these locations. First, we automatically detected these chiasmatic locations using a simple heuristic to detect locations at which axon-centered spheres intersected more than three times with the axon [Fig. 2H, n = 55,161 chiasmata; for approaches to detect such locations by machine learning, see (56, 57)]. Then, we positioned the user queries at a certain distance from the chiasma location pointing inward (Fig. 2H) and used a set of case distinctions to query a given chiasma until its configuration had been resolved (see the supplementary materials and methods for details). Chiasma annotation consumed an additional 1132 work hours [note that the detection of endings and chiasmata was iterated eight times for axons (see the supplementary materials and methods) and that, in a final step, we also detected and queried 3-fold neurite configurations to remove remaining mergers].

Synapse detection, types of postsynaptic targets, and connectome reconstruction

Given the reconstructed pre- and postsynaptic neurites in the tissue volume, we then went on to extract their connectome. For this, we used SynEM (33) to detect synapses between the axonal presynaptic processes and the postsynaptic neurites (Fig. 2E).

We trained a dedicated interface classifier for nonspine synapses using training data containing only shaft and soma synapses (Figs. S2 to S5; see the supplementary material and methods). This classifier also used four additional texture filters compared with SynEM in (33), which originated from the voxelwise predictions of a multiclass CNN trained on synaptic junctions, vesicle clouds, mitochondria, and a background class (Fig. 2E).

Because we were interested in analyzing the subcellular specificity of neuronal innervation, we had to also classify which of the postsynaptic membranes belonged to cell bodies; to classify spiny dendrites as belonging to excitatory cells and smooth dendrites as belonging to INs; and to detect AISs and those dendrites that were likely ADs of neurons located in deeper cortical layers. We developed semiautomated heuristics to detect these subcellular compartments (Fig. 3, A to D; see the supplementary materials and methods for details).

Definition of excitatory and inhibitory axons

We used the fraction of primary spine synapses per axon (out of all synapses of that axon, only axons with at least 5 μm path length and at least 10 synapses were analyzed), which had a peak at ~80% (Fig. 4, A and B), to identify spine-preferring, likely excitatory axons with at least 50% primary spine innervations. Similarly, we identified shaft-preferring, likely inhibitory axons with <20% primary spine innervations. Together, this yielded 6449 axons with clear shaft or spine preferences. For the remaining n = 528 axons with primary spine innervations >20% and <50%, we first wanted to exclude remaining mergers between excitatory and inhibitory axons (that would yield intermediate spine innervation rates) and split these axons at possible merger locations (at least 3-fold intersections). Of these, 338 now had at least 10 synapses and spine innervation rates <20% or >50%. The remaining n = 192 axons (2.75% of all axons with at least 10 synapses) were not included in the following analyses. This together yielded n = 5894 excitatory and n = 893 inhibitory axons in our data. For additional controls, see the supplementary materials.

TC axons were defined following parameters described previously (38) (see the supplementary materials and methods).

Analysis of subcellular synaptic target preference

First, we assumed that all synapses of a given axon class have the same probability to innervate a particular postsynaptic target class (as above). We then inferred this single-hit innervation rate for all combinations of presynaptic axon classes and postsynaptic target classes by determining the probability that best explains whether an axon innervated the target class under a binomial model. The optimized binomial model was then used together with the measured number of synapses of each axon to calculate the expected distribution of target innervation rates. A one-sided Kolmogorov–Smirnov test was used to search for the existence of a subpopulation with an increased target innervation rate. To identify those axons that innervated a given target class beyond chance (Fig. 4G), we computed the probability p(t)meas,i,k of finding at least the measured fraction of synapses onto target t for each axon i from axon class k. The p values were also calculated for the expected distribution of target innervation rates and combined with p(t)meas,i,k to estimate the p-value threshold (t)k at which the false discovery rate q (39) crosses 20%. Eighty percent of the axons with p(t)meas,i,k < p̂(t)k innervate target t with a rate above the single-hit innervation probability and are thus considered to be t preferring.

For the analysis of second-order innervation preference (Fig. 4, H and I), we reported the fraction of synapses onto target τ by t-preferring axons of class k after removal of synapses onto t. This innervation rate was compared against the fraction of synapses onto target τ by all axons of class k.

Geometrical predictability analysis

Peters’ rule (1) stipulates that synapses between classes of axons and dendrites are established in proportion to the prevalence of these classes. One variant of Peters’ rule considered (Fig. 5B) makes the prediction that the fraction of synapses from axon class A onto target class T is the product of pA and qT, where pA is the proportion of axonal path length made up by class A, and qT is the proportion of dendritic path length (excluding spines) made up by class T. The measured synapse fractions were compared against the predictions by calculating the ratio of observed to predicted synapse fractions.

Other formulations evaluate these predictions independently for each axon class (Fig. 5A) or each dendrite class (Fig. 5C).

Finally, to assess the effect of incorporating explicit knowledge about the synapse densities of different axon and dendrite classes, a fourth variant of Peters’ rule (Fig. 5D) was considered in which the predicted synapse fraction from axon class A onto target class T is the product of pA and qT, the overall fractions of synapses originating from A and innervating T, respectively.

How much additional information about the neuropil composition around an axon helps to predict its postsynaptic targets was assessed as follows. For each axon, we determined the total surface area of the target classes that were contained within the cylinder of radius rpred around the axon (Fig. 5E) and compared it with the actually innervated target fraction of each axon (Fig. 5, E and F). We then analyzed the correlation between the availability of the target surfaces and the established synapses on these target classes (Fig. 5G).

We then computed R2 using the following model. For all axons of a given type, we used the fraction of target innervations and fractional surface availabilities in a surround of radius rpred to find the optimal multivariate linear regression parameters. To estimate best-case geometric predictability, we then calculated R2 as 1 minus the ratio between the squared residuals of the regression and the synaptic variance on the same axons used for parameter optimization. Here, we corrected for the variance introduced by the finite number of synapses per axon: we used the axons’ fractional surface availabilities within rpred and absolute synapse numbers to calculate the expected binomial variance and subtracted it from the squared residuals.

This analysis made several assumptions that were in favor of a geometrical explanation of synaptic innervation [therefore, the conclusions about a minimal predictability (Fig. 5H) are still upper-bound estimates]. It was assumed that the number of synapses for a given axon was already known; in most settings, only average synapse rates are known for a given circuit. It also assumed that a precise knowledge of the axonal trajectory and the surrounding target surface fractions were available; again, this is usually only available as an average on the scale of rpred of several tens of micrometers.

To relax the assumption of complete knowledge about target availabilities, we repeated the above R2 analysis for a model in which the predicted fractional innervation of a target is the fractional surface availability of that target.

The computational routine used can be found at https://gitlab.mpcdf.mpg.de/connectomics/L4dense in +connectEM/+Connectome/plotGeometricPredictability.m.

Synapse-size consistency analysis

To determine the consistency of primary spine synapses between a given axon–dendrite pair, we calculated the axon–spine interface area (ASI) (42, 43) of a synapse as the total contact area between the corresponding axon and spine head agglomerates. For axon–dendrite pairs connected by exactly two primary spine synapses, we then calculated the coefficient of variation (CV) of the ASI areas by CV = 21/2(ASI1 – ASI2)/(ASI1 + ASI2), where ASI1 and ASI2 are the larger and smaller of the two ASI areas, respectively. To avoid false same-axon, same-dendrite (AADD) pairs caused by remaining merge errors in the axon reconstruction, this analysis was performed only after splitting axons at their branch points. The measured distribution of CV values was compared against the CV values obtained by randomly drawing pairs from all AADD synapses and against the CV values of observed AADd synapse pairs and pairs from different axons onto the same dendrite (AaDD) and from different axons onto different dendrites (AaDd; Fig. 7H). To test whether AADD primary spine synapse pairs are more similar in size than pairs in the control conditions, a one-sided Kolmogorov–Smirnov test was used. We calculated the decimal logarithm of the average ASI area (in square micrometers) and the CV of the ASI areas of each synapse pair to map the size-similarity plane (Fig. 7, F and I). The kernel density estimate of the observed distribution was compared against the distribution expected from random pairs (5000 Monte Carlo samples; fig. S7C) to identify statistically significantly overrepresented regions. Contour lines show the intersection of the significance regions for p-value thresholds of 0.5% and 5% (Fig. 7, E, F, I, and fig. S7, C and D), with the convex hull around the set of all data points. The fraction of data points contained within a contour was used as the upper bound on the fraction of connections consistent with saturated Hebbian plasticity (Fig. 7J).

Statistical methods

The following statistical tests were performed (in order of presentation in the figures):

The existence of axon subpopulation with unexpectedly high synapse rate onto a given target class was tested using the one-sided Kolmogorov–Smirnov test (Fig. 4, E and F). Axons belonging to a given target-preference class were identified on the basis of the false detection rate criterion [q = 20% (39)] (Fig. 4G).

The degree to which synaptic variance is explainable by geometry-based models was evaluated using R2 (Fig. 5H). Binomial variance was corrected for by subtracting the surface fraction-based expected binomial variance from the squared residuals.

F tests were used to evaluate synaptic gradients as function of cortical depth (Fig. 6, B and D) or dendritic orientation (Fig. 6, F and G). For correlation of the TC input fraction with other synaptic input fractions along dendrites, the inhibitory input fraction and seven target-preferential inhibitory and excitatory synapse types were tested. AD-preferring inhibitory synapses were the only ones with significant and substantial correlation (Pearson’s correlation after Bonferroni’s correction for n = 8 multiple tests). The correlation was also significant at the soma level (Pearson’s correlation). Both correlations were also significant using Spearman’s rank correlation.

The four variants of Peters’ rule (Fig. 5, A to D) were compared using a likelihood-ratio test based on the following multinomial model. It was assumed that the pre- and postsynaptic classes of each synapse in the connectome were sampled either after the path-length fractions of these classes (pA and qT) or after the product of the path length and a class-specific likelihood-maximizing relative synapse density. Wilks’ theorem was used to compute the corresponding p values.

To test whether the axon–spine interface areas of a given spine–synapse pair configuration were more similar than randomly sampled pairs, a one-sided Kolmogorov–Smirnov test was used (Fig. 7, H and K).

Supplementary Materials

science.sciencemag.org/content/366/6469/eaay3134/suppl/DC1

Materials and Methods

Supplementary Text

Figs. S1 to S7

Tables S1 and S2

Captions for Movies S1 and S2

References (6072)

References and Notes

Acknowledgments: We thank Y. Hua and K. M. Harris for discussions; A. Gour, Y. Hua, P. Laserstein, and H. Schmidt for discussing unpublished observations; B. Cowgill for experimental contributions; E. Eulig, R. Hesse, C. Schramm, and M. Zecevic for code contributions; C. Guggenberger for compute infrastructure; Scalable Minds for implementation of the fast-swapping user queries in webKnossos; and annotators (see the supplementary materials) for reconstruction work. Funding: Funding was provided by the Max Planck Society. Author contributions: M.H. conceived, initiated, and supervised the project; K.M.B. performed experiments; P.H. provided analysis methods; A.M., M.B., K.M.B., and B.S. developed and performed analyses with contributions by all authors; and A.M. and M.H. wrote the manuscript with contributions by all authors. Competing interests: The authors declare no competing financial interests. Data and materials availability: The raw EM data are available for browsing at https://wklink.org/9276 (see data links in the figure legends). All data and code are publicly available at https://L4dense2019.brain.mpg.de and in updated form at https://gitlab.mpcdf.mpg.de/connectomics/L4dense. Less than 2% of the 3D EM image dataset using the initial image alignment was previously published (3234, 58).
View Abstract

Stay Connected to Science

Navigate This Article