Neural Circuits Original Research Article Materials and Methods Preparation of Brain Slices

clature for morphological, physiological and molecular features was agreed upon by the Petilla Interneuron Nomenclature Group (Ascoli et al., 2008). Among GABAergic interneurons, those expressing the neu-ropeptide somatostatin (SOM) are particularly heterogeneous in their molecular, morphological and electrophysiological features (Sabo and Sceniak, 2006). To better understand this variability, we took advantage of a transgenic mouse line where SOM-positive neurons are labeled with GFP (Oliva et al., 2000) to explore if sub-types of SOM interneurons could be distinguished objectively. 59 SOM-positive interneurons from mouse primary somatosensory, frontal and visual cortex were quantitatively characterized using their morphological properties measured from reconstructions of biocytin-fi lled cells and their intrinsic electrophysiological properties , measured from whole cell recordings. Unsupervised cluster analysis revealed a group comprised of the well-known Martinotti cells, as well as two other subgroups of SOM-positive neurons that, to our knowledge, have not yet been described. The electrophysi-ological classifi cation agreed well with, and thus confi rmed, the morphological classifi cation. Furthermore, signifi cant differences in axon morphology and fi ring properties between the three groups suggest that the two novel subtypes and Martinotti cells could serve different functional roles in the cortical circuit. Acute brain slices were prepared from GIN mice (Oliva et al., 2000), with an average age postnatal 14 days (range P10–18). Mice were quickly decapitated, the skin and skull were cut and the brain was removed and then immediately placed in cold sucrose cutting solution (222 mM sucrose, 2.6 mM KCl, 27 mM NaHCO 3 , 1.5 mM NaH 2 PO 4 ,

Interneurons exhibit very diverse morphology, electrophysiology, molecular content and post synaptic targets (Cauli et al., 1997;Markram et al., 2004;Yuste, 2005). Classifi cation of neocortical interneurons is a crucial step in understanding cortical circuits as each subtype of interneuron likely has a different function. In the past, classifi cations have been based on one or a combination of descriptors, often using qualitative criteria which were not standardized. Therefore, interneuron classifi cations often differ and it is unclear which set of descriptors are the most relevant to determine a neuronal class and, more generally, how many classes of interneurons actually exist. As a fi rst step, several groups (Cauli et al., 2000;Karube et al., 2004;Dumitriu et al., 2006;Ma et al., 2006;Helmstaedter et al., 2008;Karagiannis et al., 2009) have used quantitative methods, with unsupervised clustering algorithms to seek an objective classifi cation of interneurons. In addition, to facilitate the classifi cation of interneurons, a standardized nomen-2 mM CaCl 2 , 2 mM MgSO 4 , bubbled with 95% 0 2 , 5%CO 2 ). Slices 300-µm thick were cut using a Vibratome and transferred to a chamber at room temperature with oxygenated ACSF (126 mM NaCl, 3 mM KCl, 2 mM MgSO 4 , 2 mM CaCl 2 , 1 mM NaH 2 PO 4 , 26 mM NaHCO 3 and 10 mM glucose, bubbled with 95% 0 2 , 5%CO 2 ).

MICE STRAIN
GIN mice reliably label a subset of SOM-positive cells with expression remaining consistent for more than fi ve generations (Oliva et al., 2000). Two recent studies have repeated immunohistochemical analysis of this original paper, demonstrating that GFP expressing cells are indeed SOM-positive. One group reported colocalization in 95.9% of layer 2/3 cells, 100% of layer 4 cells and 97.4% of layer 5/6 cells (Ma et al., 2006) and the other reported co-localization of 99% for cells in all layers (Halabisky et al., 2006). In this study, the percentage of SOM + cells that were GFP labeled was 34.8% in layer 2/3, 26.5% in layer 4 and 10.8% in layer 5/6. To our knowledge no one line of transgenic mice labels all SOM + cells and the most often used GIN line labels a higher percentage of SOM + cells than other lines (Ma et al., 2006).

ELECTROPHYSIOLOGY RECORDINGS
Slices were placed in a recording chamber at room temperature with fl owing oxygenated ACSF. Pipettes of 3-7 MΩ resistance were pulled from borosilicate glass. SOM positive cells were identifi ed by expression of GFP. Whole cell recordings were taken in current clamp mode. Only cells with healthy resting membrane potential (between −55 and −80 mV) were selected for recording.

Electrophysiology analysis
Ninteen variables were measured for each neuron by analysis of the recordings in MATLAB. See Table 1 for descriptions.

Histological procedure
Neurons were fi lled with biocytin by a patch pipette. The slices were kept overnight in 4% paraformaldehyde in 0.1 M phosphate buffer (PB) at 4ºC. The slices were then rinsed three times for 5 min per rinse on a shaker in 0.1 M PB. They were placed in 30% sucrose mixture (30 g sucrose dissolved in 50 ml ddH 2 0 and 50 ml 0.24 M PB per 100 ml) for 2 h and then frozen on dry ice in tissue freezing medium. The slices were kept overnight in a −80ºC freezer. The slices were defrosted and the tissue freezing medium was removed by three 20-min rinses in 0.1 M PB while on a shaker. The slices were kept in 1% hydrogen peroxide in 0.1 M PB for 30 min on the shaker to pretreat the tissue, then were rinsed twice in 0.02 M potassium phosphate saline (KPBS) for 20 min on the shaker. The slices were then kept overnight on the shaker in Avidin-Biotin-Peroxidase Complex. The slices were then rinsed three times in 0.02 M KPBS for 20 min each on the shaker. Each slice was then placed in DAB (0.7 mg/ml 3,3′-diaminobenzidine, 0.2 mg/ml urea hydrogen peroxide, 0.0 6M Tris buffer in 0.02 M KPBS) until the slice turned light brown then immediately transferred to 0.02 M KPBS and transferred again to fresh 0.02 M KPBS after a few minutes. The stained slices were rinsed a fi nal time in 0.02 M KPBS for 20 min on a shaker. Each slice was observed under a light microscope and then mounted onto a slide using crystal mount.

Reconstruction and analysis of morphology
Successfully fi lled and stained neurons were then reconstructed using Neurolucida software (MicroBrightField). The neurons were viewed with 100× oil objective on an Olympus IX71 inverted light microscope or an Olympus BX51 upright light microscope. The Neurolucida program projects the microscope image onto a computer drawing tablet. The neuron's processes were traced manually while the program recorded the coordinates of the tracing to create a digital three dimensional reconstruction. The x and y axis form the horizontal plane of the slice, while the z axis is the depth. The user defi ned an initial reference point for each tracing. The z coordinate was then determined by adjustment of the focus. In addition to the neuron, the pia and white matter were drawn.
The Neurolucida Explorer program was used to measure 67 morphological variables of the reconstruction. See Table 2 for descriptions. Table 1 | Electrophysiological variables. Action potential properties measured from response to twice threshold, 500-ms current injection from fi rst action potential (AP1) and second action potential (AP2). AP2 variables not listed as the same measurements were made for AP2 as listed for AP1.

Variable Description
Resting membrane potential (mV) Stable membrane potential when no current applied

VARIABLES DESCRIBING LOCATION
Relative distance to pia Distance from soma centroid to pia/distance between pia and white matter Shrinkage percentage in the z-axis was on average 48.97 ± 1.95%. No shrinkage correction was applied.

Principal component analysis and cluster analysis
The dataset of neurons could now be represented as a matrix, with each row corresponding to a neuron and each column corresponding to a variable. The data was then standardized, where the standardized values of each variable are computed as the difference between the actual value and the mean divided by the standard deviation. PCA reduces the dimensionality of the dataset, while preserving as much of the variance as possible by mapping the coordinate system defi ned by the old variables to a new lower dimensional coordinate system defi ned by principal components that are orthogonal, thus uncorrelated, linear combinations of the original variables. As many of the original variables are related, the principal components produced by PCA reduce this redundancy. By PCA, a new space is generated onto which the dataset can be projected and classifi ed by cluster analysis.
The principal components were calculated using the factor analysis function in STATISTICA (StatSoft). The principal components are found by an eigenvalue decomposition of the correlation matrix of the standardized data. The eigenvectors of the correlation matrix are the principal components (PCs) and the corresponding eigenvalues give the variance preserved by each principal component. The number of PCs maintaining a significant amount of variance will always be less than the number of original variables, which is why PCA reduces the dimensionality of the dataset (Jambu, 1991).
Several criteria were used to decide how many principal components to retain for cluster analysis. First, only principal components with eigenvalues greater than one are considered. The original variables have an eigenvalue of one, as the data is standardized, thus any principal component with an eigenvalue greater than one describes more of the data's variance than an original variable. Second, the "scree test" was used (Cattell, 1966). A plot of the eigenvalues of the principal components is examined for when the decrease in eigenvalue plateaus. In all cases either 2 or 3 principal components were retained.
Cluster Analysis is then performed on each dataset in the principal component vector space using Clustan (Clustan Ltd.). Hierarchical clustering using Euclidean distance squared as the multi-dimensional distance metric and Ward's method as the linkage rule was performed. Ward's method minimizes the distance squared between any two clusters that can be formed at each step of clustering (Ward, 1963). Alternatively, cluster analysis was performed on the standardized data before PCA. The results obtained from this method did not dramatically differ from the results obtained with PCA. The only differences were reordering of cells within clusters and movement of cells between neighboring clusters contained in Group 1.

Statistical analysis of clusters
Three statistical tests were done to fi nd the clustering levels which have signifi cance at the 5% level. The Best Cut Test does an analysis of variance on the fusion values (distance at which two clusters join) at every level in the dendrogram. The realized deviates, defi ned as the standardized difference between a fusion value and the mean fusion value, are compared to fi nd signifi cantly large realized deviates. A large realized deviate is signifi cant if the fusion value is at least 1.96 standard deviations from the mean. The upper tail test applies the upper tailed t-test to the fusion values using t-statistic of realized deviate*(n-1) over n-2 degrees of freedom to fi nd the signifi cant increases in fusion value. The bootstrap test performs 1000 trials of randomizing the data, cluster analysis on the randomized data and then compares the random trees to the actual tree. The data matrix was randomized by randomly permuting the order of entries in each column independently. This scrambles the values of each variable for a given cell, thus disturbing any correlations between variables for each cell, but not changing the distribution of values for each variable. The random trees were compared to the actual tree by plotting fusion value vs. number of clusters for the actual data and the distribution of the randomized data with a confidence interval of 1 standard deviation around the mean. Signifi cant departures of the actual fusion values from random, determined by the t-statistic, indicate a number of clusters that are statistically signifi cant. The fi rst level at which signifi cance was found by all three tests is the partition used (indicated in Figures 1A, 2A and 3C-E by bracketing on the bottom of the dendrogram and cut-off linkage distance by a horizontal line).
Finally, in the morphological or physiological clustering we did not detect any obvious bias towards the individual patching the cell, individual reconstructing the cell, or age of the mouse.

K-means clustering
K-means clustering was performed for each dataset using Clustan's Focal Point clustering algorithm, with K equal to the number of signifi cant clusters from hierarchical clustering. Focal Point runs the K-means algorithm several times with different initial case orders, which is an important permutation as the K-means algorithm is sensitive to case order. It ranks these solutions based on the minimization of Euclidean distance squared.

SILHOUETTE ANALYSIS
Silhouette analysis measures the quality of the clustering by examining the within cluster distances and between cluster distances (Rousseeuw, 1987). The silhouette value of a data point i in cluster A is where a(i) is the average Euclidean distance between i and all other datapoint in cluster A, b(i) is the smallest average Euclidean distance between i and all datapoints in any other clusters (the average Euclidean distance between i and all datapoints in cluster B, the nearest neighbor to cluster A). The value of s(i) is between −1 and 1. The silhouette width of a clustering is defi ned as the average of s(i) for all datapoints. Large positive values indicate that the cluster is compact and distinct from other clusters. In addition to calculation of the silhouette width of each clustering, the silhouette width of randomized clustering for each dataset was calculated. The data was randomized as described above for the bootstrap test. The same procedure for cluster analysis was performed on 500 randomized matrices of each dataset and the silhouette width of each clustering was computed. S(random) is the average of the 500 silhouette widths.

Statistical analysis of variables
To determine which variables are important in distinguishing between the groups identifi ed by cluster analysis, the group mean and standard error was calculated for every variable and the groups were compared for signifi cant differences (Mann Whitney U Test). Additionally, the correlation between the PCs and the original variables was calculated as another measure of which variables are discriminating. Variables highly correlated with the principal components used for cluster analysis are more discriminating. Most of the variables with signifi cant differences between the groups were highly correlated with the PCs (Tables 6 and 7)

MORPHOLOGICAL CLASSIFICATION OF SUBTYPES OF SOMATOSTATIN INTERNEURONS
To identify subtypes of SOM-positive neurons, morphology and electrophysiology were quantitatively characterized. GFP-positive neurons were patched, recorded from in current clamp mode, fi lled with biocytin and processed for morphological analysis. 19 electrophysiological variables were measured from current clamp recordings. In addition, 67 morphological variables describing the soma, dendrites, axon and somatic location (relative to the pial surface) were measured from the Neurolucida reconstructions of We fi rst explored the statistical structure of the morphological dataset. Cluster analysis using morphological variables was performed with the 39 cells with complete reconstructions, and it revealed fi ve clusters of neurons, grouped into three major branches ( Figure 1A, labeled purple, orange and green). One major group encompassed approximately half of the neurons, and could be subdivided into three clusters (purple cells, clusters a, b and c), and was separated by a large Euclidian distance from the other two clusters (orange-d and green-e). In addition, clusters showed no bias with respect to age of the mouse, cortical area, individual patching the cell, or individual reconstructing the cell (not shown).
We further explored the morphological variability of the neurons by plotting their position in principal component space, using the fi rst two principal components as Cartesian axes ( Figure 1B). In this analysis one seeks to draw lines that separate each group in a different section of the graph. In this PCA space, the fi rst major group of the cluster analysis (purple cells, corresponding to clusters a, b and c) was clearly separated from the other two branches (orange-cluster d and green-cluster e cells). Within this fi rst group (purple cells), no clear subgroups were observed in PCA space. Clusters d (orange) and e (green) were separated from each other by a second line, with one exception (orange neuron on top). Because of these results in PCA space, although the morphological grouping was signifi cant at the 5 cluster level (black brackets in Figure 1A; see Materials and Method and Statistical Validation section below), we chose to group the data into three classes (named groups 1, 2 and 3, and colored purple, orange and green, for the rest of the study).  The fi rst 2 principal components were retained for cluster analysis by the morphology variables (C) and electrophysiological variables (D). The fi rst 3 principal components were retained for cluster analysis by both the electrophysiology and morphology variables (E). Again the color scheme is based on the clustering by morphology (see Figure 1).
Examining the morphologies of the cells revealed that group 1 (n = 21) corresponded to morphologies traditionally associated with Martinotti cells, a well-characterized subtype of SOM interneurons (Wahle, 1993;Kawaguchi and Kubota, 1996;De Felipe, 2002;Wang et al., 2004). Martinotti cells are found in layers 2/3, 5 and 6 and have a characteristic ascending axon that sends several collaterals to layer 2/3 and layer 1. The axon branches extensively in these layers, particularly in layer 1 where the collaterals branch horizontally, for distances as long as 400 µm (see Figure 1C, for representative examples; Figure S1A in Supplementary Material, for all neurons) On the other hand, group 2 (n = 11) and 3 (n = 7) neurons were quite distinct from Martinotti cells, with very different axonal morphologies. In group 2 and 3 cells, the axon ascended in a very direct course with minimal branching and, in most cases, did not enter layer 1 (Figure 1C and Figures S1B,C in Supplementary Material). Axons generally had only one main ascending collateral, or at most, two or three. Initially, we interpreted this sparse axon as a fi lling or reconstruction artifact, but closer examination of the cells found no evidence of incomplete staining or cut processes. Indeed, endings of their axonal processes did not taper and had apparently normal terminations in the middle of the section (see below). Moreover, group 2 and 3 neurons all had healthy-looking dendrites, without beading, with group 2 cells displaying multipolar morphology and group 3 cells bipolar morphology. Group 2 and 3 cells were found in layers 2/3, 4 and 5.
The combined cluster/PCA analysis indicated that the morphological variance of our database of neurons could be well captured by the hypothesis that they belonged to three major subtypes of cells.

PHYSIOLOGICAL CLASSIFICATION OF SUBTYPES OF SOMATOSTATIN INTERNEURONS
We then sought to use the electrophysiological database to test whether the groups defi ned by the morphological variables were correlated to electrophysiological groups. For this analysis we used neurons with complete physiological data, regardless of their anatomical reconstruction. Therefore, the two datasets, anatomical and physiological, were not identical but partly overlapped.
When the 36 cells with recordings were clustered by electrophysiology variables, we found signifi cant grouping at the 5 cluster level, by several independent statistical tests (Figure 2A; black brackets; clusters a to e, all neurons with reconstructions colored according to morphological groups described above). The fi rst two clusters (a and b) constituted a group (group 1, purple), in which, remarkably, all of its reconstructed neurons had been previously identifi ed as belonging to the group 1 from the morphological clustering ( Figure 1A, purple cells). Besides this fi rst group, the remaining 3 clusters (c-e) contain all neurons previously identifi ed as belonging to the morphological groups 2 and 3 and, moreover, separated these cells in good agreement with these 2 groups. More specifi cally, physiological cluster c corresponded to morphological group 3 (green) and physiological clusters d and e correspond to morphological group 2 (orange) with one outlier neuron in cluster e that belonged to the fi rst morphological group (purple).
In addition to the cluster analysis, we again explored the statistical structure of the physiological data by plotting the dataset in principal component space ( Figure 2B). This PCA analysis revealed a less compact grouping than the morphological dataset. Nevertheless, with two orthogonal lines, one could clearly separate neurons belonging to groups 1, 2 and 3, again, with only one exception, the previously mentioned outlier neuron from cluster e ( Figure 2B, purple neuron indicated by arrow). This neuron, interestingly, was close in PCA space to the border between groups 1 and 2.
Inspection of the physiological responses of these three groups of neurons was in agreement with the statistical clustering. Group 1 cells (n = 24) displayed a regular spiking fi ring pattern, nonstuttering and non-fast spiking, with varying degrees of spike frequency adaptation from small to moderate increases in interspike interval (Figure 2C, a and b; average spike frequency adaptation 1.674 ± 0.092; see Ascoli et al., 2008 for defi nition of electrophysiological variables according to the Petilla nomenclature). For the majority of neurons the action potential (AP) amplitude decayed only slightly over the course of a train (average AP drop 2.619 ± 0.479 mV). Group 1 cells also tended to have a lower rheobase (36.042 ± 6.363 pA) than group 2 or 3 cells (49.286 ± 10.025 pA and 55.000 ± 12.845 pA, respectively) and a signifi cantly more depolarized resting membrane potential ( Table 3). On the other hand, group 2 cells (n = 7) displayed both regular spiking and stuttering fi ring patterns. There was a greater spike frequency adaptation (2.244 ± 0.385) and amplitude accommodation than group 1 cells. Group 3 cells (n = 5) were regular spiking, with a degree of frequency adaptation similar to group 2 cells (2.657 ± 0.493) and signifi cantly narrower action potentials ( Table 3, AP1 half-width, AP2 half-width). One cell (cell 3, see Figure S2 in Supplementary Material) had an unusual early offset fi ring pattern of only three to four action potentials per train, and then was silent at a depolarized membrane potential for the remainder of the current pulse. Of the four other group 3 cells, three showed marked amplitude drops in the fi rst few spikes and the other showed amplitude accommodation similar to group 2 cells ( Figure S2).
The joint interpretation of the morphological and physiological clustering and PCA analysis was consistent with a simple model of three basic subgroups of SOM neurons. Group 1 is clear: it is composed of Martinotti cells and is clearly delineated in the morphological and physiological clustering dendrograms. Group 2 and group 3 are also distinct in both morphological and physiological clustering, although they are closer to each other in statistical space. Specifi cally, the close proximity and some overlap between groups 2 and 3 in the principal component space of the morphology dataset versus their complete separation in the principal component space of the electrophysiology dataset suggests that groups 2 and 3 have morphological similarities, although they are physiologically quite different. One could further subdivide the morphological group 1 into three different subgroups (a, b, c; based on the morphological clustering), or the physiological group 3 into 2 subgroups (d and e; based on the electrophysiological clustering). Nevertheless, for simplicity, for the rest of the study we followed the basic three groups (1, 2 and 3), fi rst apparent from the morphological clustering.
These results display a remarkably good agreement between the morphological and physiological classifi cations, which because they are based on completely independent datasets, confi rm each other. In fact, with this basic classifi cation into three groups every neuron is classifi ed into the same groups in both morphological and physiological analysis, with one exception, the outlier purple cell that clusters with group 3 in the physiology, but which is actually borderline in PCA space. Given the tight correlation between morphological group and electrophysiological group for the 16 cells that are in both morphology and electrophysiology datasets, we propose that morphological groups 1, 2 and 3 correspond to electrophysiological groups 1, 2, 3. We do recognize that for a cell belonging to only one of the datasets it is possible that it does not belong to the corresponding group of the other dataset. However, we fi nd this unlikely given that for 15/16 cells the correspondence between morphological and electrophysiological groups was accurate.

STATISTICAL VALIDATION OF CLUSTERS
We then explored how robust the classifi cation based on cluster analysis actually was, using a variety of statistical analyses. We fi rst used three tests (best cut, upper tail and bootstrap) that are routinely applied to test whether the dataset has a non-uniform structure, an important validation as cluster analysis will always return a clustering tree, regardless of the distribution of the data. Specifi cally, the best cut test distinguishes the clustering profi le of a dataset containing groups from the profi le of a dataset with a continuum of points by analyzing the fusion values. Cluster analysis performed on a dataset not containing groups would return arbitrary divisions and thus the fusion values would be continuous, with no large jumps. However, if the dataset does contain groups then at the partitions between the groups large fusion values should be seen. The best cut test fi nds if there is a linkage cut-off distance at which this fusion value profi le occurs. The upper tail test also locates signifi cantly large fusion values. The bootstrap test determines the linkage distance at which the results signifi cantly deviate from random. Indeed, using these three tests, the clusters (as marked by brackets in Figures 1A and 2A) were found to be statistically signifi cant at 5% level (See Materials and Methods for details). Thus, three independent signifi cance tests all show that the dataset does have a non-uniform structure, and thus could be classifi ed into meaningful clusters. Following this, we sought to confi rm whether the hierarchical classifi cation using Ward's method was reliable, using K-means clustering with K specifi ed as the number of signifi cant clusters from hierarchical clustering. K-means clustering is a top-down clustering method, in contrast to hierarchical clustering which is a bottom-up method. (MacQueen, 1967). While hierarchical classifi cation is infl exible in that cells initially linked stay linked, resulting in a sometimes local optimization rather than global optimization, K-means can reassign cells at a later stage of clustering if it will result in a better optimization. Indeed, the reliability of Ward's method hierarchical clustering for this dataset was proven by the agreement between the two methods ( Table 4). Only a small percentage of cells switched clusters and always stayed within the same group, for example in the 39 morphology cluster cell 27 switched from cluster a to cluster b, both of which are included in Group 1 ( Table S3 in Supplementary Material). Such stability between two different algorithms demonstrates that the groups are robust.
In addition, to assess the quality of the hierarchical clusters obtained by Ward's method, we computed the silhouette width of each clustering (Figures 3A,B) (Rousseeuw, 1987, See Materials and Methods). Silhouette width is a measure of the appropriateness of each assignment of a datapoint to a cluster. A silhouette width around 0 indicates the datapoint is intermediate between two clusters, silhouette width close to 1 indicates that the point is well assigned and a silhouette width close to -1 indicates that the point is poorly assigned. Constructing a plot of silhouette widths helps reveal the structure of the dataset, and thus to identify natural groups from artifi cially imposed groups. Using this analysis, we found that the silhouette widths of the clustering by morphology and clustering by electrophysiology confi rmed that a natural structure has been found. Such a silhouette width indicates that the distances within clusters are small in comparison to the distances between clusters. Thus, the dataset has a structure containing discrete and compact clusters with clear separation between clusters. In fact, visually inspecting the silhouette plot for each cluster revealed that the majority of cells are well classifi ed.
As a fi nal additional measure of statistical signifi cance, the silhouette width was computed for 500 trials of randomized data, following the same clustering procedure ( Table 5). The average silhouette widths for these randomized trials was drastically lower than for the actual datasets, again demonstrating that the datasets have a strong structure which was disrupted by the randomization.

QUALITY OF CLUSTERING COMBINING DIFFERENT GROUPS OF VARIABLES
To determine the quality of classifi cation based on electrophysiological variables, morphological variables and all variables combined, we performed a separate analysis of the 16 cells with both complete reconstructions and recordings (Figures 3C-E). This dataset was used to control for differences in sample size and nonoverlapping data points between the 36 cell electrophysiology and 39 cell morphology datasets. The clusters found by morphology only, electrophysiology only and both morphology and electrophysiology were in good agreement with the three groups previously identifi ed. Each set of variables had one outlier in the clustering therefore clustering using the combined sets of variables was not an entirely better method. Judging the strength of structure of the three datasets by silhouette width revealed that the clustering by morphological variables has the strongest structure, by electrophysiological variables the next strongest, and by all variables the weakest. The low silhouette width of the combined morphological and electrophysiological variables dataset could be due to the strong correlations within clusters between morphological variables and between electrophysiological variables. With all variables combined, the clusters had a less compact shape due to the multimodality of variables. Thus when clustering by electrophysiology variables only or morphological variables only can discriminate between groups as in our case, we found no particular advantage to clustering by morphological and electrophysiological variables combined.

MORPHOLOGICAL CHARACTERISTICS DEFINING SOM SUBGROUPS
After establishing the reliability and validity of the classifi cation based on three major groups of SOM interneurons, we explored more carefully the physiological and morphological characteristics of each of them. First, we concentrated on the analysis of the morphological variables that were statistical predictors of each of the three subtypes of neurons (Table 3). Morphologically, the axon was the most important feature in discriminating between the group 1 and groups 2 and 3, as can be seen by direct inspection of the morphologies. Group 1 cells had typical Martinotti cell axons that extend through layers above the soma and branch  extensively, particularly at the pial surface and in the home layer. Group 2 and 3 cells also had an ascending axon, but their axons were much sparser and were contained within a narrow vertical column. Statistically, Group 1 cells had principal component 1 (PC1) values distinct from groups 2 and 3 ( Figure 1B) and, indeed, the fi rst principal component of the dataset is highly correlated with 12 axon variables ( Table 6). Six of these variables, axonal node total, total axonal length, total surface area of axon, convex hull axon area, convex hull axon perimeter and convex hull axon surface area, describe the size of the axon in both two and three dimensions. The other variables, i.e. highest order axon segment, axonal Sholl length and node densities and axonal node density describe the degree of axon branching and k-dim axon describes the dimensionality of the axon. Given these correlations PC 1 can be interpreted as a measure of axon size and complexity. Accordingly, there were signifi cant differences for several parameters describing the size, shape and dimensionality of the axon between groups 2 and 3 and group 1 ( Table 3). The layers targeted by the axon also differed between group 1 and groups 2 and 3. Group 1 cells extensively branched in layer 1, where Martinotti cells are known to synapse on the apical dendrite tufts of pyramidal cells Silberberg and Markram, 2007). In contrast, 13 of the 18 group 2 and 3 cells (72.2%) had axons that avoided layer 1. Some had axons that suddenly terminated at layer 1, while others made an abrupt turn at the layer 1 border and either hooked or continue tangentially along the border ( Figure 4A). Interestingly, there was a striking direction preference of this bend at the layer 1 border; in every single case axons turned medially in the slice. Examination of the axons of these cells showed axonal boutons in layer 2/3 ( Figure 4B). As mentioned, these axons did not taper and were not cut, terminating in an apparently normal fashion, away from the sectioned surface of the slice (Figure S3B in Supplementary Material). In contrast, the termination of cut axons had larger, more roundly swollen ends and darker ends, and always ended at the very edge of the slice (Figure S3C in Supplementary Material). The size and shape of the dendritic arbor anatomically distinguished groups 2 and 3, which had overall similar axon morphology. While group 1 was separated from groups 2 and 3 on the PC1 axis, groups 2 and 3 were separated from each other on the PC2 axis ( Figure 1B). Examining correlations between the original variables and PC2, revealed that PC2 was highly correlated with several dendrite variables describing the size of the dendrites in two and three dimensions ( Table 6). This possibility is supported by signifi cant differences between the two groups for several variables describing size of the dendritic processes (Table 3). Group 3 cells had the largest dendritic arbor of the three groups, while group 2 cells had the smallest dendritic arbor, by measures of length and convex hull size. Additionally, group 3 cells had bitufted dendrites, while group 2 cells have multipolar dendrites (Figure S1 in Supplementary Material). Group 1 dendrites were varied, with some multipolar, while others show a preference for extending only upward or downward.

ELECTROPHYSIOLOGICAL CHARACTERISTICS DEFINING SOM SUBGROUPS
We then focused on the electrophysiological parameters, analyzing both passive and active physiological properties, such as resting membrane state and fi ring pattern (Table 2, Figure 5). There were signifi cant differences among the passive properties of the three groups ( Figure 5A). Group 2 and 3 cells were signifi cantly (p < 0.025) more hyperpolarized at rest (−70.617 ± 2.385 mV and −70.236 ± 2.228 mV respectively) than Group 1 cells (−65.336 ± 0.438 mV). Group 3 cells had signifi cantly (p < 0.025) lower input resistance (256.400 ± 46.427 MΩ) than both group 1 and 2 cells (469.583 ± 35.836 MΩ and 475.429 ± 57.726 MΩ respectively). In response to hyperpolarizing current, only group 1 cells consistently demonstrated a sag response. This variable was not used for clustering as only 32 of the 36 cells (23 of 24 group 1, 5 of 7 group 2 and 4 of 5 group 3 cells) were given a series of hyperpolarizing current steps. The slope of the regression line fi tted to the sag vs. membrane potential was used as the variable for comparison, with sag defi ned as the difference between the most negative membrane potential during a 1s hyperpolarizing current step and the steady membrane potential in the last 100ms of the step. Group 1 signifi cantly differed (p < 0.005 and p < 0.01 respectively) from groups 2 and 3 ( Table 3). Some cells showed a depolarized rebound after the hyperpolarizing step, but only one cell in group 1 spiked as a result of the rebound depolarization. As Martinotti cells are known to display a sag response, these results further support that groups 2 and 3 are distinct from Martinotti cells.
We found that the fi rst principal component was highly correlated with the variables measuring the duration, half-width, rise time and fall time of AP1 and AP2 (Table 7). Therefore it is not surprising that these variables distinguished all three groups ( Figure 5B). Specifi cally, the groups had distinct AP duration, with the values of group 1 (AP1 duration 3.924 ± 0.198 ms) being intermediate between groups 3 (AP1 duration 3.000 ± 0.202 ms, p < 0.05) and 2 (AP1 duration 5.914 ± 0.734 ms, p < 0.025). In addition, groups 1 and 3 both had faster half-width (AP1 halfwidth 1.649 ± 0.083 ms and p < 0.005, 1.356 ± 0.071 ms respectively) rise time (AP1 rise time p < 0.025, 1.070 ± 0.036 ms and FIGURE 5 | Analysis of electrophysiology variables. Comparison of physiological variables that distinguish the three groups with regards to several variables. All data shown mean ± standard error (A) Groups 2 and 3 are distinct from group 1 with respect to RMP, input resistance, AP amplitude and difference in amplitude between fi rst and second AP of a train. (B) Each group has a distinct time course of AP, with group 1 intermediate between groups 2 and 3. (*p ≤ 0.05, **p ≤ 0.025, ***p ≤ 0.005, Mann Whitney U Test) p < 0.005 0.880 ± 0.037 respectively), and fall time (AP1 fall time p < 0.025, 2.853 ± 0.168 ms and p < 0.005, 2.120 ± 0.166 ms respectively) than group 2 (AP1 half width 2.271 ± 0.291 ms, AP1 rise time 1.400 ± 0.102 ms and AP1 fall time 4.514 ± 0.656 ms) while groups 1 and 2 had signifi cantly slower (p < 0.025) rise rates (AP1 rise rate 54.265 ± 3.733 mV/ms and 48.079 ± 9.906 mV/ms respectively) and fall rates (AP1 fall rate 21.812 ± 1.891 mV/ms and 17.304 ± 4.749 mV/ms respectively) than group 3 (AP1 rise rate 80.337 ± 4.179 mV/ms, AP1 fall rate 33.921 ± 2.836 mV/ms). The slower rise and fall rate of group 1 was more likely due to group 1 having smaller AP amplitude, while the difference in rise and fall rate between group 2 and 3 resulted from the difference in AP rise and fall time.
Most cells were regular spiking, but three group 2 cells (75, 76, 77) were stuttering, defi ned as the fi ring of groups of action potentials that are interrupted by silent periods of varying length Ascoli et al., 2008). With only an 8.3% incidence of stuttering cells in the dataset, this is too low a number to conclude that all stuttering cells are group 2 cells, but does indicate that stuttering is a distinctive fi ring pattern of some group 2 cells. It can additionally be noted that group 2 cells had greater variability in fi ring patterns, by plotting the dataset in the PC1-PC2 plane ( Figure 2B). Groups 1 and 3 were restricted to a much smaller region of the PC1 values, than are Group 2 cells indicating that Group 2 cells have a larger and separate range of values for those variables correlated with PC1 (AP duration, half-width, rise time and fall time).

THREE SUBTYPES OF SOMATOSTATIN-POSITIVE INTERNEURONS
In this study we used unsupervised cluster analysis to quantitatively classify SOM interneurons based on morphological and electrophysiological characteristics. Using different methods, we fi nd a statistically consistent structure which can be captured well by a simple model of three basic subtypes of cells, already apparent in the fi rst morphological clustering tree ( Figure 1A). This does not mean that there are not further subtypes of SOM neurons, but that with the size of our database and methods used in this study, and within the GFP expression pattern of this particular mouse line, at least these three distinct subtypes are apparent. Specifi cally, the good agreement between clustering based on two independent sets of variables, electrophysiological and morphological, provided strong validation that the groups found by cluster analysis are meaningful and distinct.
Because of the sparseness of the group 2 and 3 cells' axons we suspected that their axonal morphology was due to incomplete fi lling or sectioned processes. However, as discussed, we did not fi nd clear evidence of cut axonal terminations in these neurons. While we also did not observe in group 2 and 3 cells any clear signs of incomplete fi lling, such as dye blown out of the cells or tapering and spotty fi lling at the end of visible processes, it is still possible that these cells were incompletely fi lled without these tell tale signs. Another possibility is that these axons were projecting away from the region and we only revealed its local portion. As with other neuronal classes, the projecting axon could be myelinated, and may be diffi cult to visualize in its entirety. In any case, given that our data was collected from brain slices, we are the fi rst to admit that in vivo evidence is required to conclusively demonstrate the true morphological nature of the axons from groups 2 and 3 cells. Nevertheless, whether their axons are complete or not does not negate our conclusion that they represent different cell types.
Specifi cally, there are several pieces of evidence in support of group 2 and 3 being real subtypes. First, the axons of group 2 and 3 neurons bend medially, something diffi cult to reconcile with the possibility that they represented cut group 1 cells, whose axons do not show a characteristic medial bend.
Second, some of the statistical differences in the axonal structures of group 2 and 3 vs. 1 cannot be due to sectioned processes or incomplete fi lling. For example, differences in the k-dim of axons, or the standard deviation of the tortuosity of axonal nodes, or the axon node densities are not reliant on the size or length of the axon, but instead measure degree of branching and space fi lling qualities of the axon. If group 2 and 3 cells were the same as group 1, but with cut axons, one would expect these values not to be signifi cantly different across the three groups, because they refl ect differences in the topology of the axonal arbor other than size. To test this directly we performed an additional cluster analysis. We excluded all variables related to the size or length of the axon, recalculated principal components and performed hierarchical cluster analysis using these principal components. The following variables were excluded: axonal node total, total axonal length, total surface are of axon, ratio of axonal length to surface area, highest order axon segment, all axonal Sholl variables (6), and all axonal convex hull variables (4). Based on the three statistical tests, the results are signifi cant at the 6 cluster level, and agree well with our results using all morphological variables. There are two outliers, both cells from group 2. One clustered with group 1, and the other with group 3 (Figure S3 in Supplementary Material). Third, we excluded four cells from our dataset because 3 had cut processes and 1 was not completely fi lled. These cells resembled Martinotti cells with the ascending portion of the axon cut. Analysis of their axon parameters agreed with this interpretation. The cut Martinotti cells had a mean total axon length signifi cantly larger than groups 2 and 3, but slightly smaller than group 1 (p < 0.005 for group 2 and p < 0.01 for group 3). These cut cells were also signifi cantly different from group 2 (p < 0.005) and group 3 (p < 0.01) with respect to axonal node total, total axonal length, total surface area of axon, and was not signifi cantly different from group 1 with respect to the same variables. The cut Martinotti cells also had the same axon branching properties as group 1 cells, highly branched and space fi lling, distinct from the sparse axon branching of groups 2 and 3 (again, not signifi cantly different from group 1 with respect to axonal node total, k-dim axon, stdev of tortuosity of axonal node, axon node density, p < 0.025 for groups 2 and 3 with respect to those same variables).
Finally, in some slices more than one GFP cell was patched, fi lled and reconstructed. In most cases, all cells were Martinotti neurons. However in two slices there was both a Martinotti cell and a group 2 cell (cell 28 and 29 together and cell 52 and 53 together). Since the axon of the Martinotti cell was large and complete, and the cells were located in a similar section of the slice, it is unlikely that only the axon of the group 2 cell was selectively cut. Moreover, cells 52 and 53 even had overlapping axonal arbors in some portions and there was not evidence of incomplete fi lling.
Further evidence that group 2 and 3 are real comes from their electrophysiological properties. Specifi cally, clustering by electrophysiology agreed with clustering by morphology as both separated the Martinotti, group 2 and group 3 cells, providing an independent control for the possibility that the group 2 and 3 cells were poorly fi lled or reconstructed group 1 cells. The AP durations of the cells are long but consistent with literature values for recordings from juvenile mouse slices at room temperature. We attribute these long APs observed to two reasons. Firstly, juvenile animals have slower action potentials than adults and secondly APs are slower at room temperature than at 35-37ºC (Ali et al. 2007). SOM + cells in P13-P16 rat slices recorded at room temperature had AP1 duration of 3.69 ± 0.53 ms and AP2 duration of 4.32 ± 0.55 ms, similar to our Martinotti cells . The AP half-width of our cells is also in line with room temperature recordings in juvenile rat slices from bitufted interneurons, some identifi ed as Martinotti cells (Ali, 2007), and from regular spiking interneurons (Cauli et al., 1997). AP durations reported for groups 2 and 3 are not unreasonable given the ranges in the literature for recordings in slice from juvenile animals at room temperature. The RMP and input resistance are in a healthy and normal range for all cells. Finally, as mentioned, group 1 and groups 2 and 3 differed statistically in their sag.
The fact that Martinotti cells comprise a subtype of SOM interneurons was already known, thus the separation of the Martinotti cells is an external validation of the results. Finally, the distinction between group 2 and 3 was preserved in both morphological and electrophysiological clustering, suggesting that the non-Martinotti cells are comprised of at least two subtypes. These two groups do not appear to be any of the SOM subtypes previously identifi ed in mouse cortex Halabisky et al., 2006;Ma et al., 2006), although they do have some similarities to neurons from the X94 somatostatin positive line, such as avoidance of layer 1 (groups 2 and 3) and stuttering properties (group 2). However, X94 cells have highly branched axons, contained mostly in layer 4. Thus, X94 cells and our group 2 and 3 cells target different layers. Additionally X94 cells are quasi-fast spiking, with much shorter ISI than group 2 or 3 cells. For the four electrophysiological subtypes identifi ed by Halabisky et al., two of the subtypes described have a similar level of spike frequency adaptation to our cells. In addition, like group 3 cells, all four subtypes they identifi ed have a low input resistance with mean input resistance for the four subtypes between 205 and 224 MΩ. However, in this study the sEPSC kinetics in addition to passive and active membrane properties were important in distinguishing between the four subtypes and morphologies of the cells were not recovered. Therefore, comparison of our results to these subtypes is diffi cult without sEPSC or morphological comparison. In futures studies it will be interesting to determine the molecular phenotype of these SOM-positive interneurons and identify any differentially expressed calcium binding proteins or neuropeptides.
The use of juvenile mice (P10-P18) raises the question of whether the subtypes identifi ed in this study are representative of cell types in adult mice. A comparison of the properties of layer 4 neocortical interneurons in slice from juvenile rats (P18-22) and adults rats (Ali et al., 2007) indicates that subtypes identifi ed in juveniles likely are representative of adult cell types. This study found that there was no signifi cant difference in axonal arbor size between the juvenile and adult. Correlations between fi ring properties and synaptic properties were also conserved from juvenile to adult. The only signifi cant difference identifi ed was that in adult animals, action potential, EPSP and IPSP half-widths were narrower and AP-AHPs were deeper than in juvenile animals. Given that the relationship between fi ring properties and synaptic properties is not dependent on age, fi ndings in juvenile animals are relevant. In particular the identity of different classes of interneurons described in the study remains the same between juveniles and adults. For example the bitufted interneurons had broad action potentials relative to other interneuron types in both juvenile and adult animals. As all cell types had slower events in juvenile tissue, the differences between groups remain unaffected by age.

DEVELOPMENTAL ORIGIN OF SOM INTERNEURONS
Given the increasingly sophisticated understanding of SOM interneuron development, it is imperative to discuss whether there is any evidence for different subtypes of SOM interneurons based on this knowledge. Of particular relevance is the clear relation between the neuronal cell type and spatio-temporal anlage of fate determination (Wonders and Anderson, 2006;Dasen and Jessell, 2009). SOM interneurons, along with the majority of neocortical interneurons, originate from the medial ganglion eminence (MGE) (Xu et al., 2004;Wonders and Anderson, 2006), although it has been so far diffi cult to pinpoint exactly the precise spatio-temporal profi le of their precursors. By using inducible genetic fate mapping, recent studies have addressed whether different interneuron types originate at different points in development. Based on molecular content and electrophysiology, distinct groups have been reported (Miyoshi et al., 2007;Sousa et al., 2009). By fate mapping Olig-2 expressing precursors (restricted to MGE where Olig-2 expression is high) at different times in development, it was found that in mouse somatosensory cortex barrel fi eld SOM-positive, calretin (CR) negative interneurons make up 30% of early populations (precursors labeled at E9.5-10.5), but their generation drops as development continues, comprising almost none of the late population (precursors labeled at E15.5) (Miyoshi et al., 2007). In an analogous study fate mapping Nkx2.6 precursors, a similar temporal pattern was observed for SOM-positive CR-negative interneurons. In addition they found that SOM-positive CR-positive interneurons make up the majority of the late population (precursors labeled after E12.5; Sousa et al., 2009). The SOM positive CR positive interneurons have multipolar dendrite morphology and regular spiking patterns, while the SOM positive CR negative interneurons are identifi ed as Martinotti cells. Therefore it is possible that the three groups identifi ed in this study originated at different developmental points, or perhaps in different spatial positions. We would therefore venture the hypothesis that the lack of consensus in identifying the origin of SOM interneurons refl ects their basic heterogeneity, so their classifi cation into three (or potentially more) subgroups could help better delineate the spatio-temporal origin of each subgroup.
The very distinct axon morphology of group 2 and 3 cells may be the result of a common developmental response to molecular cues, and is thus congruous with the groups originating at different times or from different locations in the MGE. The avoidance of layer 1 by group 2 and 3 cells, while group 1 cells have extensive branching in layer 1 could be due to either opposite responses to the layer 1 molecular cues or attraction of group 2 and 3 cells to axon guidance cues in layer 2/3. Such a feature can be seen on some of the layer 4 short axon cells drawn by Lorente de No (1922).

POSSIBLE FUNCTIONAL CONSEQUENCES
The distinct morphology and physiology of each of the three types implies that each could have a different functional role within the circuit. First, the differently shaped dendritic arborizations mean that the groups have different laminar inputs. In addition, differing axon morphologies imply that downstream targets of the groups are distinct. In particular, Group 2 and 3 cells' axons avoid layer 1 while Martinotti cells branch extensively in this layer. The modulation of dendritic excitability in layer 1 is important as connections from higher cortical areas are received in layer 1 (Coogan and Burkhalter, 1993). Through a connected network of pyramidal cells and Martinotti cells, excitation of a single pyramidal neuron can recruit Martinotti cells which in turn cause inhibition of the apical tuft dendrites of neighboring pyramidal cells (Kozloski et al., 2001;Kapfer et al., 2007;Silberberg and Markram, 2007). This pathway is extremely sensitive, with a supralinear increase in recruited Martinotti cells as the number of pyramidal cells excited increases (Kapfer et al., 2007).
The medial directional preference of group 2 and 3 axons suggests that this asymmetry could have a functional purpose. This asymmetry resembles other incidences of neurons with asymmetric processes with potential functions. For example, retinal OFF direction selective ganglion cells have asymmetric dendrites with a dorsal to ventral direction preference. This asymmetry translates to an asymmetric and selective receptive fi eld which allows these cells to detect only upwards movement (Kim et al., 2008). These retinal ganglion cells have a distinct molecular identity and "tile" the retina, such that each distinct ganglion cell type covers the retina with minimal overlap (De Vries and Baylor, 1997).
In the neocortex there are also examples of asymmetric dendrites and axons. In macaque monkeys, Meynert cells in layer 6 of primary visual cortex have asymmetric basal dendrites, extending up to 600 µm and aligned parallel to the ocular dominance column (Winfi eld 1983;Livingstone, 1998). These cells have a direction selective receptive fi eld based in part on the dendritic conduction delays due to asymmetry of the basal dendrites (Livingstone, 1998). In primary somatosensory barrel neurons in layer 4 have asymmetric dendrites and axons that are almost always directed toward the center of the barrel (Woolsey 1975;Egger et al., 2008). It is postulated that this allows the dendrites to receive maximal input from thalamocortical afferents, while the localization of the axon provides minicolumns of recurrent excitation within the barrel. Distinct minicolumns could each process a receptive fi eld for a different quality of whisker movement. A few barrel neurons have axons directed to the neighboring barrel, providing interbarrel connections and allowing for multiwhisker receptive fi elds (Egger et al., 2008). Such examples of structure so clearly relating to a particular function could hold true for group 2 and 3 cells whose axons appear to have a very selective target, given the sparsity and directional preference of the axon.

CLUSTER ANALYSIS AS A TOOL FOR CLASSIFICATION OF NEURONAL SUBTYPES
The checkered history of neuronal classifi cation based on qualitative descriptors has led to a multitude of classifi cations. Unsupervised classifi cation methods appear ideal to objectively identify the subtypes of neurons present in biological circuits using quantitative criteria obtained by blind algorithms. At the same time, it is possible that supervised methods, using available information in the dataset, could outperform clustering algorithms for certain neuronal classifi cation tasks (L. Guerra, L. M. McGarry, P. Larrañaga, V. Robles, C. Bielza and R. Yuste, in prep.), so further methodological research appears necessary.
As cluster analysis becomes a more widely used tool for defi ning cell types, it is important to note that it should used in such a way as to optimize its power of objective and quantitative classifi cation, but while also controlling for its weaknesses. One problem of cluster analysis is that it always returns a classifi cation of the input data. Thus it is important that the results of cluster analysis are thoroughly validated for statistical signifi cance and evidence of natural groups, rather than artifi cial or arbitrary divisions.
For example, the quality of the clusters should also be examined with silhouette analysis as well as agreement between classifi cations based on multiple different independent sets of parameters. Additionally the weaknesses of an algorithm being used can bias the results. Thus, comparison to an alternate clustering algorithm, preferably a complementary algorithm such as bottom up hierarchical clustering paired with top down K-means clustering, tests for robustness of the classifi cation. All of these approaches were used in this study to assess our cluster analysis results, allowing us to derive more convincing evidence for our conclusions. Finally, in this study with Wards' hierarchical clustering, we have encountered no particular advantage to use a multimodal cluster analysis (i.e., one that included physiological and morphological features) over using just one (for example, the morphology, Figure 3). While this could refl ect the structure of our data, it should be kept in mind as a warning note for future studies using Wards method of multimodal datasets.
Using cluster analysis raises a fundamental question: what should be considered a distinct cell type? As cluster analysis will generate continual division, it becomes necessary to defi ne an endpoint that is both meaningful and useful in addition to being statistically signifi cant. Part of the problem is the uncertainty about the number of neocortical interneuron types, with estimates of as many as 100 cell types (Lorente de No, 1922;Stevens, 1998). There is still much controversy with regards to defi nition of a cell type, but several recommendations made by recent reviews have begun to address the issue, proposing some objective criteria, such as tiling of surface (Mott and Dingledine, 2003;Masland, 2004;Migliore and Shepherd, 2005) or their transcription factor expression (Yuste, 2005). The eventual goal is to defi ne functionally distinct groups, most often identifi ed through a combination of characteristic morphology, electrophysiology or molecular profi le. The continuation of quantitative classifi cation of interneuron subtypes, coupled with new fi ndings on both the developmental origin and functional role of these subtypes will make this goal come closer to realization. Such a complete description of interneuron subtypes is necessary to achieve the fi nal goal of mapping out the cortical circuitry.