We consider synapse detection as a classification of interfaces between neuronal processes as synaptic or non-synaptic (Figure 2a; see also Mishchenko et al., 2010, Kreshuk et al., 2015, Huang et al., 2016). This approach relies on a volume segmentation of the neuropil sufficient to provide locally continuous neurite pieces (such as provided by SegEM, Berning et al., 2015, for SBEM data of mammalian cortex), for which the contact interfaces can be evaluated.

Figure 2

Download asset Open asset

Figure 2—source data 1

Source data for plot in panel 2d.

https://doi.org/10.7554/eLife.26414.006

Download elife-26414-fig2-data1-v2.xlsx

The unique features of synapses are distributed asymmetrically around the synaptic interface: presynaptically, vesicle pools extend into the presynaptic terminal over at least 100–200 nm; postsynaptically, the PSD has a width of about 20–30 nm. To account for this surround information our classifier considers the subvolumes adjacent to the neurite interface explicitly and separately, unlike previous approaches (Kreshuk et al., 2015; Huang et al., 2016), up to distances of 40, 80, and 160 nm from the interface, restricted to the two segments in question (Figure 2b; the interface itself was considered as an additional subvolume). We then compute a set of 11 texture features (Table 1, this includes the raw data as one feature), and derive 9 simple aggregate statistics over the texture features within the 7 subvolumes. In addition to previously used texture features (Kreshuk et al., 2011, Table 1), we use the local standard deviation, an intensity-variance filter and local entropy to account for the low-variance (‘empty’) postsynaptic spine volume and presynaptic vesicle clouds, respectively (see Figure 2c for filter output examples and Figure 2d for filter distributions at an example synaptic and non-synaptic interface). The ‘sphere average’ feature was intended to provide information about mitochondria, which often impose as false positive synaptic interfaces when adjacent to a plasma membrane. Furthermore, we employ 5 shape features calculated for the border subvolume and the two subvolumes extending 160 nm into the pre- and postsynaptic processes, respectively. Together, the feature vector for classification had 3224 entries for each interface (Table 1).

We developed and tested SynEM on a dataset from layer 4 (L4) of mouse primary somatosensory cortex (S1) acquired using SBEM (dataset 2012-09-28_ex145_07x2, Boergens et al., unpublished; the dataset was also used in developing SegEM, Berning et al., 2015). The dataset had a size of 93 × 60 × 93 µm³ imaged at a voxel size of 11.24 × 11.24 × 28 nm³. The dataset was first volume segmented (SegEM, Berning et al., 2015, Figure 2a, see Figure 2e for a SynEM workflow diagram). Then, all interfaces between all pairs of volume segments were determined, and the respective subvolumes were defined. Next, the texture features were computed on the entire dataset and aggregated as described above. Finally, the shape features were computed. Then, the SynEM classifier was implemented to output a synapse score for each interface and each of the two possible pre-to-postsynaptic directions (Figure 3a–c). The SynEM score was then thresholded to obtain an automated classification of interfaces into synaptic / non-synaptic (θ in Figure 3a). Since the SynEM scores for the two possible synaptic directions at a given neurite-to-neurite interface were rather disjunct in the range of relevant thresholds, we used the larger of the two scores for classification (Figure 3b; θ_s and θ_nn refer to the SynEM thresholds optimized for single synapse or neuron-to-neuron connectome reconstruction, respectively, see below).

Figure 3 with 4 supplements see all

Download asset Open asset

Figure 3—source data 1

Source data for plots in panels 3a, 3b, 3d, 3e, 3f.

https://doi.org/10.7554/eLife.26414.009

Download elife-26414-fig3-data1-v2.xlsx

We obtained labels for SynEM training and validation by presenting raw data volumes of (1.6 × 1.6 × 0.7–1.7) µm³ that surrounded the segment interfaces to trained student annotators (using a custom-made annotation interface in Matlab, Figure 3—figure supplement 1). The raw data were rotated such that the interface was most vertically oriented in the image plane presented to the annotators; the two interfacing neurite segments were colored transparently for identification (this could be switched off by the annotators when inspecting the synapse, see Materials and methods for details). Annotators were asked to categorize the presented interface as either non-synaptic, pre-to-postsynaptic, or post-to-presynaptic (Figure 3c, Figure 3—figure supplement 1). The synaptic labels were then verified by an expert neuroscientist. A total of 75,383 interfaces (1858 synaptic, 73,525 non-synaptic) were annotated in image volumes drawn from 40 locations within the entire EM dataset (Figure 3—figure supplement 2). About 80% of the labels (1467 synaptic, 61,619 non-synaptic) were used for training, the remaining were used for validation.

Initially, we interpreted the annotator’s labels in an undirected fashion: irrespective of synapse direction, the label was interpreted as synaptic (and non-synaptic otherwise, Figure 3c, ‘Undir.’). We then augmented the training data by including mirror-reflected copies of the originally presented synapses, maintaining the labels as synaptic (irrespective of synapse direction) and non-synaptic (Figure 3c, ‘Augmented’). Finally, we changed the labels of the augmented training data to reflect the direction of synaptic contact: only synapses in one direction were labeled as synaptic, and non-synaptic in the inverse direction (Figure 3c, ‘Directed’).

Figure 3d shows the effect of the choice of features, aggregate statistics, classifier parameters and label types on SynEM precision and recall. Our initial classifier used the texture features from Kreshuk et al. (2011) with minor modifications and in addition the number of voxels of the interface and the two interfacing neurite segmentation objects (restricted to 160 nm distance from the interface) as a first shape feature (Table 1). This classifier provided only about 70% precision and recall (Figure 3d). We then extended the feature space by adding more texture features capturing local image statistics (Table 1) and shape features. In particular, we added filters capturing local image variance in an attempt to represent the ‘empty’ appearance of postsynaptic spines, and the presynaptic vesicle clouds imposing as high-frequency high-variance features in the EM images. Also, we added more subvolumes over which features were aggregated (see Figure 2b), increasing the dimension of the feature space from 603 to 3224. Together with additional aggregate statistics, the classifier reached about 75% precision and recall. A substantial improvement was obtained by switching from an ensemble of decision-stumps (one-level decision tree) trained by AdaBoostM1 (Freund and Schapire, 1997) as classifier to decision stumps trained by LogitBoost (Friedman et al., 2000). In addition, the directed label set proved to be superior. Together, these improvements yielded a precision and recall of 87% and 86% on the validation set (Figure 3d).

We then evaluated the best classifier from the validation set (Figure 3d, ‘Direct and Logit’) on a separate test set. This test set was a dense volume annotation of all synapses in a randomly positioned region containing dense neuropil of size 5.8 × 5.8 × 7.2 µm³ from the L4 mouse cortex dataset. All synapses were identified by two experts, which included the reconstruction of all local axons, and validated once more by another expert on a subset of synapses. In total, the test set contained 235 synapses and 20319 non-synaptic interfaces. SynEM automatically classified these at 88% precision and recall (Figure 3e, F1 score of 0.883). Since the majority of synapses in the cortex are made onto spines we also evaluated SynEM on all spine synapses in the test set (n = 204 of 235 synapses, 87%, Figure 3e). On these, SynEM performed even better, yielding 94% precision and 89% recall. (Figure 3e, F1 score of 0.914).

We next compared SynEM to previously published synapse detection methods (Figure 3f, Mishchenko et al., 2010; Kreshuk et al., 2011, 2014; Becker et al., 2012; Roncal et al., 2015; Dorkenwald et al., 2017). Other published methods were either already shown to be inferior to one of these approaches (Perez et al., 2014; Márquez Neila et al., 2016) or developed for specific subtypes of synapses, only (Jagadeesh et al., 2014; Plaza et al., 2014; Huang et al., 2016); these were therefore not included in the comparison. SynEM outperforms the state-of-the-art methods when applied to our SBEM data acquired at 3537 nm³ voxel size (Figure 3f, Figure 3—figure supplement 3). In addition, we applied SynEM to a published 3D EM dataset acquired at more than 10-fold smaller voxel size (3 × 3 × 30 = 270 nm³) using automated tape-collecting ultramicrotome-SEM imaging (ATUM, Kasthuri et al., 2015). SynEM also outperforms the method developed for this data (VesicleCNN, Roncal et al., 2015; Figure 3f and Figure 3—figure supplement 4), indicating that SynEM is applicable to EM data of various modalities and resolution.

It should furthermore be noted that for connectomics, in addition to the detection of the location of a synapse, the two neuronal partners that form the synapse and the direction of the synapse have to be determined. The performance of the published methods as reported in Figure 3f only include the synapse detection step. Interestingly, the recently published method (Dorkenwald et al., 2017) reported that the additional detection of the synaptic partners yielded a drop of performance of 2% precision and 9% recall (F1 score decreased by about 5% from 0.906 to 0.849) compared to synapse detection alone (Figure 3f, see Dorkenwald et al., 2017). This indicates that the actual performance of this method on our data would be lower when including partner detection. SynEM, because of the explicit classification of directed neurite interfaces, in contrast, explicitly provides synapse detection, partner detection and synapse directionality in one classification step.

Figure 4a shows examples of correct and incorrect SynEM classification results (evaluated at θ_s). Typical sources of errors are vesicle clouds close to membranes that target nearby neurites (Figure 4a, FP), Mitochondria in the pre- and/or postsynaptic process, very small vesicle clouds and/or small PSDs (Figure 4a, FN), and remaining SegEM segmentation errors. To estimate the effect of segmentation errors on SynEM performance, we investigated all false positive and false negative detections in the test set and checked for the local volume segmentation quality. We found that, in fact, 26 of the 28 FNs and 22 of the 27 FPs were at locations with a SegEM error in proximity. Correcting these errors also corrected the SynEM errors in 22 of 48 (46%) of the cases. This indicates that further improvement of volume segmentation can yield an even further reduction of the remaining errors in SynEM-based automated synapse detection.

Figure 4

Download asset Open asset

Figure 4—source data 1

Source data for plot in panel 4b.

https://doi.org/10.7554/eLife.26414.017

Download elife-26414-fig4-data1-v2.xlsx

We then asked which of the SynEM features had highest classification power, and whether the newly introduced texture and shape features contributed to classification. Boosted decision-stump classifiers allow the ranking of features according to their classification importance (Figure 4b). 378 out of 3224 features contributed to classification (leaving out the remaining features did not reduce accuracy). The 10 features with highest discriminative power (Table 2) in fact contained two of the added texture filters (int-var and local entropy) and a shape feature. The three most distinctive subvolumes (Figure 4b) were the large presynaptic subvolume, the border and the small postsynaptic subvolume. This suggests that the asymmetry in pre- vs. postsynaptic aggregation volumes in fact contributed to classification performance, with a focus on the presynaptic vesicle cloud and the postsynaptic density.

Finally, SynEM is sufficiently computationally efficient to be applied to large connectomics datasets. The total runtime on the 384592 μm³ dataset was 2.6 hr on a mid-size computational cluster (480 CPU cores, 16 GB RAM per core). This would imply a runtime of 279.9 days for a large 1 mm³ dataset, which is comparable to the time required for current segmentation methods, but much faster than the currently required human annotation time (10⁵ to 10⁶ hr, Figure 1c). Note that SynEM was not yet optimized for computational speed (plain matlab code, see git repository posted at https://gitlab.mpcdf.mpg.de/connectomics/SynEM).

We so far evaluated SynEM on the basis of the detection performance of single synaptic interfaces. Since we are interested in measuring the connectivity matrices of large-scale mammalian cortical circuits (connectomes) we obtained a statistical estimate of connectome error rates based on synapse detection error rates. We assume that the goal is a binary connectome containing the information whether pairs of neurons are connected or not. Automated synapse detection provides us with weighted connectomes reporting the number of synapses between neurons, from which we can obtain binary connectomes by considering all neuron pairs with at least γ_nn synapses as connected (Figure 5a). Synaptic connections between neurons in the mammalian cerebral cortex have been found to be established via multiple synapses per neuron pair (Figure 5b, Feldmeyer et al., 1999, 2002, 2006; Frick et al., 2008; Markram et al., 1997, range 1–8 synapses per connection, mean 4.3 ± 1.4 for excitatory connections, Supplementary file 2). The effect of synapse recall R_s on recall of neuron-to-neuron connectivity R_nn can be estimated (Figure 5c) for each threshold γ_nn given the distribution of the number of synapses per connected neuron pair n_syn. For connectomes in which neuron pairs with at least one detected synapse are considered as connected (γ_nn = 1), a neuron-to-neuron connectivity recall R_nn of 97% can be achieved with a synapse detection recall R_s of 65.1% (Figure 5c, black arrow) if synapse detection is independent between multiple synapses of the same neuron pair. SynEM achieves 99.4% synapse detection precision P_s at this recall (Figure 3e).

Figure 5 with 2 supplements see all

Download asset Open asset

Figure 5—source data 1

Source data for plots in panels 5b, 5c, 5d, 5e.

https://doi.org/10.7554/eLife.26414.020

Download elife-26414-fig5-data1-v2.xlsx

The resulting precision of neuron-to-neuron connectivity P_nn then follows from the total number of synapses in the connectome N_syn = N² × c_r×<n_syn>, with c_r the pairwise connectivity rate, about 20% for local excitatory connections in cortex (Feldmeyer et al., 1999), <n_syn> the mean number of synapses per connection (4.3 ± 1.4, Figure 5b), and N² the size of the connectome. A fraction R_s of these synapses is detected (true positive detections, TPs). The number of false positive (FP) synapse detections was deduced from TP and the synapse precision P_s as FP=TP×(1-P_s)/P_s, yielding R_s×N_syn×(1-P_s)/P_s false positive synapse detections. These we assumed to be distributed randomly on the connectome and estimated how often at least γ_nn synapses fell into a previously empty connectome entry. These we considered as false positive connectome entries, whose rate yields the binary connectome precision P_nn (see Materials and methods for details of the calculation). At R_nn of 97.1%, SynEM yields a neuron-to-neuron connection precision P_nn of 98.5% (Figure 5d, black arrow, Figure 5e; note that this result is stable against varying underlying connectivity rates c_re = 5%..30%, see indicated ranges in Figure 5e).

For the treatment of inhibitory connections, we followed the notion that synapse detection performance could be optimized by restricting classifications to interfaces established by inhibitory axons (as we had analogously seen for restricting analysis to spine synapses above, Figure 3e). For this, we evaluated SynEM on a test set of inhibitory axons for which we classified all neurite contacts of these axons (171 synapses, 9430 interfaces). While the precision and recall for single inhibitory synapses is lower than for excitatory ones (75% recall, 82% precision, Figure 5—figure supplement 1, SynEM⁽ⁱ⁾_s), the higher number of synapses per connected cell pair (n⁽ⁱ⁾_syn is on average about 6, Supplementary file 3, Gupta et al. (2000); Markram et al. (2004); Koelbl et al. (2015); Hoffmann et al. (2015)) still yields substantial neuron-to-neuron precision and recall also for inhibitory connectomes (98% recall, 97% precision, Figure 5e, Figure 5—figure supplement 1, SynEM⁽ⁱ⁾_nn; this result is stable against varying underlying inhibitory connectivity rates c_ri = 20%..80%, see ranges indicated in Figure 5e). Error rates of less than 3% for missed connections and for wrongly detected connections are well below the noise of synaptic connectivity so far found in real biological circuits (e.g., Helmstaedter et al., 2013; Bartol et al., 2015), and thus likely sufficient for a large range of studies involving the mapping of cortical connectomes.

In summary, SynEM provides fully automated detection of synapses, their synaptic partner neurites and synapse direction for binary mammalian connectomes up to 97% precision and recall, a range which was previously prohibitively expensive to attain in large-scale volumes by existing methods (Figure 5e, Figure 5—figure supplement 2).

We applied SynEM to a sparse local cortical connectome between 104 axons and 100 postsynaptic processes in the dataset from L4 of mouse cortex (Figure 6a, neurites were reconstructed using webKnossos (Boergens et al., 2017) and SegEM as previously reported (Berning et al., 2015)). We first detected all contacts and calculated the total contact area between each pair of pre- and postsynaptic processes (‘contactome’, Figure 6b). We then classified all contacts using SynEM (at the classification threshold θ_nn (Table 3) yielding 98.5% precision and 97.1% recall for excitatory neuron-to-neuron connections and 97.3% precision and 98.5% recall for inhibitory neuron-to-neuron connections) to obtain the weighted connectome C_w (Figure 6c). The detected synapses were clustered when they were closer than 1500 nm for a given neurite pair. This allowed us to concatenate large synapses with multiple active zones or multiple contributing SegEM segments into one (Figure 6—figure supplement 1). To obtain the binary connectome we thresholded the weighted connectome at γ_nn = 1 for excitatory and at γ_nn = 2 for inhibitory neuron-to-neuron connections (Figure 6d). The resulting connectome contained 880 synapses distributed over 536 connections.

Figure 6 with 1 supplement see all

Download asset Open asset

Figure 6—source data 1

Source data for plots in panels 6b, 6c, 6d.

https://doi.org/10.7554/eLife.26414.026

Download elife-26414-fig6-data1-v2.xlsx

Finally, to check whether SynEM-detected synapses matched previous reports on synapse frequency and size, we applied SynEM to half of the entire cortex dataset used for this study (i.e. a volume of 192296 µm³). SynEM detected 195644 synapses, i.e. a synapse density of 1.02 synapses per µm³, consistent with previous reports (Merchán-Pérez et al., 2014).

We then measured the size of the axon-spine interface of SynEM detected synapses in the test set (Figure 7a,b). We find an average axon-spine interface size of 0.263 ± 0.206 µm² (mean ± s.d.; range 0.033–1.189 µm²; n = 181), consistent with previous reports (de Vivo et al., 2017: (SW) 0.297 ± 0.297 µm² (p=0.518, two-sample two-tailed t-test on the natural logarithm of the axon-spine interface size), (EW) 0.284 ± 0.275 µm² (p=0.826, two-sample two-tailed t-test on the natural logarithm of the axon-spine interface size). This indicates that, first, synapse detection in our lower-resolution SBEM data (in-plane image resolution about 11 nm, section thickness about 26–30 nm) yields similar synapse size distributions as in the higher-resolution data in de Vivo et al. (2017) (in-plane image resolution 5.9 nm; section thickness about 50 nm) and secondly, that SynEM-based synapse detection has no obvious bias towards larger synapses.

Figure 7

Download asset Open asset

Figure 7—source data 1

Source data for plots in panels 7a, 7b.

https://doi.org/10.7554/eLife.26414.030

Download elife-26414-fig7-data1-v2.xlsx

Neuropil composition (Figure 1b) was considered as follows: Neuron density of 157,500 per mm³ (White and Peters, 1993), axon path length density of <strike>4</strike> km per mm³ and dendrite path length density of 1 km per mm³ (Braitenberg and Schüz, 1998), spine density of about 1 per µm dendritic shaft length, with about 2 µm spine neck length per spine (thus twice the dendritic path length), synapse density of 1 synapse per µm³ (Merchán-Pérez et al., 2014) and bouton density of 0.1–0.25 per µm axonal path length (Braitenberg and Schüz, 1998). Annotation times were estimated as 200–400 hr per mm path length for contouring, 3.7–7.2 h/mm path length for skeletonization (Helmstaedter et al., 2011, 2013; Berning et al., 2015), 0.6 h/mm for flight-mode annotation (Boergens et al., 2017), 0.1 h/µm³ for synapse annotation by volume search (estimated form the test set annotation) and an effective interaction time of 60 s per identified bouton for axon-based synapse search. All annotation times refer to single-annotator work hours, redundancy may be increased to reduce error rates in neurite and synapse annotation in these estimates (see Helmstaedter et al., 2011).

SynEM was developed and tested on a SBEM dataset from layer 4 of mouse primary somatosensory cortex (dataset 2012-09-28_ex145_07x2, K.M.B. and M.H., unpublished data, see also Berning et al., 2015). Tissue was conventionally en-bloc stained (Briggman et al., 2011) with standard chemical fixation yielding compressed extracellular space (compare to Pallotto et al., 2015).

The image dataset was volume segmented using the SegEM algorithm (Berning et al., 2015). Briefly, SegEM was run using CNN 20130516T204040_8,3 and segmentation parameters as follows: r_se = 0; θ_ms = 50; θ_hm = 0.39; (see last column inTable 2 in Berning et al., 2015). For training data generation, a different voxel threshold for watershed marker size θ_ms = 10 was used. For test set and local connectome calculation the SegEM parameter set optimized for whole cell segmentations was used (r_se = 0; θ_ms = 50; θ_hm = 0.25, see Table 2 in Berning et al., 2015).

Interfaces between a given pair of segments in the SegEM volume segmentation were extracted by collecting all voxels from the one-voxel boundary of the segmentation for which that pair of segments was present in the boundary’s 26-neighborhood. Then, all interface voxels for a given pair of segments were linked by connected components, and if multiple connected components were created, these were treated as separate interfaces. Interface components with a size of 150 voxels or less were discarded.

To define the subvolumes around an interface used for feature aggregation (Figure 2b), we collected all voxels that were at a maximal distance of 40, 80 and 160 nm from any interface voxel and that were within either of the two adjacent segments of the interface. The interface itself was also considered as a subvolume yielding a total of 7 subvolumes for each interface.

Eleven 3-dimensional image filters with one to 15 instances each (Table 1) were calculated as follows and aggregated over the 7 subvolumes of an interface using 9 summary statistics, yielding 3224 features per directed interface. Image filters were applied to cuboids of size 548 × 548 × 268 voxels, each, which overlapped by 72,72 and 24 voxels in x,y and z dimension, respectively, to ensure that all interface subvolumes were fully contained in the filter output.

Gaussian filters were defined by evaluating the unnormalized 3d Gaussian density function

at integer coordinates (x, y, z) $\in$ U = {-f_x,-f_x-1, … f_x} x {-f_y,-f_y-1, … f_y} x {-f_z,-f_z-1, … f_z} for a given standard deviation σ = (σ_x, σ_y, σ_z) and a filter size f = (f_x, f_y, f_z) and normalizing the resulting filter by the sum over all its elements

First and second order derivatives of Gaussian filters were defined as

and analogously for the other partial derivatives. Normalization of g_σ and evaluation of derivatives of Gaussian filters was done on U as described above. Filters were applied to the raw data I via convolution (denoted by ${\displaystyle \ast }$) and we defined the image’s Gaussian derivatives as

and analogously for the other partial derivatives.

Gaussian smoothing was defined as ${\displaystyle I\ast g_{\sigma }}$.

Difference of Gaussians was defined as (${\displaystyle I\ast g_{\sigma }-I\ast g_{k\sigma }}$), where the standard deviation of the second Gaussian filter is multiplied element-wise by the scalar k.

Gaussian gradient magnitude was defined as

Laplacian of Gaussian was defined as

Structure tensor S was defined as a matrix of products of first order Gaussian derivatives, convolved with an additional Gaussian filter (window function) g_σw:

and analogously for the other dimensions, with standard deviation σ_D of the image’s Gauss derivatives. Since S is symmetric, only the diagonal and upper diagonal entries were determined, the eigenvalues were calculated and sorted by increasing absolute value.

The Hessian matrix was defined as the matrix of second order Gaussian derivatives:

and analogously for the other dimensions. Eigenvalues were calculated as described for the Structure tensor.

The local entropy feature was defined as

where p(L) is the relative frequency of the voxel intensity in the range {0, …, 255} in a given neighborhood U of the voxel of interest (calculated using the entropyfilt function in MATLAB).

Local standard deviation for a voxel at location (x, y, z) was defined by

for the neighborhood U of location (x, y, z) with |U| number of elements and calculated using MATLABs stdfilt function.

Sphere average was defined as the mean raw data intensity for a spherical neighborhood U_r with radius r around the voxel of interest, with

where Z³ is the 3 dimensional integer grid; x,y,z are voxel indices; z anisotropy was approximately corrected.

The intensity/variance feature for voxel location (x, y, z) was defined as

for the neighborhood U of location (x, y, z).

The set of parameters for which filters were calculated is summarized in Table 1.

11 shape features were calculated for the border subvolume and the two 160 nm-restricted subvolumes, respectively. For this, the center locations (midpoints) of all voxels of a subvolume were considered. Shape features were defined as follows: The number of voxel feature was defined as the total number of voxels in the subvolumes. The voxel based diameter was defined as the diameter of a sphere with the same volume as the number of voxels of the subvolumes. Principal axes lengths were defined as the three eigenvalues of the covariance matrix of the respective voxel locations. Principal axes product was defined as the scalar product of the first principal components of the voxel locations in the two 160 nm-restricted subvolumes. Voxel based convex hull was defined as the number of voxels within the convex hull of the respective subvolume voxels (calculated using the convhull function in MATLAB).

Interfaces were annotated by three trained undergraduate students using a custom-written GUI (in MATLAB, Figure 3—figure supplement 1). A total of 40 non-overlapping rectangular volumes within the center 86 × 52 × 86 μm³ of the dataset were selected (39 sized 5.6 × 5.6 × 5.6 μm³ each and one of size 9.6 × 6.8 × 8.3 μm³). Then, all interfaces within these volumes were extracted as described above. Interfaces with a center of mass less than 1.124 µm from the volume border were not considered. For each interface, a raw data volume of size (1.6 × 1.6 × 0.7–1.7) μm³, centered on the center of mass of the interface voxel locations was presented to the annotator. When the center of mass was not part of the interface, the closest interface voxel was used. The raw data were rotated such that the second and third principal components of the interface voxel locations (restricted to a local surround of 15 x 15 × 7 voxels around the center of mass of the interface) defined the horizontal and vertical axes of the displayed images. First, the image plane located at the center of mass of the interface was shown. The two segmentation objects were transparently overlaid (Figure 3—figure supplement 1) in separate colors (the annotator could switch the labels off for better visibility of raw data). The annotator had the option to play a video of the image stack or to manually browse through the images. The default video playback started at the first image. An additional video playback mode started at the center of mass of the interface, briefly transparently highlighted the segmentation objects of the interface, and then played the image stack in reverse order to the first plane and from there to the last plane. In most cases, this already yielded a decision. In addition, annotators had the option to switch between the three orthogonal reslices of the raw data at the interface location (Figure 3—figure supplement 1). The annotators were asked to label the presented interfaces as non-synaptic or synaptic. For the synaptic label, they were asked to indicate the direction of the synapse (see Figure 3—figure supplement 1). In addition to the annotation label interfaces could be marked as ‘undecided’. Interfaces were annotated by one annotator each. The interfaces marked as undecided were validated by an expert neuroscientist. In addition, all synapse annotations were validated by an expert neuroscientist, and a subset of non-synaptic interfaces was cross-checked. Together, 75,383 interfaces (1858 synaptic, 73,525 non-synaptic) were labeled this way. Of these, the interfaces from eight label volumes (391 synaptic and 11906 non-synaptic interfaces) were used as validation set; the interfaces from the other 32 label volumes were used for training.

The target labels for the undirected, augmented and directed label sets were defined as described in the Results (Figure 3c). We used boosted decision stumps (level-one decision trees) trained by the AdaBoostM1 (Freund and Schapire, 1997) or LogitBoost (Friedman et al., 2000) implementation from the MATLAB Statistical Toolbox (fitensemble). In both cases the learning rate was set to 0.1 and the total number of weak learners to 1500. Misclassification cost for the synaptic class was set to 100. Precision and recall values of classification results were reported with respect to the synaptic class. For validation, the undirected label set was used, irrespective of the label set used in training. If the classifier was trained using the directed label set then the thresholded prediction for both orientations were combined by logical OR.

To obtain an independent test set disjunct from the data used for training and validation, we randomly selected a volume of size 512 × 512 × 256 voxels (5.75 × 5.75 × 7.17 μm³) from the dataset that contained no soma or dominatingly large dendrite. One volume was not used because of unusually severe local image alignment issues which are meanwhile solved for the entire dataset. The test volume had the bounding box [3713, 2817, 129, 4224, 3328, 384] in the dataset. First, the volume was searched for synapses (see Figure 1d) in webKnossos (Boergens et al., 2017) by an expert neuroscientist. Then, all axons in the volume were skeleton-traced using webKnossos. Along the axons, synapses were searched (strategy in Figure 1e) by inspecting vesicle clouds for further potential synapses. Afterwards the expert searched for vesicle clouds not associated with any previously traced axon and applied the same procedure as above. In total, that expert found 335 potential synapses. A second expert neuroscientist used the tracings and synapse annotations from the first expert to search for further synapse locations. The second expert added eight potential synapse locations. All 343 resulting potential synapses were collected and independently assessed by both experts as synaptic or not. The experts labeled 282 potential locations as synaptic, each. Of these, 261 were in agreement. The 42 disagreement locations (21 from each annotator) were re-examined jointly by both experts and validated by a third expert on a subset of all synapses. 18 of the 42 locations were confirmed as synaptic, of which one was just outside the bounding box. Thus, in total, 278 synapses were identified. The precision and recall of the two experts in their independent assessment with respect to this final set of synapses was 93.6%, 94.6% (expert 1) and 97.9%, 98.9% (expert 2), respectively.

Afterwards all shaft synapses were labeled by the first expert and proofread by the second. Subsequently, the synaptic interfaces were voxel-labeled to be compatible with the methods by Becker et al. (2012) and Dorkenwald et al. (2017). This initial test set comprised 278 synapses, of which 36 were labeled as shaft/inhibitory.

Next, all interfaces between pairs of segmentation objects in the test volume were extracted as described above. Then, the synapse labels were assigned to those interfaces whose border voxels had any overlap with one of the 278 voxel-labeled synaptic interfaces. Afterwards, these interface labels were again proof-read by an expert neuroscientist. Finally, interfaces closer than 160 nm from the boundary of the test volume were excluded to ensure that interfaces were fully contained in the test volume. The final test set comprised 235 synapses out of which 31 were labeled as shaft/inhibitory. With this we obtained a high-quality test set providing both voxel-labeled synapses and synapse labels for interfaces, to allow the comparison of different detection methods.

For the calculation of precision and recall, a synapse was considered detected if at least one interface that had overlapped with the synapse was detected by the classifier (TPs); a synapse was considered missed if no overlapping interface of a given synapse was detected (FNs); and a detection was considered false positive (FP) if the corresponding interface did not overlap with any labeled synapse.

The labels for inhibitory-focused synapse detection were generated using skeleton tracings of inhibitory axons. Two expert neuroscientists used these skeleton tracings to independently detect all synapse locations along the axons. Agreeing locations were considered synapses and disagreeing locations were resolved jointly by both annotators. The resulting test set contains 171 synapses. Afterwards, all SegEM segments of the consensus postsynaptic neurite were collected locally at the synapse location. For synapse classification all interfaces in the dataset were considered that contained one SegEM segment located in one of these inhibitory axons. Out of these interfaces all interfaces were labeled synaptic that were between the axon and a segment identified as postsynaptic. The calculation of precision and recall curves was done as for the dense test set (see above) by considering a synapse detected if at least one interface overlapping with it was detected by the classifier (TPs); a synapse was considered missed if no interface of a synapse was detected (FNs); and a detection was considered false positive (FP) if the corresponding interface did not overlap with any labeled synapse.

The approach of Becker et al. (2012) was evaluated using the implementation provided in Ilastik (Sommer et al., 2011). This approach requires voxel labels of synapses. We therefore first created training labels: an expert neuroscientist created sparse voxel labels at interfaces between pre- and postsynaptic processes and twice as many labels for non-synaptic voxels for five cubes of size 3.4 × 3.4 × 3.4 μm³ that were centered in five of the volumes used for training SynEM. Synaptic labels were made for 115 synapses (note that the training set in Becker et al. (2012)) only contained 7–20 synapses). Non-synaptic labels were made for two training cubes first. The non-synaptic labels of the remaining cubes were made in an iterative fashion by first training the classifier on the already created synaptic and non-synaptic voxel labels and then adding annotations specifically for misclassified locations using Ilastik. Eventually, non-synaptic labels in the first two training cubes were extended using the same procedure.

For voxel classification all features proposed in (Becker et al., 2012) and 200 weak learners were used. The classification was done on a tiling of the test set into cubes of size 256 × 256 × 256 voxels (2.9 × 2.9 × 7.2 μm³) with a border of 280 nm around each tile. After classification, the borders were discarded, and tiles were stitched together. The classifier output was thresholded and morphologically closed with a cubic structuring element of three voxels edge length. Then, connected components of the thresholded classifier output with a size of at least 50 voxels were identified. Synapse detection precision and recall rates were determined as follows: A ground truth synapse (from the final test set) was considered detected (TP) if it had at least a single voxel overlap with a predicted component. A ground truth synapse was counted as a false negative detection if it did not overlap with any predicted component (FN). To determine false positive classifications, we evaluated the center of the test volume (shrunk by 160 nm from each side to 484 × 484 × 246 voxels) and counted each predicted component that did not overlap with any of the ground truth synapses as false positive detection (FP). For this last step, we used all ground truth synapses from the initial test set, in favor of the Becker et al. (2012) classifier.

For comparison with (Kreshuk et al., 2014) the same voxel training data as for (Becker et al., 2012) was used. The features provided by Ilastik up to a standard deviation of 5 voxels for the voxel classification step were used. For segmentation of the voxel probability output map the graph cut segmentation algorithm of Ilastik was used with label smoothing ([1, 1, 0.5] voxel standard deviation), a voxel probability threshold of 0.5 and graph cut constant of λ = 0.25. Objects were annotated in five additional cubes of size 3.4 × 3.4 × 3.4 μm³ that were centered in five of the interface training set cubes different from the one used for voxel prediction resulting in 299 labels (101 synaptic, 198 non-synaptic). All object features provided by Ilastik were used for object classification. The evaluation on the test set was done as for (Becker et al., 2012).

For comparison with (Dorkenwald et al., 2017) six of the 32 training cubes used for interface classification with a total volume of 225 μm³ were annotated with voxel labels for synaptic junctions, vesicle clouds and mitochondria. The annotation of vesicle clouds and mitochondria was done using voxel predictions of a convolutional neural network (CNN) trained on mitochondria, vesicle clouds and membranes. The membrane predictions were discarded and the vesicle clouds and mitochondria labels were first proofread by undergraduate students and then twice by an expert neuroscientist. The voxels labels for synaptic junctions were added by an expert neuroscientist based on the identified synapses in the interface training data. Overall 310 synapses were annotated in the training volume. A recursive multi-class CNN was trained on this data with the same architecture and hyperparameter settings as described in (Dorkenwald et al., 2017) using the ElektroNN framework. For the evaluation of synapse detection performance only the synaptic junction output was used. The evaluation on the test set was done as for (Becker et al., 2012) with a connected component threshold of 250 voxels.

The image data, neurite and synapse segmentation from (Kasthuri et al., 2015) hosted on openconnecto.me (kasthuri11cc, kat11segments, kat11synapses) was used (downloaded using the provided scripts at https://github.com/neurodata-arxiv/CAJAL). The segmentation in the bounding box [2432, 7552; 6656, 10112; 769, 1537] (resolution 1) was adapted to have a one-voxel boundary between segments by first morphologically eroding the original segmentation with a 3-voxel cubic structuring element and running the MATLAB watershed function on the distance-transform of the eroded segmentation on a tiling with cubes of size [1024, 1024, 512] voxels. Since the Kasthuri et al. (2015) segmentation in the selected bounding box was not dense, voxels with a segment id of zero in the original segmentation whose neighbors at a maximal distance of 2 voxels (maximum-distance) also all had segment ids zero were set to segment id zero in the adapted segmentation. All segments in the adapted segmentation that were overlapping with a segment in the original segmentation were set to the id of the segment in the original segmentation. The bounding box [2817, 6912; 7041, 10112; 897, 1408] of the resulting segmentation was tiled into non-overlapping cubes of [512, 512, 256] voxels. For all synapses in the synapse segmentation the pre- and postsynaptic segment of the synapse were marked using webKnossos (Boergens et al., 2017) and all interfaces between the corresponding segments at a maximal distance of 750 nm to the synapse centroid that were also overlapping with an object in the synapse segmentation were associated to the corresponding synapse and assigned a unique group id. Only synapses labeled as ‘sure’ in Kasthuri et al. (2015) were evaluated. All interfaces with a center of mass in the region ac3 with the bounding box [5472, 6496; 8712, 9736; 1000, 1256] were used for testing. All interfaces with a center of mass at a distance of at least 1 μm to ac3 were used for training if there was no interface between the same segment ids in the test set. Interfaces between the same segment ids as an interface in the test set were only considered for training if the distance to ac3 was above 2 μm. For feature calculation the standard deviation of Gaussian filters was adapted to the voxel size 6 × 6 × 30 nm of the data (i.e. s in Table 1 was set to 12/6 in x- and y-dimension and 12/30 in z-dimension). The directed label set approach was used for classification. The calculation of precision recall rates was done as described above (‘test set generation and evaluation’).

The neuron-to-neuron connection recall was calculated assuming an empirical distribution p(n) of the number of synapses n between connected excitatory neurons given by published studies (see Supplementary file 2, Feldmeyer et al., 1999, 2002, 2006; Frick et al., 2008; Markram et al., 1997). For inhibitory connections we used a fixed value of 6 synapses (see Supplementary file 3, Koelbl et al., 2015; Hoffmann et al., 2015; Gupta et al., 2000; Markram et al., 2004). We further assumed that the number of retrieved synapses is given by a binomial model with retrieval probability given by the synapse classifier recall R_s on the test set:

Where γ_nn is the threshold on the number of synapses between a neuron pair to consider it as connected (see Figure 5a). This equates to the neuron-to-neuron recall: R_nn = P(k ≥ γ_nn | R_s).

To compute the neuron-to-neuron precision, we first calculated the expected number of false positive synapse detections (FP_s) made by a classifier with precision P_s and recall R_s:

where N_syn is the total number of synapses in a dataset calculated from the average number of synapses per connected neuron pair <n_syn> times the number of connected neuron pairs N_con and c_r is the connectivity ratio given by N_con/N² with N the number of neurons in the connectome.

We then assumed that these false positive synapse detections occur randomly and therefore are assigned to one out of N² possible neuron-to-neuron connections with a frequency FP_s/N².

We then used a Poisson distribution to estimate the number of cases in which at least γ_nn FP_s synapses would occur in a previously zero entry of the connectome, yielding a false positive neuron-to-neuron connection (FP_nn).

Finally, the true positive detections of neuron-to-neuron connections in the connectome TP_nn are given in terms of the neuron-to-neuron connection recall R_nn by

Together, the neuron-to-neuron connection precision P_nn is given by

The connectivity ratio was set to c_r = 0.2 (Feldmeyer et al., 1999) for excitatory and to 0.6 for inhibitory connections (Gibson et al., 1999; Koelbl et al., 2015).

For determining the local connectome (Figure 6) between 104 pre- and 100 postsynaptic processes, we used 104 axonal skeleton tracings (traced at 1 to 5-fold redundancy) and 100 dendrite skeleton tracings. 10 axons were identified as inhibitory and are partially contained in the inhibitory test set. All volume objects which overlapped with any of the skeleton nodes were detected and concatenated to a given neurite volume. Then, all interfaces between pre- and postsynaptic processes were classified by SynEM. The area of each interface was calculated as in (Berning et al., 2015) and the total area of all contacts between all neurite pairs was calculated (Figure 6b). To obtain the weighted connectome C_w (Figure 6c), we applied the SynEM scores threshold θ_nn (Table 3) for the respective presynaptic type (excitatory, inhibitory). Detected synaptic interfaces were clustered using hierarchical clustering (single linkage, distance cutoff 1500 nm) if the interfaces were between the same pre- and postsynaptic objects. To obtain the binary connectome C_bin (Figure 6d) we thresholded the weighted connectome at the connectome threshold γ_nn = 1 for excitatory and γ_nn = 2 for inhibitory connections (Table 3). The overall number of synapses in the dataset was calculated by considering all interfaces above the score threshold for the best single synapse performance (θ_s) as synaptic. To obtain the final synapse count the retrieved synaptic interfaces were clustered using hierarchical clustering with single linkage and a distance cutoff between the centroids of the interfaces of 320.12 nm (this distance cutoff was obtained by optimizing the synapse density prediction on the test set).

For the evaluation of axon-spine interface area (ASI) all spine synapses in the test set were considered for which SynEM had detected at least one overlapping neurite interface (using θ_s for spine synapses, Figure 3e). The ASI of a detected synapse was calculated by summing the area of all interfaces between segmentation objects that overlapped with the synapse. For comparison to ASI distributions obtained at higher imaging resolution in a recent study (spontaneous wake (SW) and enforced wake (EW) conditions reported in Table S1 in de Vivo et al. (2017)), it was assumed that the ASI distributions are lognormal (see Figure 2B in de Vivo et al., 2017). Two-sample two-tailed t-tests were performed for comparing the natural logarithmic values of the SynEM-detected ASI from the test set (log ASI −1.60 ± 0.74, n = 181; mean ± s.d.) with the lognormal distributions for SW and EW from de Vivo et al. (2017), (log ASI −1.56 ± 0.83, n = 839, SW; −1.59 ± 0.81, n = 836, EW; mean ± s.d.), p=0.5175 (SW) and p=0.8258 (EW).

All code used to train and run SynEM are available as source code and also at https://gitlab.mpcdf.mpg.de/connectomics/SynEM under the MIT license. A copy is available at https://github.com/elifesciences-publications/SynEM. To run SynEM, please follow instructions in the readme.md file. Data used to train and evaluate SynEM is available at http://synem.brain.mpg.de.

The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.

We thank Jan Gleixner for first test experiments on synapse detection and fruitful discussions in an early phase of the project, Alessandro Motta for comments on the manuscript, Christian Guggenberger for excellent support with compute infrastructure, Raphael Jacobi, Raphael Kneissl, Athanasia Natalia Marahori, and Anna Satzger for data annotation and Elias Eulig, Robin Hesse, Martin Schmidt, Christian Schramm and Matej Zecevic for data curation. We thank Heiko Wissler and Dalila Rustemovic for support with illustrations.

Animal experimentation: All animal experiments were performed in accordance with the guidelines for the Use of Laboratory Animals of the Max Planck Society and approved by the local authorities Regierungspräsidium Oberbayern, AZ 55.2-1-54-2532.3-103-12.

© 2017, Staffler et al.

This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.