Abstract
In all species examined, with the exception of rodents, the axons of neocortical neurons form boutons in multiple separate clusters. Most descriptions of clusters are anecdotal, so here we developed an objective method for identifying clusters. We applied a meanshift clusteralgorithm to threedimensional reconstructions of 39 individual neurons and three thalamic afferents from the cat primary visual cortex. Both spiny (20 of 26) and smooth (7 of 13) neurons formed at least two distinct ellipsoidal clusters (range, 2–7). For all cell types, cluster formation is heterogenous, but is regulated so that cluster size and the number of boutons allocated to a cluster equalize with increasing number of clusters formed by a neuron. The bouton density within a cluster is inversely related to the spatial scale of the axon, resulting in a four times greater density for smooth neurons than for spiny neurons. Thus, the inhibitory action of the smooth neurons is much more concentrated and focal than the excitatory action of spiny neurons. The cluster with the highest number of boutons (primary cluster) was typically located around or above the soma of the parent neuron. The distance to the next cluster was proportional to the diameter of the primary cluster, suggesting that there is an optimal distance and spatial focus of the lateral influence of a neuron. The lateral spread of clustered axons may thus support a spokelike network architecture that routes signals to localized sites, thereby reducing signal correlation and redundancy.
Introduction
Studies of the local circuits in the neocortex have concentrated primarily on the laminar distribution of the different types of cortical neurons and on their interlaminar connections (Gilbert and Wiesel, 1979, 1983; Martin and Whitteridge, 1984). However, axonal and dendritic trees are structures of considerable spatial complexity (Gilbert and Wiesel, 1981; Binzegger et al., 2004). In the cat, tree shrew, and monkey cortex, but not in rodents (Hooser et al., 2006), the axons of many spiny neurons form long laterally directed collaterals that branch to form clusters of terminal boutons (Rockland and Lund, 1982, 1983; Gilbert and Wiesel, 1983; Martin and Whitteridge, 1984). These patchy projections were first seen for the ocular dominance system in area 17 of the primate, where the thalamic afferents not only form the retinotopic map on the cortex but also segregate anatomically into left and right eye dominated clusters or “columns” (Hubel and Wiesel, 1969). Such clustering also is particularly evident in the superficial cortical layers, where bulk injections of tracers revealed patchy retrograde and anterograde labeling (Rockland and Lund, 1982, 1983). Although the importance of clusters for the emergent physiology was apparent for the ocular dominance system, it has been more difficult to discover the relationship and the role of axonal clusters to the many functional maps revealed by optical imaging of the superficial cortical layers. It is generally supposed that excitatory neurons in the superficial layers connect to other neurons that have similar receptive field properties, although the evidence for this is equivocal (Malach et al., 1993; Kisvárday et al., 1996; Bosking et al., 1997). Another view is that the clusters are a means of optimizing the “wire length” of the cortical connections (Koulakov and Chklovskii, 2001).
Although these studies suggest a relationship between the functional maps and the axonal arbors of spiny neurons, it is clear that our understanding of the significance of the clusters for the physiology remains poorly understood even for the best studied cortical area: the primary visual cortex of the cat. Until now, clusters themselves have been defined only by qualitative criteria that depend on the eye of the observer. As a consequence, we know almost nothing quantitatively about the patch organization of individual neurons and how the patch organization might vary between different neurons and between different cell types. Still, less do we understand the constraints or rules by which an individual neuron apportions out its total complement of boutons to its various patches.
Here, we make a systematic study of the threedimensional (3D) clusters formed by examples of the most common cell types in cat area 17, including smooth (GABAergic) neurons and thalamic afferents. Our aim was to understand quantitatively the clustered structure of cortical axons. To do this, we used the meanshift algorithm for objectively defining clusters and derived a quantitative description of the clusters formed by individual axons. We found generic rules by which individual neurons distributed their boutons in the 3D volume of the neuropil.
Materials and Methods
Preparation and maintenance of animals
The neurons examined in this study were obtained from anesthetized adult cats that had been prepared for in vivo intracellular recording (Martin and Whitteridge, 1984) (for details, see Douglas et al., 1991). All experiments were performed under the authorization of animal research licenses granted to Kevan Martin by the Home Office of the United Kingdom and the Cantonal Veterinary Authority of Zürich (Switzerland). We first recorded from each cell extracellularly and mapped the receptive field orientation preference, size, type, binocularity, and direction preference by hand. The mapping was repeated intracellularly and horseradish peroxidase (HRP) was then ionophoresed into the cell. Thalamic afferents were classified as X or Ylike using a battery of tests (Friedlander and Stanford, 1984; Martin and Whitteridge, 1984). After appropriate survival times, the brains were fixed. The block of tissue containing the intracellularly filled neurons was serially sectioned in the coronal plane at a thickness of 80 μm and processed to reveal the HRP and osmicated and embedded in resin to reduce differential shrinkage (Anderson et al., 1994). This processing allowed the material to be examined at both the light and the electron microscopic level.
Horizontal tissue shrinkage was estimated to be 11%. The shrinkage was measured by placing tungsten microelectrodes at known distances in the cortex in vivo, fixing the brain, and then measuring the distances between the electrodes tracks in 80 μm sections. This gave the average shrinkage resulting from fixation. The shrinkage caused by tissue preparation was assessed by measuring the thickness of single sections immediately after vibratome sectioning and through osmication to the final embedding in Epon resin. This shrinkage was found to be negligible.
Cell reconstructions
Neurons were reconstructed in three dimensions with the aid of a light microscope (Leitz Dialux 22) and drawing tube attached to an inhouse 3D reconstruction system (Botha et al., 1987). Neurons were reconstructed at 400× magnification. Checks with higher magnifications confirm that a magnification 400× allows the fine details of individual structures (boutons, axons, and dendrites) to be identified with very high confidence. The reconstructions were characterized by a list of data points and stored for additional usage. Each data point consists of a code describing the digitized structure (axon, bouton, or dendrite) and its three spatial coordinates and thickness (where relevant). Axonal and dendritic collaterals are represented by open, piecewise linear curves. The axonal arborizations are complex and often extend through many histological sections. The measurement error of the 3D position of digitized structures was estimated by measuring seven boutons (horizontal displacement <100 μm) within a single section 10 times. The SD was smaller than 0.7 μm in all three dimensions.
Laminar distribution of boutons and dendritic trees
For each section, the affiliation of the boutons, axonal, and dendritic segments in the section to one of the five cortical layers 1, 2/3 (i.e., layers 2 and 3 where combined), 4, 5, and 6 was determined using the lamina border criteria of Henry et al. (1979). The relative position of the laminar borders and the boutons, axonal, and dendritic trees of all sections was summarized in a single coronal projection (see Fig. 1). In general, the lamina borders of this summary diagram are distorted, and the diagram does not show the correct layer thickness. In particular, for those cells in which the axon extended over large (>1 mm) anteroposterior distances, it was not possible to represent in the single coronal projection both, the correct relationship between the reconstructed neuron and laminas boundaries, and the correct location of the boundaries.
Gaussian ellipsoids
For a trivariate Gaussian distribution, which was fitted to a set of points in the 3D space, we define the “αellipsoid” where C is the sample covariance matrix and μ is the sample mean of the fitted Gaussian distribution. The αellipsoid is centered at μ, and the axes have directions given by the three eigenvectors v_{k} of C (k = 1,2,3) and the length of each halfaxis is where λ_{k} is the eigenvalue of v_{k}. The volume of the αellipsoid is , and the elongation of the ellipsoid is expressed as the ratio of the smallest and largest eigenvalues. For points approximately distributed according to a trivariate Gaussian distribution, ∼75% of the points will be contained in the 2ellipsoid.
Bouton cluster identification using the meanshift algorithm
The threedimensional arborization pattern of an axon is typically heterogenous, composed of spatially separated regions with intense axonal arborizations and bouton formation (see Fig. 1). These “patches” have a high bouton density relative to the surrounding zone, and the density distribution therefore resembles a mountainlike landscape, where the mountain peaks indicate the regions of high bouton density, and the different regions are delineated by low density valleys. In this paper, we study these regions of high bouton density, referring to them as clusters.
A smooth density landscape was determined by convolving the bouton locations with a spherical Gaussian kernel of width (or SD) h (Scott, 1992). Each local maximum of the density landscape is a peak of a mountain and the set of boutons forming the mountain defines the cluster. To determine which boutons belong to the mountain, boutons are moved along the local gradient until a local maximum is reached and they can easily be identified. The well established “meanshift” algorithm is an iterative procedure that performs these steps without the need to calculate explicitly the density landscape and the gradient (Fukunaga and Hostetler, 1975; Cheng, 1995). We implemented an accelerated version of this algorithm (CarreiraPerpinan, 2006) with a resolution of 5 μm, stopping threshold 10^{−8} for the entropy difference and 10^{−3} for the shift in the bouton locations. The meanshift algorithm forms the heart of the clustering procedure, which involves three major steps as described below.
Preprocessing: elimination of “linear” structures.
Not all the boutons of an axon are contained in clusters, which is not surprising, because the location of individual boutons is primarily constrained by the linear branching pattern of the axon (Binzegger et al., 2005) and almost every axon contains isolated axonal branches. These were identified in a preprocessing step. We also excluded clusters containing only a small number of axonal branches, or ringlike arrangements where the center of the cluster was bouton free. Subsets of boutons arranged in long onedimensional strings (typically connecting the cell body with a patch) (see Fig. 1) were reduced before the meanshift algorithm was applied. To do so, the bouton cloud was partitioned into small local regions where each region was classified as being part of a bigger linear set of points L if its boutons were formed by only one or two axonal branches, or the fitted 2ellipsoid was highly elongated (<0.1), or the region of the fitted 2ellipsoid was very small (<5 μm^{3}). Partitioning was done using the meanshift algorithm with a spherical Gaussian kernel with width ρ. The set L will then depend on ρ where a large ρ will produce a small set L, and a small ρ will produce a large set L containing most of the boutons in the cloud. We chose ρ so that the fractal dimension of L was close to 1 (between 0.95 and 1.05). An example set L is shown for one neuron in Figure S1 (available at www.jneurosci.org as supplemental material). The proportion of boutons in L ranged from 1 to 44% (14 ± 10%) for the different neurons.
We found that inhibitory neurons have a low L value (12 of 13 neurons; L between 1 and 10%), and excitatory neurons tended to have a high L value (17 of 26 neurons; L between 15 and 25%). Excitatory neurons tend to form long isolated axonal paths, often connecting clusters to the soma (see Fig. 1A2). Inhibitory neurons (in particular, basket cells) can also form long paths, but they give off many small twiglets (see Fig. 1C2). We regard these “bushy” paths as not linear, and they are not detected by the preprocessing step. Some excitatory neurons simply did not form many long axonal arbors (see Fig. 1D1), and those neurons also had a low L value.
Choosing the appropriate width h* for the convolution kernel.
Boutons may also display clustering at more than one spatial scale. For example, measured distributions show the characteristics of fractal sets (Binzegger et al., 2005). A thalamic afferent illustrated by Gilbert and Wiesel (1983) had two major patches, corresponding to ocular dominance columns, and additionally, within the major patch there was a subsystem of microclusters spaced ∼100 μm. A similarly spaced subsystem was observed for a small basket cell of layer 4 (Kisvárday et al., 1985). Of course, the cluster algorithm can detect these fine grain undulations, but we chose a parameter (h) that detected only the coarser grain clusters associated with the columnar systems of visual cortex (200–400 μm) (Gilbert and Wiesel, 1983; Martin and Whitteridge, 1984; Freund et al., 1985; Humphrey et al., 1985; Kisvárday et al., 1985, 1986; Gabbott et al., 1987; Martin, 1988; Kisvárday and Eysel, 1992). This choice of parameter h was done in the following way.
After elimination of L, we partitioned the remaining boutons into clusters using the meanshift algorithm width kernel width h. The choice of h controls the smoothness of the density landscape. It can be thought as the equivalent of the bin width in histograms. A large h will result in a very smooth landscape where only the gross features are represented, and the smaller h is chosen, the more local maxima appear. To select the appropriate width h* for a bouton cloud, we used the intuitive idea that for reasonably separated high density regions, there should exist an interval of h values for which the density distributions (and the resulting clusters) between consecutive h values do not differ much from each other. To find this stable regime, we applied the meanshift algorithm to a large range of h values (30–250 μm for excitatory neurons and thalamic afferents; 10–150 μm of inhibitory neurons; step size, 5 μm) and so determined for each h value the resulting partition of the boutons into clusters. We then used an index to quantify the similarity between consecutive partitions. The similarity index is a value between 0 and 1 with 1 indicating equal partitions. It is defined as 1 − m/(u − 1), where m is the minimal number of points that had to be deleted in order for the two partitions to be equal, and u is the number of boutons (Giurcaneanu and Tabus, 2004). We stopped evaluating the similarity index if for an h value the resulting cluster contained all boutons in the cloud. It is trivially true that for all larger values of h, the similarity index will be constant 1. The width h* was then selected from the interval of h values where the similarity index did not change considerably (Figs. S2, S3, available at www.jneurosci.org as supplemental material).
We chose the first such interval of length 20 μm for which the similarity index was >0.99. For only 5 of the 39 cells (4 of a total of 11 basket cells, and the double bouquet cell), no such stable regime existed. Smaller intervals still existed (10–15 μm), which indicates that some heterogeneity still existed in the bouton cloud. Inspection of the 3D bouton clouds using different viewing angles confirmed this interpretation. We decided to include those five neurons in the sample using the h* from the smaller interval. Their inclusion did not change the results in any significant way.
Of course, the length of the interval that was considered as a stable regime had to be selected conservatively, because bouton clouds can display patches at more than one spatial scale. The length of 20 μm for a stable regime (or 10–15 μm for the five exceptions) was long enough to prevent the selection in clusters systems much smaller in scale than is usually observed in columnar systems of the primary visual cortex (200–400 μm).
After every application of the meanshift algorithm, clusters with 2ellipsoids of small volume (<5 μm^{3}) and large elongation (<0.1) were ignored. In addition, if neighboring clusters were not well separated by a reasonably deep valley, they were merged. In particular, two clusters were merged if the 3ellipsoids intersected and f_{min}/f_{mean} > t(t = 0.85) where f_{min} is the minimal bouton density and f_{mean} is the mean bouton density along the straight line connecting the two cluster centers. The parameter t was selected small enough as not to allow for the identification of spurious clusters, which are separated by shallow valleys and typically occur in homogenous like point clouds (Fig. S2B, available at www.jneurosci.org as supplemental material) but large enough so that moderately separated peaks are still separated (Fig. S2A, available at www.jneurosci.org as supplemental material).
Postprocessing: excluding clusters from analysis.
Clusters obtained from the meanshift algorithm (with kernel width h*) with predominately inhomogeneous bouton arrangements were excluded from analysis. In particular, clusters containing only a couple of branches (B = 15) were excluded, as well as clusters with no boutons within the αellipsoid (α = 0.55). The latter criteria excludes “ringlike” arrangements, where no boutons are present close to the center of the cluster. For the 39 reconstructed neurons, a total of 193 clusters were identified, of which 44% were excluded. Most of these clusters did not contain many boutons, and only 9% of the total complement of all boutons (not counting the boutons that were excluded in the preprocessing step) were excluded. For a change of B by ± 5, the fraction of excluded clusters ranged between 38 and 48% (excluded boutons 7–10%). For a change of α by ± 0.2, the excluded clusters ranged between 35 and 62% (excluded boutons 5–21%).
Basic cluster measurements were based on the 2ellipsoids. The “cluster center” is the center of the 2ellipsoid. The “cluster diameter” is defined as the geometric mean of the three diameters of the 2ellipsoid [i.e., d = (d_{1}d_{2}d_{3})^{1/3} with d_{k} = 4λ^{1/2}. It is also the diameter of the sphere with the volume of the 2ellipsoid. The “cluster bouton density” is measured by the ratio of the number of boutons in the 2ellipsoid and the volume of the 2ellipsoid. The “relative diameter” of a cluster is the ratio of the cluster diameter and the summed diameters of all the clusters formed by the axon. The “cluster weight” is the ratio of the boutons in the cluster and the summed bouton number in all clusters. In the cat, the cortical layers of area 17 as seen in transverse sections are strongly curved on the apex of the gyrus, where most of the neurons were recorded (see Fig. 1). To investigate the lateral arrangement of the clusters, we corrected for this curvature and applied a transformation to each 2ellipsoid, which maintained its relative position and orientation in the flattened layer (Fig. S4, available at www.jneurosci.org as supplemental material). For simplicity, any curvature along the dorsoventral axis was ignored.
An alternative approach is to apply the flattening transformation to the boutons and perform the clustering on the flattened distribution. However, this would have introduced artificial distortions, which we wanted to avoid (the distance between some boutons is stretched, whereas for others it is reduced). No correlation was found between cluster diameter and layer width, so we also refrained from rescaling the bouton distributions to account for variations in cortical thickness.
A generative model for cluster formation
We start with a neuron that forms one cluster containing B boutons (see Fig. 13A, row n = 1, black disc). If the neuron forms no additional clusters, we are done. In the case where another cluster is formed, we assume that the new cluster is formed from a fraction γ < 1 of the boutons from the existing cluster. After splitting, the original cluster contains (1 − γB) boutons, and the new cluster γB boutons. The new cluster is indicated in Figure 13A by an open disc (row n = 2). One might imagine that such a splitting procedure could involve the formation of one or more axonal paths extending from the original cluster to the new cluster. In Figure 13, this is indicated by an arrow between the discs. More generally, for a neuron forming n ≥ 3 clusters, the cluster weights are given by w_{1}=(1 − γ)γ^{0}, w_{2} = (1 − γ)γ^{1},…, w_{n}_{1} = (1 − γ)γ^{n−}^{2}, w_{n} = γ^{n−}^{1} in the “linear chain” arrangement, where newly formed patches are created in succession by splitting from the cluster that was created previously. The other extreme case is the “spoke” arrangement, where clusters are always created by splitting the root cluster. This results in the weights w_{1} = (1 − γ)^{n}^{−1}, w_{2} = γ (1 − γ)^{0}, w_{3} = γ(1 − γ)^{1},…, w_{n} = γ(1 − γ)^{n−}^{2}.
Results
We recorded from and filled 39 cells or axons with principal innervation that was area 17 of the cat. After processing, which included osmication and embedding in Epon to minimize the shrinkage associated with conventional histological methods, we digitized the entire 3D distribution of the boutons formed by the axons. Three axons were afferents of the dorsal lateral geniculate nucleus. The remaining axons originated from examples of cortical pyramidal cells (18 of 39), spiny stellate cells (5 of 39), basket cells (12 of 39), and one double bouquet cell. These 3D data were the raw material to which we applied a cluster algorithm, which allowed us to define salient characteristics of the bouton patches formed by the different cell types. Examples of patchy bouton clouds and the identified clusters are given in Figure 1 (for a summary diagram showing the clusters of all neurons, see Figure 2). Our neurons exhibited between one and seven clusters, confirming the few previous reports on patch number in the literature that report between four and eight and patches per neuron for superficial pyramidal cells (Kisvárday and Eysel, 1992), up to seven clusters for spiny stellate cells (Martin and Whitteridge, 1984), and between one and four patches per thalamic afferent axon (Freund et al., 1985).
Clustering is a fundamental property of spatial bouton organization
Figure 1A shows the spatial distribution for the boutons of a pyramidal cell of layers 2 and 3 (p2/3) in coronal view and top view. This neuron formed 7069 boutons, most of them in layers 2 and 3, which is the primary layer of innervation for all pyramidal neurons that we analyzed in the superficial cortical layers. The spatial arrangement of the boutons appears highly inhomogeneous, featuring several regions of high bouton density in layers 2, 3 and in layer 5. We used the meanshift algorithm to identify those high density regions, and we refer to them here as clusters (see Materials and Methods). For this layer 2/3 pyramidal cell, the algorithm found six clusters. The association of the boutons to a cluster is indicated by different colors, where the color codes for the relative number of boutons each cluster contains (i.e., black indicates the cluster with the highest number of boutons, followed by red, green, and blue, etc.). We also refer to the cluster with the highest number of boutons as the rank 1 cluster, the cluster with the second highest bouton number the rank 2 cluster, and so on.
Although the meanshift algorithm identified high density regions of any shape, we focused on the clusters that had an ellipsoidal appearance (Fig. 1). These clusters were typically elongated rather than spherical, with the ratio of shortest to longest diameter being approximately constant (0.37) (Fig. S5, available at www.jneurosci.org as supplemental material). For simplicity, we will always describe the size of the cluster by the geometric mean of the three diameters and refer to it as the cluster diameter (see Materials and Methods).
It is often difficult to grasp the 3D organization of the boutons, and the 2D projection may be misleading. For example, the group of boutons at the top right of Figure 1A2 appears as a patch formed by the convergence of three long axonal arbours. A close inspection using other viewing angles reveals that the three arbours target the same vertical column but at quite different vertical positions. Similarly, in Figure 1A1, the cluster algorithm identified a small patch (purple), whereas seemingly ignoring other, similar patch like structures in layer 5. Again, inspection of the bouton cloud from different viewing angles confirmed that the cluster labeled purple is a compact structure of high density, whereas the grouping of boutons in layer 5 is only present in one particular view (compare with Fig. 1A2).
Although multiple cluster formation is primarily associated with superficial layer pyramidal cells, our data suggest that the arrangement of boutons into several discrete clusters of high bouton density is a characteristic of most cortical axons. Almost all neurons (32 of 39) formed at least two clusters (3.2 ± 1.4; range, 2–7) (Fig. 3B,C). For the remaining seven neurons, the algorithm identified one cluster. In some of these cases, a significant number of the boutons were not assembled into a cluster but were diffusely distributed. The identified cluster was the only dense ellipsoidal set of boutons, whereas the diffuse part was eliminated by the preprocessing and postprocessing stages of the cluster algorithm (Fig. 1E). For other neurons, the cluster contained almost all boutons of the cloud, in which case the word cluster may be less appropriate. However, for all those axons, the bouton cloud was highly compact (Fig. 1D), and we think it is justified to call them cluster and analyze them together with the other clusters.
Three quarters of the boutons in the cloud typically belonged to a cluster (range, 37–98%) (Fig. 3A). Our sample of clusters (n = 109) had diameters ranging from 87–941 μm (373 ± 156 μm) (Fig. 3E), and each cluster contained between 70–8265 boutons (1143 ± 1392) (Fig. 3D) with cluster weights ranging from 0.02 to 1.0 (0.4 ± 0.3).
Excitatory neurons formed on average slightly more clusters than inhibitory neurons (3.0 ± 1.8, 1–7 vs 2.5 ± 1.0, 1–5). Although the formation of several distinct clusters could be observed for several inhibitory neurons, for five inhibitory neurons (two basket cells from layer 2/3, two basket cells from layer 4, and the double bouquet cell) the cluster algorithm could not segment the landscape well, indicating that these inhibitory neurons formed less distinct clusters. Superficial pyramidal cells were among the neurons with the most clusters, but interestingly, the cluster number ranged from 2 to 6 (n = 6), indicating that even for neurons of the same morphological type, there is heterogeneity in their axonal distribution. All reconstructed superficial layer pyramidal cells were located in layer 3, and no clear relationship was found between cluster number and radial position of cell body. Similar heterogeneity was also found for the layer 6 pyramidal cells [p6(L4); n = 6], which formed between one and five clusters (Figs. 1B,D, 3B). Neurons forming more than four clusters (8 of 39 neurons) also included a spiny stellate cell [ss4 (L4)], thalamic afferents (X/Y), and a superficial basket cell (b2/3) (Fig. 1C). However, in our set of neurons, the more common finding was that neurons formed two or three clusters (23 of 39 neurons).
Cluster number determines cluster weight
For most neurons, there is a marked difference in the numbers of boutons contained in rank 1 versus rank 2 clusters (Fig. 4). The rank 1 cluster of an excitatory neuron forms on average 2036 ± 1592 (376–8265) boutons. The rank 2 cluster typically contains already four times less boutons (590 ± 281; 195–1252). For the inhibitory neurons, the difference between rank 1 and rank 2 clusters is even more dramatic, falling from 3336 ± 1273 (1096–5201) boutons per cluster to 386 ± 303 (70–1034) boutons per cluster (i.e., an almost 10fold reduction). More generally, the mean number of boutons in a given cluster rank decays with increasing rank like a power function (Fig. 5C,D).
The weight distribution ranged from neurons that had one dominant patch, which contained most of the boutons, to neurons in which the boutons were distributed more equitably among multiple patches. To understand how the clusters of a given rank compared between individual neurons, we plotted for each neuron the bouton number per cluster as a function of rank (Fig. 5A). No simple pattern emerged. For example, the rank 1 cluster of neurons with two clusters (Fig. 5A, black lines) and five clusters (stippled dark gray lines) overlapped considerably. Indeed, for both the excitatory and inhibitory neurons, the correlation between rank 1 clusters and cluster number was abysmal (r^{2} < 0.3).
Because the total number of boutons in a cloud varied greatly between neurons (Fig. 3A), we tested whether a pattern emerges when the cluster weight (the proportion of the total boutons contained in a given cluster) of a given rank, rather than the absolute bouton number, was compared between the neurons (Fig. 5B). The correlation of rank 1 weights with cluster number indeed turned out to be statistically significant (all neurons r^{2} = 0.69, p < 0.01; excitatory neurons r^{2} = 0.78, p < 0.01). For the excitatory neurons in particular, we observed that the decay of the weights slows down with increasing cluster number per neuron (Fig. 5B, comparing again the two groups with two and five clusters).
We quantified these observations by fitting a power function A_{n}k^{−Bn} to the cluster weight distributions of all neurons with n = 1,…,7 clusters (Fig. 6A). The fits were constrained to power functions for which = 1, because we are fitting cluster weights, of which sum over all ranks equals 1. This means that the coefficient A_{n} is simply the rank 1 cluster weight for neurons forming n clusters, whereas the exponent B_{n} determines how fast the function decays with increasing cluster rank (the bigger A_{n}, the faster the decay). As anticipated, the coefficient decreases (linearly) for increasing number of clusters per neuron (Fig. 6B). This means that the more clusters a neuron forms, the lower is the rank 1 cluster of a neuron. The exponent B_{n} decreases (also linearly) with cluster number per neuron (excluding the trivial case n = 1) (Fig. 6C) (i.e., the decay is slower for neurons with high cluster numbers), and boutons are allocated more equally between the different clusters.
We conclude that, independently of the morphological cell type of a neuron, the number of boutons in its clusters tends to be determined by the number of clusters (n) and the total number of boutons (B). If k is the cluster rank, the bouton number in the rank k cluster is given by coefficient A_{n} and the exponent B_{n}, which depends linearly on n. Figure 6D makes a direct comparison of the observed bouton number per cluster with the predicted bouton number per cluster, indicating a significance correspondence (r^{2} = 0.95; p < 0.01). An equivalent view is that the number of boutons in the rank 1 cluster, together with the cluster number, specifies the number of boutons in the remaining clusters. This follows from u_{k}/u_{1} = w_{k}/w_{1,} where w_{k} is the weight and u_{k} the bouton number of cluster k.
Cluster number determines relative cluster diameter
Although by definition the bouton number per cluster decays with rank, the dependence of the cluster diameter with rank is a priori not clear. We found that the cluster diameter is typically larger for rank 1 clusters than for rank 2 clusters (Fig. 7). For excitatory neurons, the diameter dropped from 541 ± 154 μm (308–941 μm) to 402 ± 99 μm (225 ± 576 μm), and for inhibitory neurons from 377 ± 78 μm (258 ± 571 μm) to 218 ± 108 μm (87 ± 419 μm). In general, the rankwise average diameter decayed for both classes like a power function, with the rankwise average of the excitatory neurons always ∼1.5 times greater than that of the inhibitory neurons (Fig. 8C,D).
As for the observation made in the last section that the correlation between rank 1 clusters and cluster number was poor, a comparison of the diameter between clusters of different neurons also did not reveal a clear pattern (Fig. 8A). The correlation between first rank diameter and cluster number was again poor for both the excitatory and the inhibitory neurons (r^{2} < 0.15). However, an inspection of the cluster systems formed by the different neurons showed that they operate on quite different spatial scales (Fig. 2), which suggests that it may be more appropriate to compare the relative cluster diameter between the different neurons (Fig. 8B). When we did this, we found a significant correlation between relative diameter and cluster number (all neurons, r^{2} = 0.78, p < 0.01; excitatory neurons, r^{2} = 0.79, p < 0.01), and the decay in diameter with cluster rank changed linearly as a function of cluster number per neuron (Fig. 9). Thus, we find an analogous result for the cluster diameter as we did for the bouton number per cluster: The diameter of the rank 1 cluster, together with the cluster number, tends to determine the diameter of the remaining clusters formed by the neuron.
The covariation of cluster diameter and number of boutons per cluster might be such that the bouton density is conserved for the different patches. The histogram of cluster density (Fig. 3F) shows that this is not globally true. Clusters from excitatory neurons had approximately four times lower bouton densities than clusters from inhibitory neurons [2.24 ± 1.42 per (50 μm)^{3} vs 9.63 ± 6.69 per (50 μm)^{3}], and even within cell types, cluster density varied still considerably. This variation in density is evident even in single sections [Ahmed et al. (1994), their Fig. 2]. What we observed, however, is a tendency for individual neurons to have constant cluster density, but the constant density could still vary considerably between neurons (Fig. S6, available at www.jneurosci.org as supplemental material).
Vertical and horizontal organization of clusters
It is evident that most neurons contribute their boutons to more than one layer of cortex. However, it is not known whether individual clusters are confined to a single laminar or cross lamina boundaries. To analyze the laminar distribution of the clusters, we assigned each cluster to the cortical layer where the center of gravity was located. We found that typically, a cluster with origin in a given cortical layer also had most of its boutons in that layer (93 ± 12%; 43–100%).
Figure 10 shows for each neuron the layer location of the clusters together with their cluster weights. With a few exceptions (two layer 4 pyramidal cells, p4), neurons of the same cell type had rank 1 clusters always in the same cortical layer (indicated with a shaded area). We call this layer the primary layer of axonal innervation, because neurons within a cell type form most of the axon and boutons in this layer (Binzegger et al., 2004). Rank 2 clusters (large open dots) were more variable, but still 29 of 39 neurons formed this cluster in the primary layer of axonal innervation. Pyramidal cells in layer 2/3 primarily innervate layers 2/3 but also have a prominent axonal arborization in layer 5. It is therefore not surprising to find that three of the six pyramidal cells of layer 2/3 had rank 2 clusters in this layer (Fig. 10).
Within the primary layer of axonal innervation, we also observed clear trends in the horizontal organization of cluster. Figure 2 indicates that if more than two clusters were formed in a layer, the clusters are asymmetrically distributed relative to the soma location and often tend to be aligned along a specific axis or biased toward a halfplane. However, because of the relatively small sample size, we did not attempt to quantify this observation.
The extent of the area over which clusters were distributed was clearly different for smooth and spiny neurons. Within the layer of primary innervation, the lateral displacement of clusters can extend up to 2 mm from the origin (layer 2/3 pyramidal cell), but only spiny neurons were found to form clusters with displacements of >1 mm (Fig. 11). If a smooth neuron did form a distal cluster, it was of low weight, indicating that smooth neurons do not contact many distal targets. With the exception of rank 1 clusters, cluster rank was not convincingly correlated with cluster location. For most neurons (30 of 36, excluding thalamic afferents), the rank 1 cluster had a horizontal displacement from soma that was <250 μm, indicating that neurons allocate a significant proportion of their boutons locally, either around the soma, or in another layer above the soma (Figs. 2, 11). For those neurons, the rank 2 cluster was located more distally with a horizontal distance ranging widely from 0.5–2 mm. Thus, for more distal patches, there is no suggestion of a trend to decrease the number of boutons per cluster with increasing horizontal displacement.
Figure 2 indicates that neurons operate on different spatial scales. For example, when comparing the horizontal organization of basket cells and pyramidal cells of layer 2/3, it is easy to recognize that the area over which the pyramidal cells distribute their clusters is much bigger than that of the basket cells. Correlated with this increase in area is also a greater cluster diameter and a greater centertocenter distance of the pyramidal cell clusters. In contrast, inhibitory neurons tend to have adjoining clusters. We quantified this observation by measuring for each neuron the diameter of the rank 1 cluster and correlated this diameter with the horizontal displacement to the next closest cluster in the primary layer of axonal innervation (if such a cluster existed). The correlation was significant (p < 0.01) (Fig. 12).
Generative cluster model
We explored the possibility that the observed distribution of weights across neurons could be explained by changes in a small number of parameters of a single simple growth rule (see Materials and Methods) (Fig. 13A). The model supposes that a neuron has a fixed pool of boutons that it can form, and the way it distributes these boutons depends on its branching topology. This dependence arises out of the constraint that each daughter cluster receives a constant (for the axon) fraction of the boutons present in its parent cluster. Thus, an axon that grows successive clusters by chainlike extension of its single trunk will have a very different distribution of bouton number than an axon that grows all its clusters by spokelike branching (Fig. 13B,C). Although the model is simple, it is able to predict the observed weight distributions successively in the case of the spokelike arrangement (Fig. 13D).
Discussion
Because the remarkable single neuron labeling experiments of Gilbert and Wiesel (1979, 1983), many papers have reported patchy or latticelike patterns of axonal arborizations in the neocortex. These data have usually been presented in the form of light micrographs of single tangential sections. We have now investigated such bouton distributions quantitatively and have defined the axonal formation of clusters of boutons by means of an objective clustering algorithm. Consequently, we can now make a distinction between a “patch,” which denotes a cloud of boutons identified by eye, and a “cluster,” which denotes an algorithmically detected cloud. Our source data were detailed 3D reconstructions of the entire axon and bouton locations of single intracellularly labeled neurons and thalamic afferents in the cat primary visual cortex.
Previous descriptions have emphasized the patchy arborizations of the thalamic afferents and superficial layer pyramidal cells. We were surprised to find by our quantitative analysis that clustering of boutons is a much more universal feature of cortical neurons than previously appreciated. We found that 80% of all neurons, including both excitatory (20 of 26) and inhibitory (12 of 13) neurons, formed clusters. This striking result immediately raises questions concerning the performance of the clustering algorithm, on which our measurements of clustering depends. We will thus address this technical question before discussing our results further in detail.
Identifying patches using a cluster algorithm
Many algorithms are available for partitioning heterogeneous clouds of points into clusters (Buhmann, 1998; Xu and Wunsch, 2005). Because clustering is a form of induction, we cannot be assured of a unique solution and method of validation. At best, one can expect a method that detects stable clusters that are natural to the data. In the case of the bouton data, “natural” means that the algorithm is able to detect clusters that even an unbiased observer would agree with, and by “stable,” we mean that the algorithm detects the same set of clusters over a reasonable range of clustering parameters. We tested a few candidates against these properties, including expectationmaximization (EM) of Gaussian mixtures (McLachlan and Basford, 1988), before settling on the meanshift algorithm, because it gave consistently the most meaningful results over our entire data set. We also verified the performance of the meanshift clustering procedure using synthetic data drawn from mixture of trivariate Gaussian distributions (Fig. S2A, available at www.jneurosci.org as supplemental material). The algorithm successfully recovered all components, so long as they did not overlap too much. Importantly, for this study, the meanshift algorithm was able to resolve clusters with considerably different volume, elongation, and point densities and also with unequal distances between them.
Rules of cluster formation
By our analysis, 80% of neurons have more than one cluster. When these multiple clusters are ranked according to the number of boutons that they contain, the dominant, or primary, cluster is almost always (90%) located in the same radial column (but not necessarily in the same layer) as its parent soma. That is, the horizontal ellipsoid that is tangential to the border of the cluster will almost always enclose the source soma. This is not true of the secondary clusters.
We found that if clusters are ranked by the number of boutons they contained of the total number of boutons formed (the “cluster weight”), these weights decay as a power function of rank (Fig. 5). The exponent of decay correlated with the number of clusters (Figs. 6, 9). The greater the number of clusters, the smaller the exponent, and the more similar are the numbers of boutons between the clusters. We were surprised to find such an orderly process, and so looked for additional support from published pictures of axonal trees of spiny neurons. A brief survey suggests (without the benefit of direct measurement) that such a range of patch weight distributions are common (Gilbert and Wiesel, 1983; Martin and Whitteridge, 1984; Kisvárday et al., 1986; Gabbott et al., 1987; Martin, 1988; Kisvárday and Eysel, 1992).
A similar relationship holds for the cluster diameter (Fig. 8), indicating that the density of boutons is constant across the clusters of an individual neuron (Fig. S6, available at www.jneurosci.org as supplemental material). Moreover, bouton density is inversely related to the spatial scale of the cluster system of a neuron, so that superficial pyramidal cells with spatially extended cluster systems have lower bouton densities than inhibitory neurons, the clusters of which are smaller and closer to one another. As a class, the bouton density of smooth cells is approximately four times greater than that of excitatory neurons, a feature that is apparent even in light micrographs [Ahmed et al. (1994), their Fig. 2].
We consider this relationship between spatial scale and cluster density to be a fundamental principle of cortical organization. If all neurons generate clusters with the same constant density, independent of the spatial scale of the axonal arborization of a neuron, and about the same total bouton number, then compact arborization patterns are not possible. In this case, each small volume of cortex will contain boutons from the different cell types in direct proportion to their frequency. Instead, we find that different clusters have different densities, but that different cell types have approximately the same total number of boutons. As a result, neurons with compact clusters will have a much higher representation of boutons and perhaps a disproportionate local influence within a given small volume of cortex than neurons with axons that have a larger spatial extent. This anatomical organization principle supports fast, intensive signal innervation of cortical sites, which would not be possible by, say, simply increasing the cell density of classes of neurons with large axonal extent to strengthen their influence within a given small volume of cortex. The smooth cells are inhibitory, and although they form only ∼15% of the complement of cortical neurons, their dense axonal arbors permit them nevertheless to have a strong focal effect [Douglas and Martin (2004), their Fig. 4].
The predictable order in cluster properties suggests that there is a fundamental rule governing axonal growth and its associated bouton formation. We tested a possible model (see Materials and Methods) (Fig. 13), which forms bouton clusters with weights depending on the axonal branching topology. Our model predicts that the observed cluster organization can occur if successive clusters are formed in a spokelike manner emerging directly from the primary patch, but not if they are formed successively as chains. An inspection of the axonal topology confirms a predominance for a spokelike organization of the clusters.
Cluster organization as a means of reducing signal correlation
We discovered that there is a distinct order in the distance between clusters (Fig. 2). The distance between the bouton cluster nearest to the soma (which is usually also the primary cluster) and the next closest cluster is about one and onehalf diameters of the primary cluster (Fig. 12). This relationship is true for all our clustered neurons, suggesting that there is an optimal distance and spatial focus of the lateral effect of a cortical neuron. Interestingly, this relationship matches closely the relationship between the diameters of extracellularly labeled patches and their interpatch distance, as observed across various cortical areas and animal species (Fig. 12, black crosses). It seems then that this separation distance between clusters may be a general property of neocortical construction. It will be interesting to investigate how the clustering principles relate to the functional representation of the visual space in the primary visual cortex. For example, we observed that the horizontal cluster arrangement of the superficial and deep pyramidal cells tends to be correlated with the orientation preference of the neuron, but the low sample numbers did not permit a systematic analysis.
One might argue that this regular separation of clusters reflects the periodic organization of cortex, for example, the spacing between isoorientation domains, but this cannot be the whole answer, because our results show that the distance between clusters is not constant across all neurons (as expected for a single lattice), but rather that the intercluster distance scales with primary cluster size. What could the reason for this possibly be? We speculate that this distance may be required to reduce correlations between populations of neurons engaged in successive stages of a spike timedependent computation. The intuition is as follows: because of the concentration of effect in the primary cluster, it is likely that the spikes of the neurons within the primary cluster population will become correlated via their recurrence with respect to one another (Douglas et al., 1995). If the results (output spike trains) of these neurons (call them population A) should now be combined with statistically independent output from population B by converging onto population C, then that combination should occur at a “safe” distance from both A and B. That means that the output clusters of A and B should converge onto C, but C needs to be at a sufficient distance to avoid statistical cross talk between the source populations A and B.
Selective routing
Overall, a spokelike construction method, together with a regular separation distance between primary and subsidiary clusters, indicates that cortical neurons in higher mammals have a predilection for constructing regular latticelike patterns of connections with nearby neighborhoods. Such a coherent pattern of clustered connections, which we have called a “daisy” (Douglas and Martin, 2004) would have to arise from the collective connectivity properties of a population of local neurons. In the context of daisylike periodic structures formed by local neurons, one could consider two fundamental kinds of individual neuronal connection pattern: one in which each neuron projects to most or all of its possible lattice neighborhoods; or one in which each neuron projects to only a single or very few neighborhoods. These organizations imply very different computational architectures. The first, shared output, corresponds to a degenerate or redundant architecture in which a given neuron broadcasts its output to neighboring lattice regions. In contrast, a segregated output suggests selective routing of output information. Our results support this second, selective routing, case.
Conclusions
The functional role for the clusters of boutons is still speculative, but the results of this study indicate that the clusters form a network with connectivity structure that is less regular than is usually associated with the visual cortex. Nearby neurons may route signals to entirely different, well separated sites, allowing for information to be distributed and mixed without much redundancy and correlation. How such a heterogenous network leads to the functional architecture of the visual cortex is the next key problem.
Footnotes

This work was supported by European Union Grant FP6–2005015803. T.B. has a Research Councils United Kingdom Academic Fellowship.
 Correspondence should be addressed to Tom Binzegger, School of Biology, Henry Wellcome Building for Neuroecology, University of Newcastle upon Tyne NE2 4HH, UK. tom.binzegger{at}ncl.ac.uk