## Abstract

In the cerebral cortex, most synapses are found in the neuropil, but relatively little is known about their 3-dimensional organization. Using an automated dual-beam electron microscope that combines focused ion beam milling and scanning electron microscopy, we have been able to obtain 10 three-dimensional samples with an average volume of 180 µm^{3} from the neuropil of layer III of the young rat somatosensory cortex (hindlimb representation). We have used specific software tools to fully reconstruct 1695 synaptic junctions present in these samples and to accurately quantify the number of synapses per unit volume. These tools also allowed us to determine synapse position and to analyze their spatial distribution using spatial statistical methods. Our results indicate that the distribution of synaptic junctions in the neuropil is nearly random, only constrained by the fact that synapses cannot overlap in space. A theoretical model based on random sequential absorption, which closely reproduces the actual distribution of synapses, is also presented.

## Introduction

In the cerebral cortex, most synapses (90–98%) are established in the neuropil (Alonso-Nanclares et al. 2008), which is made of dendrites, axons, and glial processes. One major issue in cortical circuitry is to ascertain how synapses are distributed and whether synaptic connections are specific or not (DeFelipe et al. 2002). Therefore, to understand the anatomical design principles of the cortical circuits, it is essential to analyze the ultrastructure of all components of the neuropil and, in particular, the number and spatial distribution of synapses. There are 2 major morphological types of synapses: Asymmetric (or Gray's type I) and symmetric (or Gray's type II; Gray 1959; Colonnier 1968). The major sources of asymmetric synapses are spiny neurons (pyramidal and spiny nonpyramidal cells) and extrinsic cortical afferents, whereas the vast majority of symmetric synapses are of intrinsic origin and are formed by the population of aspiny or sparsely spiny interneurons. Since spiny neurons and the major cortical afferent systems are excitatory in function, and virtually all aspiny or sparsely spiny interneurons are γ--aminobutyric (GABA acid)-ergic, it is assumed that axon terminals forming asymmetric synapses are excitatory and those forming symmetric synapses are inhibitory (Ascoli et al. 2008). Furthermore, synaptic size plays an important role in the functional properties of synapses (Schikorski and Stevens 1997; Takumi et al. 1999; Lüscher et al. 2000; Tarusawa et al. 2009). Thus, numerous researchers have been trying to find simple and accurate methods for estimating the distribution, size, and number of synapses. To this end, 2 sampling procedures are currently available: One is based on serial reconstructions and the other on single sections. Clearly, serial reconstruction should be the method of choice for the challenging task of unraveling the extraordinary complexity of the nervous system. Indeed, serial sectioning transmission electron microscopy is a well established and mature technology for obtaining 3-dimensional data from ultrathin sections of brain tissue (Stevens et al. 1980; Harris et al. 2006; Hoffpauir et al. 2007; Mishchenko et al. 2010; Bock et al. 2011). It is based on imaging ribbons of consecutive sections with a conventional transmission electron microscope. However, the major limitation is that obtaining long series of ultrathin sections is extremely time-consuming and difficult, often making it impossible to reconstruct large volumes of tissue. Hence, the recent development of automated electron microscopy techniques represents another crucial step in the study of neuronal circuits (e.g. Briggman and Denk 2006; Knott et al. 2008; Merchán-Pérez, Rodriguez, Alonso-Nanclares, et al. 2009).

In the present study, we use a new dual-beam electron microscope that combines a focused ion beam (FIB) column and a scanning electron microscope (SEM). The FIB column directs a gallium (Ga^{+}) ion beam toward the sample. The collision of Ga^{+} ions with substrate atoms results in the removal of these atoms from the substrate. Since the FIB can be focused and controlled on a nanometer scale, this effect can be used to mill the specimen, removing a thin layer of material. A backscattered electron image is then obtained from the freshly milled surface using the SEM. The milling/imaging processes are automatically repeated in order to obtain a long series of images that represent a 3-dimensional sample of tissue. Using this methodology, we have previously shown that virtually all cortical synapses can be accurately identified as asymmetric (excitatory) or symmetric (inhibitory) synaptic junctions when they are visualized, reconstructed, and quantified from large 3-dimensional tissue samples obtained in an automated manner (Merchán-Pérez, Rodriguez, Alonso-Nanclares, et al. 2009). Since the segmentation, quantification, and analysis of synaptic junctions are a labor-intensive procedure, we used Espina software, a specifically developed tool that greatly facilitates and accelerates these processes (Morales et al. 2011). For the segmentation of synaptic junctions, Espina makes use of the fact that presynaptic and especially postsynaptic densities are electron-dense structures. It uses a gray-level threshold to extract all the voxels that fit the gray levels of the synaptic junction. Asymmetric synaptic junctions have a thicker postsynaptic density, so they can be distinguished from symmetric ones, whose postsynaptic densities are thinner and similar to the presynaptic densities. As a result, each synaptic junction (already tagged as asymmetric or symmetric) is segmented as a flattened irregular volume and its geometrical features, including its spatial position, are determined by the same software.

Points in space—the synaptic junctions in the present case—can be distributed according to 3 basic patterns: Random, clustered, and regular. A random pattern follows the basic reference model of a point process, the so-called complete spatial randomness (CSR) or homogeneous spatial Poisson point process (Seidel 1876). It is characterized by 2 fundamental properties: (1) The number of points falling in a given region follow a Poisson distribution with a constant mean number of points per unit volume; and (2) each point is equally likely to occur at any position and is not affected by the locations of the other points (Illian et al. 2008; Gaetan and Guyon 2009). In a clustered distribution, the points tend to concentrate in groups, leaving other regions of space empty or almost devoid of points. Finally, in a regular or dispersed pattern, every point is located as far as possible from its neighbors, resulting in a lattice-like distribution of points. Although there are no clear-cut limits between these 3 basic patterns, CSR represents a boundary condition between clustered and dispersed spatial processes, so our strategy was first to test for CSR, trying to identify deviations from it and to build a model that accounted for these deviations (Illian et al. 2008; Gaetan and Guyon 2009).

## Materials and Methods

### Tissue Preparation and 3-Dimensional Electron Microscopy

Male Wistar rats sacrificed in postnatal day 14 were used for this study. Animals were administered a lethal intraperitoneal injection of sodium pentobarbital (40 mg/kg) and were intracardially perfused with 2% paraformaldehyde and 2.5% glutaraldehyde in 0.1 M phosphate buffer. The brain was then extracted from the skull and processed for electron microscopy according to a previously described protocol (Merchán-Pérez, Rodriguez, Alonso-Nanclares, et al. 2009). All animals were handled in accordance with the guidelines for animal research set out in the European Community Directive 2010/63/EU, and all procedures were approved by the local ethics committee of the Spanish National Research Council (CSIC).

Three-dimensional brain tissue samples were obtained using a combined FIB/SEM (Crossbeam^{®} Neon40 EsB, Carl Zeiss NTS GmbH, Oberkochen, Germany). This instrument combines a high-resolution field emission SEM column with a focused gallium ion beam, which can mill the sample surface, removing thin layers of material on a nanometer scale. After removing each slice (20 nm of thickness), the milling process was paused, and the freshly exposed surface was imaged with a 1.8-kV acceleration potential using the in-column energy selective backscattered (EsB) electron detector. A 30-mm aperture was used and the retarding potential of the EsB grid was 1500 V. The milling and imaging processes were sequentially repeated, and long series of images were acquired through a fully automated procedure, thus obtaining a stack of images that represented a 3-dimensional sample of the tissue (Merchán-Pérez, Rodriguez, Alonso-Nanclares, et al. 2009). Image resolution in the *xy* plane was 3.7 nm/pixel. Resolution in the *z* axis (section thickness) was 20 nm.

Ten different samples of the neuropil of layer III of the somatosensory cortex were obtained from 3 different animals (Table 1). To select the exact location of the samples, we first obtained plastic semithin sections (2 µm thick) from the block surface and stained them with toluidine blue to identify cortical layers. These sections were then photographed with a light microscope. The last of these light microscope images (corresponding to the section immediately adjacent to the block face) was then collated with SEM photographs of the surface of the block. In this way, it was possible to accurately identify the regions of the neuropil to be studied. All samples selected for FIB milling/SEM imaging were located at mid-depth of layer III. After applying a 3-dimensional unbiased counting frame and correcting for tissue shrinkage (Howard and Reed 2005; Merchán-Pérez, Rodriguez, Alonso-Nanclares, et al. 2009; Morales et al. 2011), the volume of tissue analyzed from each sample ranged from 149.13 to 247.58 µm^{3} (mean 180.33 µm^{3}; see Supplementary Table 1). Synaptic junctions within these volumes were visualized, automatically segmented, and reconstructed in 3 dimensions with Espina software (Morales et al. 2011; Fig. 1). This software also provided data such as the volume of the unbiased counting frame, the number of synaptic junctions inside the frame, the spatial position of the centroids or centers of gravity of the synaptic junctions, and an estimation of their sizes by the Feret's diameter (the diameter of the smallest sphere circumscribing the synaptic junction; Fig. 1 and SupplementaryVideos 1 and 2). As previously described in detail (Merchán-Pérez, Rodriguez, Alonso-Nanclares, et al. 2009), synaptic junctions with either a prominent or thin postsynaptic density were classified as asymmetric or symmetric synaptic junctions, respectively (Supplementary Video 1).

Sample no. and animal identification | Number of synapses per cubic micron | Mean distance to nearest neighbor (nm) ± SD | Mean Feret's diameter of synaptic junctions (nm) ± SD |
---|---|---|---|

1 (W31) | 0.9857 | 519.55 ± 136.35 | 377.19 ± 159.63 |

2 (W31) | 0.6936 | 594.07 ± 192.28 | 462.18 ± 177.52 |

3 (W33) | 0.9279 | 537.43 ± 159.20 | 437.62 ± 168.04 |

4 (W33) | 1.0088 | 537.39 ± 157.70 | 414.22 ± 169.04 |

5 (W33) | 0.9474 | 597.30 ± 174.02 | 466.03 ± 215.91 |

6 (W33) | 0.9399 | 533.21 ± 163.29 | 423.38 ± 169.83 |

7 (W33) | 0.9881 | 487.17 ± 172.30 | 397.29 ± 168.22 |

8 (W35) | 0.7997 | 568.21 ± 178.51 | 427.90 ± 168.15 |

9 (W35) | 1.1267 | 501.38 ± 156.97 | 378.35 ± 166.60 |

10 (W35) | 1.0178 | 523.74 ± 150.36 | 405.43 ± 175.62 |

All samples | 0.9399 | 535.78 ± 166.81 | 417.06 ± 175.97 |

Sample no. and animal identification | Number of synapses per cubic micron | Mean distance to nearest neighbor (nm) ± SD | Mean Feret's diameter of synaptic junctions (nm) ± SD |
---|---|---|---|

1 (W31) | 0.9857 | 519.55 ± 136.35 | 377.19 ± 159.63 |

2 (W31) | 0.6936 | 594.07 ± 192.28 | 462.18 ± 177.52 |

3 (W33) | 0.9279 | 537.43 ± 159.20 | 437.62 ± 168.04 |

4 (W33) | 1.0088 | 537.39 ± 157.70 | 414.22 ± 169.04 |

5 (W33) | 0.9474 | 597.30 ± 174.02 | 466.03 ± 215.91 |

6 (W33) | 0.9399 | 533.21 ± 163.29 | 423.38 ± 169.83 |

7 (W33) | 0.9881 | 487.17 ± 172.30 | 397.29 ± 168.22 |

8 (W35) | 0.7997 | 568.21 ± 178.51 | 427.90 ± 168.15 |

9 (W35) | 1.1267 | 501.38 ± 156.97 | 378.35 ± 166.60 |

10 (W35) | 1.0178 | 523.74 ± 150.36 | 405.43 ± 175.62 |

All samples | 0.9399 | 535.78 ± 166.81 | 417.06 ± 175.97 |

Note: Densities of synapses were calculated by dividing the actual number of synaptic junctions within an unbiased 3-dimensional counting frame by the volume of the counting frame. Distances to the nearest neighbor were calculated between the centroids of the synaptic junctions. The Feret's diameter is the diameter of the smallest sphere circumscribing the synaptic junction; it is used here as an estimate of the size of synaptic junctions.

### Spatial Point Pattern Analysis

To analyze the spatial distribution of synaptic junctions, we compared the actual positions of the centroids of synaptic junctions with 2 theoretical point processes: CSR and random sequential adsorption (RSA). CSR defines a situation where a point is equally likely to occur at any location within the study volume, regardless of the locations of other points (Diggle 2003). The mathematical model underlying CSR is a homogeneous spatial Poisson point process. Thus, the number of points occurring within a finite volume, *V*, is a random variable following a Poisson distribution with mean *λ*|*V*|, where *λ* > 0 and |*V*| is the volume of *V*. The intensity *λ* represents the expected number of points per unit volume. CSR serves as a boundary condition between spatial processes that are more clustered than random and spatial processes that are more regular than random.

For modeling the spatial distribution of synapses, we should take into account that synaptic junctions cannot overlap, and thus the minimum intersynaptic distances (measured between their centroids) must be limited by the sizes of synaptic junctions themselves. A spatial point process to model this situation can be seen as a type of hard-core point process where no points can occur at a distance smaller than a specific minimum distance (Illian et al. 2008). In this respect, an RSA process (Evans 1993) is a kind of hard-core process where the pattern is constructed by placing spheres in 3-dimensional space iteratively and randomly, with their radii following a probability density function. If any newly generated sphere intersects with an existing sphere, the new sphere is rejected and another sphere with a different center and radius is generated. This process is stopped when the required number of spheres is reached (Supplementary Video 3).

In our case, the minimum interpoint distance would depend at least partially on the sizes of synaptic junctions. Thus, we obtained the empirical distribution of the Feret's diameters of synaptic junctions as an estimate of their sizes. We then performed a goodness-of-fit test to find the theoretical probability density distribution that best fitted to the empirical distribution. We found that the best fit using the Kolmogorov–Smirnov test was with a log-normal distribution with parameters *µ* = 5.828 and *σ* = 0.446 (Fig. 2). The estimates of the expectation and variance were *υ* = 375.198 nm and *τ*^{2} = 30,981.351 nm^{2}, respectively (Supplementary Methods).

For each of the 10 different samples and their corresponding CSR and RSA processes, we calculated 3 functions commonly used for spatial point pattern analysis: the *G*, *F*, and *K* functions (Ripley and Kelly 1977; Diggle 1979).

The *G* function, also called the nearest-neighbor distance cumulative distribution function or the event-to-event distribution, is, for a distance *d*, the probability that a typical point separates from its nearest-neighbor a distance of at most *d*. The empirical distribution function *G** of the observed nearest-neighbor distances is used to estimate the *G* function:

*D*is the distance of synapse

_{i}*i*to its nearest-neighbor synapse,

*N*the number of synapses, and

*I*(·) the indicator function which is one if the argument is true and is zero otherwise.

The *F* function, also called the empty space function or the point-to-event distribution, is, for a distance *d*, the probability that the distance of each point (in a regularly spaced grid of *L* points superimposed over the sample) to its nearest synapse centroid is at most *d*. The empirical distribution function *F** of the observed empty space distances on a grid of *L* locations is used to estimate the *F* function:

*D*′

_{i}is the distance of grid point

*i*to its nearest centroid. No edge correction was used for

*G*and

*F*functions.

The *K* function, also called the reduced second moment function or Ripley's function, is, for a distance *d*, the expected number of points within a distance *d* of a typical point of the process divided by the intensity *λ*. An estimation of the *K* function is given by the mean number of points within a sphere of increasing radius *d* centered on each sample point, divided by an estimation of the expected number of points per unit volume. A common formula is the Milles–Lantuejoul–Stoyan–Hanisch translation-corrected estimate (Baddeley et al. 1993), given by:

*V*| is the volume of the counting frame,

*D*″

*is the distance between centroid*

_{ij}*i*and centroid

*j*, and |

*V*| denotes the volume of the 3-dimensional subspace obtained as the intersection between

_{ij}*V*and the (translated) sphere centered at synapse

*i*and radius

*D*″

*.*

_{ij}Spatial analysis of the positions of the centroids of synaptic junctions in experimental and simulated distributions was performed with R Project, specifically with “spatstat” (Baddeley and Turner 2005) and “fitdistplus” (Venables and Ripley 2002; Vose 2008) packages and with SA3D software (Eglen et al. 2008). Additional statistical analyses were performed with SPSS (IBM Corp., New York, USA).

## Results

We obtained 10 independent samples by FIB/SEM technology from regions of the neuropil of layer III of the somatosensory cortex (hindlimb representation) from three 14-day-old rats. From each of the 10 samples of cortical tissue used in this study, we obtained the number, position, and shape of synaptic junctions inside a 3-dimensional unbiased counting frame (Howard and Reed 2005; Fig. 1 and Supplementary Video 2). No cell somata or blood vessels were present within the counting frames. In this way, we directly calculated the number of synapses per unit volume of neuropil. These densities ranged from 0.69 to 1.13 synapses/µm^{3}, with a mean of 0.94 synapses/µm^{3} (Table 1). The spatial positions of the geometric barycenters or centroids of each individual synaptic junction (Fig. 1 and Supplementary Video 2) were recorded for further analysis. As an estimation of the size of synaptic junctions, we also recorded the Feret's diameter of each synaptic junction, which is the diameter of the smallest sphere circumscribing the synaptic junction (Fig. 1, Table 1, and Supplementary Videos 1 and 2). The total tissue volume studied for the 10 samples was 1803.32 µm^{3}, where we found 1695 synapses (1555 asymmetric and 140 symmetric). Additional data on these tissue volumes and the actual counts of asymmetric and symmetric synaptic junctions per sample are shown in Supplementary Table 1.

In order to explore the spatial distribution pattern of synaptic junctions, we first performed a nearest-neighbor analysis on the centroids of the synaptic junctions found within each sampled volume. We calculated the mean distance from each centroid to its nearest neighbor within the counting frame. Centroids that were closer to the boundaries of the counting frame than to any other centroid were excluded from the calculations, since their nearest neighbor could lay outside the counting frame at an unknown distance (Baddeley et al. 1993; Illian et al. 2008). The mean distance to the nearest neighbor ranged from 487.17 to 597.30 nm (Table 1).

To further analyze the spatial distribution of nearest neighbors, we obtained the *G* function (Diggle 1979). This function is estimated using the distances from each centroid to its nearest neighbor, and plotting the fraction of points in the sample having their nearest neighbor at a given distance or less (Fig. 3). We then compared the experimentally observed *G* functions with the theoretical *G* functions that would be obtained if the same number of centroids were randomly distributed according to the CSR or homogeneous Poisson process. In the latter conditions, the function takes a sigmoid shape, as in the examples shown in Figure 3*A*,*C* (red traces) and Supplementary Figures 1–10. In a clustered pattern, the curve would be skewed to the left (shorter interpoint distances than expected), while it would be shifted to the right in a regular pattern (longer interpoint distances than expected).

The experimentally observed *G* functions showed a similar shape in the 10 different samples. Two examples of them are shown in Figure 3 (blue traces) and the rest are shown in Supplementary Figures 1–10. In all cases, when the experimentally observed *G* functions were compared with the *G* functions of the corresponding Poisson processes, the more noticeable differences were found at short interpoint distances, since all experimental samples showed an empty space around centroids where no other centroids were found (arrows in Fig. 3*A*,C, see also Supplementary Figs 1–10). Due to this empty space, the fraction of centroids having their nearest neighbor at short distances was lower than would be expected in a homogeneous Poisson process, while at longer distances the observed curves tended to overlap or slightly exceed the Poisson curves. In general, 20–40% of the centroids had their nearest neighbor at 500 nm or less, while practically all of them had their nearest neighbor at 1000 nm or less (Fig. 3*A*,*C* and Supplementary Figs 1–10).

The existence of this dead space can be explained by the fact that the centroids of synaptic junctions represent the geometric center of a physical volume (the volume of the reconstructed synaptic junction) instead of a dimensionless point, and thus their position cannot be independent from the position of their immediate neighbors. In other words, centroids cannot be too close to each other since the volumes they represent, the synaptic junctions, cannot overlap. Therefore, in order to simulate the distribution of synapses in the neuropil, the empty space around the centroids of synaptic junctions had to be accounted for by a random minimum-spacing distance that would depend, at least partially, on the sizes of synaptic junctions themselves.

To account for this minimum-spacing distance between centroids, we simulated random distributions of points in space using a sequential algorithm that first generates a random point that is the center of a sphere with a certain diameter. The diameters used in the simulations were randomly drawn from the probability distribution of the experimentally observed Feret's diameters (Table 1). This distribution was found to be log-normal (Fig. 2). If the next simulated point with its corresponding Feret's diameter does not overlap with previously generated points, then it is accepted. Otherwise, it is rejected and another random point is generated. The process terminates when the desired number of points has been reached (Supplementary Video 3). This kind of process is known as RSA (Evans 1993).

In order to compare the experimentally observed distribution of centroids with the RSA process, 100 simulated realizations of RSA were generated for each experimental sample, with the same number of points and volumes as the original samples. Fitness between the observed *G* functions and the ones obtained from simulated RSA processes was better than that previously found with homogeneous Poisson processes, although the observed curves tended to be located slightly to the left of the RSA curves (Fig. 3 and Supplementary Figs 1–10).

Further analysis of the spatial distribution of the centroids of synaptic junctions was performed by calculating the *F* and *K* functions (Baddeley et al. 1993; O'Sullivan and Unwin 2002; Gaetan and Guyon 2009). To estimate the *F* function, a regular grid is traced within the 3-dimensional bounding box that contains the centroids, the distances between each grid crossing point and its nearest neighboring centroid are measured, and the cumulative probability of having the nearest centroid at a given distance or less is plotted. In a homogeneous Poisson process, the *F* function takes the same values as the *G* function (Baddeley et al. 1993; Fig. 4 and Supplementary Figs 1–10). Values of the *F* function greater than the Poisson value suggest that there is regularity in the point pattern, while lower values suggest clustering. The *K* function is estimated as the mean number of points within a sphere of increasing radius centered on each sample point (Ripley and Kelly 1977). It is a rapidly growing curve in a homogeneous Poisson process (Fig. 5 and Supplementary Figs 1–10), with lower values in a regular pattern and higher values in a clustered pattern.

In the case of *F* functions, the Poisson and RSA curves were very similar, and the experimentally observed curves either overlapped with RSA curves or tended to be located slightly to the right of them (Fig. 4 and Supplementary Figs 1–10). *K* functions were also very similar for both the Poisson and RSA models. The experimentally observed curves of the *K* function were very similar to the Poisson and RSA curves or were located slightly to the left of them (Fig. 5 and Supplementary Figs 1–10).

In order to test the goodness-of-fit of the experimentally observed functions in the 10 different samples with those for the corresponding Poisson and RSA processes, we performed a battery of 2-sample Kolmogorov–Smirnov tests (Table 2). The tests showed that fitness was better for *F* and *K* functions than for *G* functions. Comparing the actual samples with RSA simulations yielded a clearly better fit than comparing with Poisson distributions. In fact, in the case of RSA simulations, 26 of 30 tests yielded *P*-values higher than 0.950, while for Poisson distributions 17 of 30 tests yielded *P*-values over 0.950 (Table 2). Only for sample 2, the null hypothesis of RSA was rejected for the *G* function (*P* = 0.035). For the homogeneous Poisson process, the null hypothesis was only rejected for the *G* function of sample 10 (*P* = 0.007). In these samples, however, the *G* functions showed no clear deviations toward either a clustered or regular distribution pattern (Supplementary Figs 2 and 10).

Sample no. | Two-sample Kolmogorov–Smirnov tests (P-values) | |||||
---|---|---|---|---|---|---|

Experimental samples versus homogeneous Poisson processes | Experimental samples versus simulated realizations of RSA | |||||

G | F | K | G | F | K | |

1 | 0.952 | 0.952 | 0.468 | 0.799 | 0.952 | 0.994 |

2 | 0.799 | 0.999 | 0.699 | 0.035 | 0.999 | 0.906 |

3 | 0.998 | 1.000 | 0.994 | 0.952 | 0.999 | 0.994 |

4 | 0.998 | 0.999 | 0.906 | 0.998 | 1.000 | 0.994 |

5 | 0.998 | 1.000 | 0.699 | 0.999 | 1.000 | 0.994 |

6 | 0.952 | 0.999 | 0.906 | 0.952 | 0.999 | 0.994 |

7 | 0.236 | 0.952 | 0.994 | 0.134 | 0.998 | 0.994 |

8 | 0.998 | 0.999 | 0.906 | 0.998 | 1.000 | 0.999 |

9 | 0.071 | 0.134 | 0.994 | 0.998 | 0.998 | 1.000 |

10 | 0.007 | 0.236 | 0.906 | 0.998 | 0.952 | 1.000 |

Sample no. | Two-sample Kolmogorov–Smirnov tests (P-values) | |||||
---|---|---|---|---|---|---|

Experimental samples versus homogeneous Poisson processes | Experimental samples versus simulated realizations of RSA | |||||

G | F | K | G | F | K | |

1 | 0.952 | 0.952 | 0.468 | 0.799 | 0.952 | 0.994 |

2 | 0.799 | 0.999 | 0.699 | 0.035 | 0.999 | 0.906 |

3 | 0.998 | 1.000 | 0.994 | 0.952 | 0.999 | 0.994 |

4 | 0.998 | 0.999 | 0.906 | 0.998 | 1.000 | 0.994 |

5 | 0.998 | 1.000 | 0.699 | 0.999 | 1.000 | 0.994 |

6 | 0.952 | 0.999 | 0.906 | 0.952 | 0.999 | 0.994 |

7 | 0.236 | 0.952 | 0.994 | 0.134 | 0.998 | 0.994 |

8 | 0.998 | 0.999 | 0.906 | 0.998 | 1.000 | 0.999 |

9 | 0.071 | 0.134 | 0.994 | 0.998 | 0.998 | 1.000 |

10 | 0.007 | 0.236 | 0.906 | 0.998 | 0.952 | 1.000 |

Results were similar when all synapses (asymmetric and symmetric) were studied as a single group and when only asymmetric synapses were analyzed (Supplementary Table 2). This result was expected since asymmetric synapses represented 91.74% of the total number of synapses found in layer III. However, due to the small number of symmetric synapses present in our samples (the average number of symmetric synaptic junctions per sample was 14, see Supplementary Table 1), it was not possible to calculate the *G*, *F*, and *K* functions for them.

There is a simple way, however, to test whether asymmetric and symmetric synapses are intermingled at random. If they were intermingled at random (null hypothesis), the nearest neighbor of any given synapse should be asymmetric or symmetric with a probability that would only depend on the respective proportions of both types of synapses in the sample. On the other hand, if both types of synapses were spatially segregated, their respective nearest neighbors would tend to be of the same type. Note that in this case we did not measure the distances between centroids but counted how many synapses had a nearest neighbor of the same or different type. These counts showed that 704 asymmetric synapses had an asymmetric nearest neighbor and 38 had a symmetric nearest neighbor. In the case of symmetric synapses, 43 had an asymmetric nearest neighbor, while 20 had a symmetric one. These data were used to create a 2 × 2 contingency table showing both types of synapses against the type of their nearest neighbor (Table 3). The 2-tailed Fisher's exact test applied to the contingency table rejected the hypothesis that both types of synapses are intermingled at random (*P* < 0.0001). In fact, the number of symmetric synapses that have a symmetric nearest neighbor was more than 4 times higher than expected under the null hypothesis (Table 3), possibly indicating that symmetric synapses tend to group together.

Type of synapse | Type of nearest neighbor | ||
---|---|---|---|

Asymmetric | Symmetric | Total | |

Asymmetric | 704688.54 | 3853.46 | 742742.00 |

Symmetric | 4358.46 | 204.54 | 6363.00 |

Total | 747747.00 | 5858.00 | 805805.00 |

Type of synapse | Type of nearest neighbor | ||
---|---|---|---|

Asymmetric | Symmetric | Total | |

Asymmetric | 704688.54 | 3853.46 | 742742.00 |

Symmetric | 4358.46 | 204.54 | 6363.00 |

Total | 747747.00 | 5858.00 | 805805.00 |

The observed counts of synapses (**bold**) were obtained from the 10 samples of the neuropil of layer III. The expected counts (in italics) are calculated from the marginal totals. Synapses whose nearest neighbors were not known—those that were located closer to the boundaries than to any other synapse in the sample—were excluded from the analysis. The Fisher's exact test applied to the observed counts rejected the hypothesis that asymmetric and symmetric synapses are intermingled at random (*P* < 0.0001). The most salient mismatch between the observed and expected values was the number of symmetric synapses that had a symmetric nearest neighbor, which was 4.4 times higher than expected.

## Discussion

Using FIB/SEM technology, we have been able to accurately quantify the density and spatial distribution of synapses in stacks of images from the neuropil of layer III of the young rat somatosensory cortex. In addition, virtually all synapses present in the 3-dimensional samples could be unambiguously identified since they can be visualized in consecutive serial sections and, if necessary, they can be digitally resectioned at different planes to ascertain their identity as asymmetric or symmetric (Merchán-Pérez, Rodriguez, Alonso-Nanclares, et al. 2009). In this way, the technology used in the present study eliminates the relatively large pool of synapses that were labeled as “uncharacterized” in previous reports since they could not be classified as asymmetric or symmetric from single 2-dimensional images (see for example DeFelipe et al. 1999). Unfortunately, the mean density of synapses per unit volume that we have obtained in this study (mean of 0.94 synaptic junctions/µm^{3}) and the proportion of asymmetric and symmetric synapses could not be compared with previous studies since no detailed electron microscope data are available in this particular layer of the rat somatosensory cortex at the age analyzed in the present study.

Regarding the spatial position of synapses, our results indicate that the distribution of synaptic junctions in the neuropil is nearly random, only constrained by the fact that synapses cannot overlap in space and so their geometric centers or centroids cannot be too close to their neighbors. This conclusion is further supported by the finding that the spatial distribution of the centroids of synaptic junctions can be closely modeled by a random process where overlapping is prevented by assigning a volume to each centroid. In mathematical terms, this model is known as “random sequential adsorption” (Evans 1993; Stoyan and Schlather 2000; van Lieshout 2006). None of the functions *F*, *G*, and *K* alone suffice for the characterization of a point pattern, so we have used them in combination (O'Sullivan and Unwin 2002; Gaetan and Guyon 2009). Visual inspection of all these functions showed only small deviations between the experimental data and the RSA model, which was further corroborated by a battery of goodness-of-fit tests. In our particular case, the volumes that were randomly assigned to each centroid were not arbitrary, but were drawn from the distribution of Feret's diameters of the 3-dimensionally reconstructed synaptic junctions. Although the Feret's diameter is only a rough estimate of the sizes of synaptic junctions, the RSA model not only prevents the overlapping of synaptic junctions, but also very closely mimics the actual distribution of synapses in the neuropil. Moreover, since no other constraints were imposed, the positions of synapses are otherwise nearly independent of the position of their neighbors. This and the fact that the Feret's spheres occupy only 5.85% of the total volume of the neuropil help to explain why the distributions of centroids in the samples and the simulated RSA processes are still very similar to a CSR process.

It must be stressed that our claim that the spatial distribution of synapses follows an RSA process does not necessary imply that during development synapses are actually formed and positioned in space following the same rules as the RSA. However, it is interesting to consider the idea that a synapse could be formed anywhere in space where an axon terminal and a dendritic element may touch, provided this particular spot is not already occupied by a preexisting synapse. Indeed, a similar model based on random axonal outgrowth and competition for available space has already been proposed (Kaiser et al. 2009). Nevertheless, the same end result could be achieved in different ways, including the synaptic turnover that is known to occur during cortical development (Rakic et al. 1994), provided that synapses are randomly added and/or withdrawn from the population.

Since there seems to be no limitation to the position of any synapse except the space already occupied by other synapses, at least a subpopulation of them would be located very close to one another. That is the case in our samples, since *G* functions show that approximately 20–40% of synaptic junctions have a nearest neighbor with an intercentroid distance of <500 nm. Considering the actual size of synapses (Table 1), this would mean that some of them would be located side by side, opening the possibility that the neurotransmitter released by one synapse could reach its neighbor by diffusion and influence its behavior (Fuxe et al. 2007; Rusakov et al. 2011). However, this seems to be true mostly for glutamatergic synapses (asymmetric synapses, the most abundant type), because the number of GABAergic synapses (symmetric synapses) is so low that very few of them were found near another asymmetric synapse. This is an interesting difference between the neuropil and perisomatic innervation regarding the organization of synaptic circuits. Indeed, there are numerous glutamatergic and perisomatic GABAergic axon terminals in close proximity around the soma and proximal dendrites of pyramidal cells, providing a putative basis for nonsynaptic interactions between these 2 types of synapses (Merchán-Pérez, Rodriguez, Ribak, et al. 2009).

An issue that needs to be addressed is whether the densities of synaptic junctions and their nearly random distribution in space have any predictive value to calculate the possible number of synaptic inputs that a given neuron may receive in layer III. This is related to the so called Peter's rule (Braitenberg and Schüz 1998), which represents a generalization to all types of cortical synapses based on the conclusion of Peters and Feldman (1976) regarding the synaptic connections of geniculo-cortical afferents and cortical neurons in the rat visual cortex: “Although thalamocortical afferents terminate principally in layer IV, their distribution with respect to postsynaptic targets may be essentially random, in the sense that no specific types of neurons receive the afferents. Instead, all elements in layer IV that are capable of forming asymmetric synapses may be involved” (Peters and Feldman 1976).

Our results seem to support this assertion since the positions of synapses do not show any spatial preference, so the rule would apply to the neuropil of layer III as a whole. Our findings are also in line with the distribution of swellings along pyramidal cell axon collaterals stained by the Golgi method in the mouse neocortex. In that case, the intervals between the axon swellings (characterized as putative synapses) were found to follow a Poisson distribution, with a dead space caused by the spatial extent of the swellings (Hellwig et al. 1994). However, it must be pointed out that the densities of synapses do vary locally (Table 1), and that they are randomly distributed, so local deviations from the rule will arise simply by chance. In any case, it is clear that the formulation of Peter's rule should take into account that different neuron types have different bouton distributions, project preferentially to different types of neurons and/or to different parts of these neurons (dendritic tree, soma, axon initial segment), and that both the axon terminals and dendritic trees may or may not originate in the same cortical layer (Anderson et al. 2002; Binzegger et al. 2004; Mishchenko et al. 2010; Bock et al. 2011). Thus, our data represent a heterogeneous population of synapses, and our conclusion that they are randomly distributed in space should be maintained only for this population as a whole, with the only possible exception of inhibitory synapses (see below).

Spatial randomness does not necessarily mean unspecific connections, since this does not preclude the existence of subjacent specificities based on molecular affinities, electrical activity, the developmental history of synapses, or other mechanisms not detected by our methodology. Indeed, our technique is capable of distinguishing between excitatory and inhibitory synapses, but there is no indication as to their particular neuronal type of origin. Different groups of the synapses in any given sample could come from cells located in the immediate vicinity, from other layers, from cortical areas outside the somatosensory cortex, from the contralateral cortex, or from subcortical regions. If all these different types of synapses were to form segregated groups (e.g. DeFelipe and Jones 1991), it would not be detected by our present methodology. However, we did find evidence that symmetric synapses tend to be spatially segregated from asymmetric synapses. Although we could not apply the same tests to asymmetric and symmetric synapses due to the scarcity of the latter, we found that symmetric synapses have a nearest neighbor of the same type over 4 times more than would be expected if they were intermingled at random with asymmetric synapses. A possible explanation is that these synapses may originate from one or few single axons, giving “en passant” synapses or ramifying locally, as typically occurs with GABAergic interneurons (Ascoli et al. 2008).

Another possibility is that spatial specificity in the neocortex is scale-dependent. It is well known that, at the macroscopic and mesoscopic scales, the mammalian nervous system is a highly ordered and stereotyped structure where connections are established in a highly specific and ordered way. Even at the microscopic level, it is also clear that different areas and layers of the cortex receive specific inputs. At the ultrastructural level, however, our results seem to indicate that synapses are distributed in a nearly random pattern. This could mean that, as the axon terminals reach their destination, the spatial resolution achieved by them is fine enough to find a specific cortical layer, but not sufficiently fine to make a synapse on a smaller target such as a specific dendritic branch or dendritic spine within that layer. For example, axon terminals from thalamic nuclei reach specific areas and layers of the cerebral cortex but, once there, they would form synapses randomly among their possible targets. In other words, at this ultrastructural scale, the specificity of connections would not only rely on spatial cues but also on other mechanisms such as molecular or activity-dependent cues in order to assure synaptic specificity.

## Supplementary Material

Supplementary material can be found at: http://www.cercor.oxfordjournals.org/.

## Funding

The study was supported by Centre for Networked Biomedical Research on Neurodegenerative Diseases (CIBERNED, CB06/05/0066) and Spanish Ministry of Economy and Competitiveness (grants SAF2009-09394, SAF 2010-18218, TIN2010-21289-C02-02, TIN2010-20900-C04-04, and the Cajal Blue Brain Project, Spanish partner of the Blue Brain Project initiative from EPFL). Funding to pay the Open Access publication charges for this article was provided by the Cajal Blue Brain Project.

## Notes

We thank L. Valdés and I. Fernaud for technical assistance. *Conflict of Interest*. None declared.