Quantitative classification of somatostatin-positive neocortical interneurons identifies three interneuron subtypes.

Deciphering the circuitry of the neocortex requires knowledge of its components, making a systematic classification of neocortical neurons necessary. GABAergic interneurons contribute most of the morphological, electrophysiological and molecular diversity of the cortex, yet interneuron subtypes are still not well defined. To quantitatively identify classes of interneurons, 59 GFP-positive interneurons from a somatostatin-positive mouse line were characterized by whole-cell recordings and anatomical reconstructions. For each neuron, we measured a series of physiological and morphological variables and analyzed these data using unsupervised classification methods. PCA and cluster analysis of morphological variables revealed three groups of cells: one comprised of Martinotti cells, and two other groups of interneurons with short asymmetric axons targeting layers 2/3 and bending medially. PCA and cluster analysis of electrophysiological variables also revealed the existence of these three groups of neurons, particularly with respect to action potential time course. These different morphological and electrophysiological characteristics could make each of these three interneuron subtypes particularly suited for a different function within the cortical circuit.

Introduction

Despite comprising a minority of all neocortical neurons, GABAergic interneurons appear to play an important circuit role by responding to dynamic changes in excitation, perhaps to keep the circuit responsive over a wide range of inputs, synchronize activity, control timing of pyramidal cell firing, or suppress runaway excitation (McBain and Fisahn, 2001; Pouille and Scanziani, 2001; Markram et al., 2004; Trevelyan et al., 2006; Klausberger and Somogyi, 2008). In support of these roles, GABAergic neurons have been implicated in several pathologies including epilepsy (Cossart et al., 2001; Cobos et al., 2005; Freund and Katona, 2007; Trevelyan et al., 2006, 2007), autism (Tabuchi et al., 2007), Rett syndrome (Dani et al., 2005), anxiety disorder (Powell et al., 2003; Freund and Katona, 2007), Tourette syndrome (Kalanithi et al., 2005) and schizophrenia (Lewis et al., 2005; Gonzalez-Burgos and Lewis, 2008).

Interneurons exhibit very diverse morphology, electrophysiology, molecular content and post synaptic targets (Cauli et al., 1997; Markram et al., 2004; Yuste, 2005). Classification of neocortical interneurons is a crucial step in understanding cortical circuits as each subtype of interneuron likely has a different function. In the past, classifications have been based on one or a combination of descriptors, often using qualitative criteria which were not standardized. Therefore, interneuron classifications 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 first 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 classification of interneurons. In addition, to facilitate the classification of interneurons, a standardized nomenclature 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 neuropeptide 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 subtypes 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-filled 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 sub-groups of SOM-positive neurons that, to our knowledge, have not yet been described. The electrophysiological classification agreed well with, and thus confirmed, the morphological classification. Furthermore, significant differences in axon morphology and firing properties between the three groups suggest that the two novel subtypes and Martinotti cells could serve different functional roles in the cortical circuit.

Materials and Methods

Preparation of brain slices

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 NaHCO3, 1.5 mM NaH2PO4, 2 mM CaCl2, 2 mM MgSO4, bubbled with 95% 02, 5%CO2). 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 MgSO4, 2 mM CaCl2, 1 mM NaH2PO4, 26 mM NaHCO3 and 10 mM glucose, bubbled with 95% 02, 5%CO2).

Mice strain

GIN mice reliably label a subset of SOM-positive cells with expression remaining consistent for more than five 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 co-localization 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 flowing oxygenated ACSF. Pipettes of 3–7 MΩ resistance were pulled from borosilicate glass. SOM positive cells were identified 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.

Table 1

Electrophysiological variables. Action potential properties measured from response to twice threshold, 500-ms current injection from first 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
Input resistance (MΩ)	Calculated from small hyperpolarizing or depolarizing current steps (≤10-mV deflection)
Rheobase (pA)	Threshold current
AP1 amplitude(mV)	Amplitude of the 1st action potential (AP)
AP1 duration (ms)	Time from onset of 1st AP, calculated as an increase ≥1 mV/100 ms, to offset, calculated as return to same voltage as before AP onset
AP1 half-width (ms)	Time from half-amplitude during rise to half-amplitude during fall of 1st AP
AP1 rise time(ms)	Time from onset to peak of 1st AP
AP1 fall time (ma)	Time from peak to offset of 1st AP
AP1 rise rate (mV/ms)	AP1 amplitude/AP1 rise time
AP1 fall rate (mV/ms)	AP1 amplitude/AP1 fall rime
AP drop (mV)	AP1 amplitude-AP2 amplitude
Spike frequency adaptation	t2/t1, where t1 is the time between the peaks of the first two APs (interspike interval) and t2 is the last interspike interval

Histological procedure

Neurons were filled 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 ddH20 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 final 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 filled 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 defined 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 2

Morphological variables. Variables were extracted using the Neurolucida Explorer program by MicroBrightField.

Variable	Description
VARIABLES DESCRIBING THE SOMA
Somatic perimeter (μm)	Perimeter of the soma
Somatic area (μm2)	Area of the soma
Somatic aspect ratio	Max diameter of soma/min diameter of soma
Somatic compactness	((4/π)*Area)1/2/max diameter
Somatic form factor	(4π*Area)/(Perimeter2)
Somatic roundness	(4*Area)/(π* max diameter2)
VARIABLES DESCRIBING AXON
Axonal node total	Total number of axonal nodes (branching points)
Total axonal length (μm)	Sum of lengths of all axon segments, measured along tracing (not straight line distance)
Total surface area of axon (μm2)	2πr2 + 2πrh, SA calculated by modeling axon as a cylinder with diameter defined by thickness of segment in reconstruction
Ratio of axonal length to surface area (1/μm)	Total axonal length/total surface area of axon
Highest order axon segment	Maximum number obtained after each segment is numbered by how many nodes it is removed from the initial segment
Axonal torsion ratio	Total axonal length/total axonal length of fan in diagram where the fan in diagram is 2-D projection of the neuron constructed by compiling traces swept around a vertical axis. Torsion ratio = 1 corresponds to no loss of length, values larger than 1 correspond to the factor by which the processes have decreased in the fan in diagram
K-dim of axon	Fractal dimension of the axon calculated using linear regression and the nested cubes method
Axonal polar angle average	Average of polar angles of all axonal nodes. The polar angle is the angle between the 2 lines passing through the node and the endpoints of the next segments.
Axonal polar angle standard deviation	Standard deviation of axonal polar angles
Axonal local angle average	Average of local angles of all axonal nodes. The local angle is the angle between the 2 lines passing through the node and points adjacent to the node on the two following segments
Axonal local angle standard deviation	Standard deviation of axonal local angles
Axonal spline angle average	Average of spline angles of all axonal nodes. The spline angle is the angle between the 2 lines passing through the node and smoothed points adjacent to the node when the following two segments are approximated by a cubic spline
Axonal spline angle standard deviation	Standard deviation of axonal spline angles
Average tortuosity of axonal segments	Average of tortuosities measured for each axonal segment. Segment tortuosity = distance along segment/straight line distance between segment endpoints
Standard deviation of tortuosity of axonal segments	Standard deviation of tortuosities of all axonal segments
Axonal segment length average (μm)	Total axonal length/number of segments
Axonal segment length standard deviation (μm)	Standard deviation of axonal segment length
Average tortuosity of axonal nodes	Average of tortuosities measured for each axonal node. Node tortuosity = distance along process from origin of process to node/straight line distance from origin of process to node
Standard deviation of tortuosity of axonal nodes	Standard deviation of tortuosities of all axonal nodes
Number of axonal sholl sections	Number of sholl sections (concentric spheres centered at the soma with radii at 100 μm intervals) containing axonal processes
Axonal sholl length at 100 μm	Total length of axonal segments contained in first sholl section/total axonal length
Axonal sholl length at 200 μm	Total length of axonal segments contained in second sholl section/total axonal length
Axonal sholl length at 300 μm	Total length of axonal segments contained in third sholl section/total axonal length
Axonal sholl length density (μm)	Total axonal length/number of axonal sholl sections
Axonal sholl node density	Axonal node total/number of axonal sholl sections
Convex hull axon area (μm2)	Area of the 2-D convex polygon created by connecting the distal axon segment endpoints of 2-D projection of neuron
Convex hull axon perimeter (μm)	Perimeter of the 2-D convex polygon created by connecting the distal axon segment endpoints of 2-D projection of neuron
Convex hull axon volume (μm3)	Volume of the 3-D convex polygon created by connecting the distal axon segment endpoints
Convex hull axon surface area (μm2)	Surface area of the 3-D convex polygon created by connecting the distal axon segment endpoints
Axon node density (1/μm)	Axonal node total/total axonal length
VARIABLES DESCRIBING DENDRITES
Number of dendrites	Total number of dendrites
Dendritic node total	Total number of dendritic nodes (branching points)
Total dendritic length (μm)	Sum of lengths of all dendrite segments, measured along tracing (not straight line distance)
Average length of dendrites (μm)	Total dendritic length/number of dendrites
Total surface area of dendrites (μm2)	See total surface area of dendrites
Ratio of dendritic length to surface area (1/μm)	See ratio of axonal length to surface area
Highest order dendrite segment	See highest order axonal segment
Dendritic torsion ratio	See axonal torsion ratio
K-dim dendrites	See K-dim axon
Dendritic polar angle average	See axonal polar angle average
Dendritic polar angle standard deviation	See axonal polar angle standard deviation
Dendritic local angle average	See axonal local angle average
Dendritic local angle standard deviation	See axonal local angle standard deviation
Dendritic spline angle average	See axonal spline angle average
Dendritic spline angle standard deviation	See axonal spline angle standard deviation
Average tortuosity of dendritic segments	See average tortuosity of axonal segments
Standard deviation of totuosity of dendritic segments	See standard deviation of totuosity of axonal segments
Dendritic segment length average (μm)	See axonal segment length average
Dendritic segment length standard deviation (μm)	See axonal segment length standard deviation
Average tortuosity of dendritic nodes	See average tortuosity of axonal nodes
Standard deviation of tortuosity of dendritic nodes	See standard deviation of tortuosity of axonal nodes
Number dendritic sholl sections	Number of sholl sections (concentric spheres centered at the soma with radii at 50-μm intervals) containing dendritic processes
Dendritic sholl length at 50-μm	Total length of dendritic segments contained in first sholl section/total dendritic length
Dendritic sholl length at 100-μm	Total length of dendritic segments contained in second sholl section/total dendritic length
Dendritic sholl length at 150-μm	Total length of dendritic segments contained in third sholl section/total dendritic length
Convex hull dendrite area (μm2)	See convex hull axon area
Convex hull dendrite perimeter (μm)	See convex hull axon perimeter
Convex hull dendrite volume (μm3)	See convex hull axon volume
Convex hull dendrite surface area (μm2)	See convex hull axon surface area
Dendrite node density (1/μm)	See axon node density
VARIABLES DESCRIBING LOCATION
Relative distance	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 defined by the old variables to a new lower dimensional coordinate system defined 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 classified 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 find the clustering levels which have significance 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, defined as the standardized difference between a fusion value and the mean fusion value, are compared to find significantly large realized deviates. A large realized deviate is significant 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 find the significant 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. Significant departures of the actual fusion values from random, determined by the t-statistic, indicate a number of clusters that are statistically significant. The first level at which significance 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).

Figure 1

Anatomical classification of SOM neurons. (A) Ward's method of hierarchical unsupervised clustering based on 67 morphological variables applied to 39 SOM positive interneurons. The first 2 principal components were retained for cluster analysis. The x-axis of dendrogram shows individual cells and the y-axis represents the linkage distance measured by the Euclidean distance squared. Open circles indicate the centroid cell of each cluster. Black brackets outline the clusters statistically significant at the 5% level and a horizontal line indicates this cut-off linkage distance. The groups discussed in the text are colored: group 1 in purple, group 2 in orange and group 3 in green. (B) Scatterplot of dataset in principal component space. The x-axis represents the first principal component (PC1) and the y-axis represents the second principal component (PC2). Centroids of clusters are labeled and circled in red. Orthogonal lines separating the three groups are shown. (C) Neurolucida reconstructions of representative cells of each cluster. Axons are shown in blue and dendrites in red.

Figure 2

Physiological classification of SOM neurons. (A) Ward's method of hierarchical unsupervised clustering based on 19 electrophysiological variables applied to 36 SOM positive interneurons. The first 2 principal components were retained for cluster analysis. The x-axis of dendrogram shows individual cells and the y-axis represents the linkage distance measured by Euclidean distance squared. Open circles indicate the centroid cell of each cluster. Black brackets outline the clusters statistically significant at the 5% level and a horizontal line indicates this cut-off linkage distance. Group 1 (Martinotti cells) are shown in purple, group 2 in orange and group 3 in green. This color scheme is based on the clustering by morphological variables (See Figure 1) and will be preserved for all figures. Black cells are those without morphological reconstruction. (B) Scatterplot of dataset in principal component space. The x-axis represents the first principal component (PC1) and the y-axis represents the second principal component. Centroids of clusters are labeled and circled in red. Orthogonal lines separating the three groups are shown. The outlier in cluster e is indicated with an arrow. (PC2). (C) Current clamp recording of response to twice threshold current pulse for the centroid cells of each cluster.

Figure 3

Silhouette analysis. (A) Plot of silhouette values for clustering of 39 cells by morphological variables (see Figure 1A for dendrogram of this dataset). On the y-axis cells in each cluster are ordered by decreasing silhouette value. The silhouette value can range between −1 and 1. Large positive values indicate clusters are distinct with greater intra-cluster similarity than between cluster similarity (See Materials and Methods for further explanation). The x-axis represents the silhouette value (See Materials and Methods). (B) Plot of silhouette values for clustering of 36 cells by electrophysiolgical variables (see Figure 2A for dendrogram of this dataset). (C–E) Plot of silhouette values (left) and dendrogram (right) for clustering of the 16 cells with both recordings and reconstructions by morphological variables (C), by electrophysiological variables (D) and by electrophysiological and morphological variables (E). The first 2 principal components were retained for cluster analysis by the morphology variables (C) and electrophysiological variables (D). The first 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).

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 significant 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

s(i) = (b(i)−a(i))/max[a(i),b(i)]

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 defined 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 identified by cluster analysis, the group mean and standard error was calculated for every variable and the groups were compared for significant 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 significant differences between the groups were highly correlated with the PCs (Tables 6 and 7)

Table 6

Correlations between morphological variables and the first 2 principal components of the 39-cells morphology variables dataset.

	PC 1	PC 2
Somatic perimeter	0.06917	−0.244601
Somatic area	0.37328	0.274426
Somatic aspect ratio	0.02352	−0.294285
Somatic compactness	0.06039	0.242742
Somatic form factor	0.01698	0.401195
Somatic roundness	0.05458	0.204208
Axonal node total	−0.92976	0.050520
Total axonal length	−0.93096	−0.022486
Total surface area of axon	−0.94186	0.018395
Ratio of axonal length to surface area	0.36009	0.044226
Highest order axon segment	−0.90042	−0.095688
Axonal torsion ratio	−0.67806	0.090679
k-dim axon	−0.84822	−0.057876
Axonal planar angle avg	0.26554	0.108213
Axonal planar angle stdv	0.15675	0.061679
Axonal local angle avg	0.62546	−0.125786
Axonal local angle stdv	0.25285	−0.371129
Axonal spline angle avg	0.64852	0.038303
Axonal spline angle stdv	0.28264	−0.050636
Avg tortuosity of axonal segments	−0.08317	−0.384847
Stdv of tortuosity of axonal segments	−0.45397	−0.338419
Axonal segment length avg	0.58497	0.180537
Axonal segment length stdv	0.53072	0.210869
Avg tortuosity of axonal nodes	−0.66627	−0.170654
Stdv of tortuosity of axonal nodes	−0.71848	−0.157837
Number axonal sholl sections	−0.60628	−0.026166
Axonal sholl length at 100 μm	0.27203	−0.050482
Axonal sholl length at 200 μm	0.04591	0.166692
Axonal sholl length at 300 μm	0.21472	0.270899
Axonal sholl length density	−0.90404	−0.037149
Axonal sholl node density	−0.89645	0.024540
Convex hull axon area	−0.83877	0.087294
Convex hull axon perimeter	−0.82856	0.042246
Convex hull axon volume	−0.68435	0.264920
Convex hull axon surface area	−0.86393	0.125719
Axon node density	−0.75277	−0.043953
Number of dendrites	−0.02817	0.287928
Dendritic node total	−0.05756	0.886395
Total dendritic length	0.12102	0.811526
Avg length of dendrites	0.23202	0.629327
Total surface area of dendrites	−0.29875	0.743853
Ratio of dendritic length to surface area	0.23126	0.039592
Highest order dendritic segment	0.03966	0.707728
Dendritic torsion ratio	−0.31617	0.056369
k-dim dendrites	0.16448	0.545213
Dendritic planar angle avg	−0.10496	−0.162805
Dendritic planar angle stdv	0.13594	−0.176268
Dendritic local angle avg	0.52598	−0.471847
Dendritic local angle stdv	0.57963	−0.379845
Dendritic spline angle avg	0.48650	−0.341205
Dendritic spline angle stdv	0.36733	0.046085
Avg tortuosity of dendritic segments	0.01934	−0.375656
Stdv of tortuosity of dendritic segments	0.03568	−0.090445
Dendritic segment length avg	0.06082	−0.577189
Dendritic segment length stdv	0.13379	−0.267916
Avg tortuosity of dendritic nodes	−0.12717	−0.241970
Stdv of tortuosity of dendritic nodes	−0.09880	−0.147266
Number dendritic sholl sections	0.14490	0.136173
Dendritic sholl length at 50μm (fraction)	−0.16094	−0.106816
Dendritic sholl length at 100μm (fraction)	−0.03877	−0.064122
Dendritic sholl length at 150μm (fraction)	−0.09192	0.450425
Convex hull dendrite area	0.23731	0.784130
Convex hull dendrite perimeter	0.27421	0.651568
Convex hull dendrite volume	−0.26638	0.770294
Convex hull dendrite surface area	0.18344	0.819703
Dendrite node density	−0.09162	0.684138
Relative distance to pia	−0.46770	0.075928

Correlations > 0.7 are highlighted in red. These 2 principal components were retained for cluster analysis.

Table 7

Correlations between electrophysiological variables and the first 2 principal components of the 36-cells electrophysiology variables dataset.

	PC1	PC2
Resting membrane potential (mV)	0.279283	−0.713146
AP1 Amplitude (mV)	−0.619046	0.700534
AP1 duration(ms)	0.966063	0.149507
AP1 half-width (ms)	0.989002	0.055855
AP1 rise time(ms)	0.934992	0.169495
AP1 fall time(ms)	0.949161	0.142030
AP2 Amplitude(mV)	−0.620920	0.654053
AP2 duration(ms)	0.929245	0.324156
AP2 half-width(ms)	0.942784	0.237551
AP2 rise time(ms)	0.870369	0.269250
AP2 fall time(ms)	0.908103	0.323843
AP Drop(mV)	−0.348181	0.606961
Spike frequency adaptation	−0.412130	0.517381
Input resistance (MΩ)	0.387743	−0.214216

Correlations > 0.7 are highlighted in red. These 2 principal components were retained for cluster analysis.

Results

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, filled 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 the patched cells. From a dataset of 59 SOM-positive interneurons from mouse primary somatosensory, visual and frontal cortex, 39 cells had complete reconstructions, 36 cells had high quality electrophysiological recordings, and 16 cells had both complete reconstructions and recordings.

We first 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 five 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 first 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 first 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 first 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 significant 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).

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 filling 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 defined 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 significant 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 first two clusters (a and b) constituted a group (group 1, purple), in which, remarkably, all of its reconstructed neurons had been previously identified as belonging to the group 1 from the morphological clustering (Figure 1A, purple cells). Besides this first group, the remaining 3 clusters (c–e) contain all neurons previously identified as belonging to the morphological groups 2 and 3 and, moreover, separated these cells in good agreement with these 2 groups. More specifically, 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 first 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 firing pattern, non-stuttering 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 definition 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 significantly more depolarized resting membrane potential (Table 3). On the other hand, group 2 cells (n = 7) displayed both regular spiking and stuttering firing 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 significantly 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 firing 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 first few spikes and the other showed amplitude accommodation similar to group 2 cells (Figure S2).

Table 3

Average values ± standard error for the morphological and electrophysiological variables with significant differences between groups.

Parameter	1	2	3	*	**
Axonal node total	258.52 ± 28.16	14.55 ± 2.31	12.71 ± 6.10		1,2; 1,3
Total axonal length (μm)	15561.014 ± 1427.684	1441.91 ± 181.77	1638.16 ± 5674.09		1,2; 1,3
Total surface area of axon (μm2)	11732.37 ± 1162.32	761.23 ± 119.50	751.92 ± 278.47		1,2; 1,3
Highest order axon segment	28.81 ± 1.58	7.82 ± 0.78	7.00 ± 1.59		1,2; 1,3
K-dim of axon	1.308 ± 0.019	1.118 ± 0.030	1.072 ± 0.033		1,2; 1,3
Standard deviation of tortuosity of axonal nodes	1.78180 ± 0.20810	0.33380 ± 0.07713	0.59388 ± 0.29117		1,2; 1,3
Axonal sholl length density	2003.05 ± 206.73	389.90 ± 60.51	271.20 ± 57.43		1,2; 1,3
Axonal sholl node density	33.953 ± 4.051	4.497 ± 1.081	1.947 ± 0.712		1,2; 1,3
Axon convex hull area (μm2)	382893.18 ± 36652.35	54929.40 ± 12429.32	55076.94 ± 133820.97		1,2; 1,3
Axon convex hull perimeter (μm)	2489.32 ± 148.76	973.75 ± 130.80	1031.27 ± 152.76		1,2; 1,3
Axon convex hull surface area (μm2)	835531.38 ± 81482.02	112537.67 ± 25745.01	113458.51 ± 27038.39		1,2; 1,3
Axon node density (μm−1)	0.017 ± 0.001	0.009 ± 0.001	0.006 ± 0.001		1,2; 1,3
Dendritic node total	23.95 ± 4.30	20.73 ± 0.94	30.29 ± 3.037		2,3
Total dendritic length (μm)	2243.148 ± 226.407	2241.845 ± 144.953	3097.971 ± 319.028		2,3
Dendrite convex hull area (μm2)	55682.88 ± 7889.45	55123.70 ± 7962.79	101859.29 ± 10149.84		1,3; 2,3
Dendrite convex hull volume (μm3)	1286065.99 ± 345715.10	641800.64 ± 89721.61	1208996.01 ± 11724.36		2,3
Dendrite convex hull surface area (μm2)	119939.05 ± 16552.37	113281.71 ± 15994.98	207289.43 ± 20315.47		1,3; 2,3
Resting membrane potential (mV)	−65.336 ± 0.438	−70.617 ± 2.385	−70.236 ± 2.228		1,2; 1,3
AP1 amplitude (mV)	55.407 ± 2.548	61.602 ± 7.865	70.245 ± 2.425	1,2	1,3
AP1 duration (ms)	3.924 ± 0.198	5.914 ± 0.734	3.000 ± 0.202	1,3	1,2
AP1 half-width (ms)	1.649 ± 0.083	2.271 ± 0.291	1.356 ± 0.071		2,3
AP1 rise time (ms)	1.070 ± 0.036	1.400 ± 0.102	0.880 ± 0.037	1,2	1,3
AP1 fall time (ms)	2.853 ± 0.168	4.514 ± 0.656	2.120 ± 0.166		1,2; 2,3
AP2 amplitude (mV)	52.788 ± 2.367	56.781 ± 6.782	64.014 ± 2.458		1,3
AP2 duration (ms)	4.282 ± 0.212	6.829 ± 0.762	3.580 ± 0.263		1,2; 2,3
AP2 half-width (ms)	1.829 ± 0.096	2.729 ± 0.321	1.560 ± 0.125		1,2; 2,3
AP2 rise time (ms)	1.170 ± 0.044	1.586 ± 0.116	1.040 ± 0.051		1,2; 2,3
AP2 fall time (ms)	3.112 ± 0.176	5.243 ± 0.699	0.2540 ± 0.220		1,2; 2,3
AP drop (mV)	2.619 ± 0.479	4.822 ± 1.424	6.230 ± 0.911		1,2; 1,3
Input resistance (MΩ)	469.583 ± 35.836	475.429 ± 57.726	256.400 ± 46.427		1,3; 2,3
AP1 rise rate (mV/ms)	54.265 ± 3.733	48.079 ± 9.906	80.337 ± 4.179	2,3	1,3
AP1 fall rate (mV/ms)	21.812 ± 1.893	17.304 ± 4.749	33.921 ± 2.836	2,3	1,3
AP2 rise rate (mV/ms)	47.308 ± 3.178	38.509 ± 7.135	62.224 ± 4.383	2,3	1,3
AP2 fall rate (mV/ms)	18.807 ± 1.584	12.850 ± 2.944	26.062 ± 2.692	1,3	2,3
Sag vs. Vm slope	0.201 ± 0.022	0.065 ± 0.050	0.024 ± 0.012		1,2; 1,3

Significance level marked in the last two columns, *p < 0.05, **p < 0.01. p-values calculated by the Mann Whitney U Test.

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. Specifically, 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), first apparent from the morphological clustering.

These results display a remarkably good agreement between the morphological and physiological classifications, which because they are based on completely independent datasets, confirm each other. In fact, with this basic classification into three groups every neuron is classified 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 find 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 classification based on cluster analysis actually was, using a variety of statistical analyses. We first 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. Specifically, the best cut test distinguishes the clustering profile of a dataset containing groups from the profile 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 finds if there is a linkage cut-off distance at which this fusion value profile occurs. The upper tail test also locates significantly large fusion values. The bootstrap test determines the linkage distance at which the results significantly deviate from random. Indeed, using these three tests, the clusters (as marked by brackets in Figures 1A and 2A) were found to be statistically significant at 5% level (See Materials and Methods for details). Thus, three independent significance tests all show that the dataset does have a non-uniform structure, and thus could be classified into meaningful clusters.

Following this, we sought to confirm whether the hierarchical classification using Ward's method was reliable, using K-means clustering with K specified as the number of significant 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 classification is inflexible 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.

Table 4

Agreement of hierarchical clustering (HC) using Ward's method with K-means clustering for K = number of clusters determined significant by HC.

Clustering Dataset	Number of cells in different cluster	Agreement between HC and K-means
39 cells, morphology	4	89.74%
36 cells, electrophysiology	1	97.22%
16 cells, morphology and electrophysiology	2	87.5%
16 cells, electrophysiology	1	93.75%
16 cells, morphology	0	100%

Sixteen cells refer to the cells with both complete reconstructions and recordings.

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 artificially imposed groups. Using this analysis, we found that the silhouette widths of the clustering by morphology and clustering by electrophysiology confirmed 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 classified.

As a final additional measure of statistical significance, 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.

Table 5

Comparison of silhouette widths of each dataset to the mean silhouette width of 500 trials of clustering randomized data.

Dataset	S(actual)	S(random) ± SD
39 morphology	0.547213	−0.23095 ± 0.180533
36 electrophysiology	0.490572	−0.18934 ± 0.097226
16 morphology and electrophysiology	0.276167	−0.08121 ± 0.173016
16 morphology	0.668405	−0.08059 ± 0.186636
16 electrophysiology	0.463107	0.004352 ± 0.179617

Quality of clustering combining different groups of variables

To determine the quality of classification 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 non-overlapping 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 identified. 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 classification 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 first 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 significant 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 (Wang et al., 2004; 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).

Figure 4

Axonal features of group 2 and 3 cells. (A) Axons of group 2 and 3 cells avoid layer 1. The 9 cells with axons that turn or hook are shown, all oriented with respect to the medial-lateral axis above. Note the preference of direction medially. (B) Light microscope images of section of two cells’ axons in the area boxed in the reconstructions below. Swellings in the axon indicated by arrow heads are possible boutons, found along the axons in layer 2/3.

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 significant 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 firing pattern (Table 2, Figure 5). There were significant differences among the passive properties of the three groups (Figure 5A). Group 2 and 3 cells were significantly (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 significantly (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 fitted to the sag vs. membrane potential was used as the variable for comparison, with sag defined 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 significantly 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.

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 first 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)

We found that the first 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). Specifically, 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 half-width 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 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 significantly 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.

In addition, all three groups showed spike frequency adaptation, with group 2 and 3 cells showing a greater degree of adaptation (spike frequency adaptation values group 1: 1.674 ± 0.092, group 2: 2.244 ± 0.385, group 3: 2.657 ± 0.493). Group 2 and 3 cells had greater action potential amplitude as measured from the first two action potentials in a train elicited by twice threshold current injection (AP1 amplitude p < 0.05, 61.602 ± 7.865 mV and p < 0.025, 70.245 ± 2.425 mV respectively) than group 1 (AP1 amplitude 55.407 ± 2.548 mV) and also had a significantly larger (p < 0.025, p < 5 × 10−3) drop in amplitude between these first two action potentials (AP drop 4.822 ± 1.424 mV and 6.230 ± 0.911 mV respectively) than Group 1 (AP drop 2.619 ± 0.479 mV; Figure 4A, Table 3).

Most cells were regular spiking, but three group 2 cells (75, 76, 77) were stuttering, defined as the firing of groups of action potentials that are interrupted by silent periods of varying length (Markram et al., 2004; 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 firing pattern of some group 2 cells. It can additionally be noted that group 2 cells had greater variability in firing 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).

Discussion

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 find a statistically consistent structure which can be captured well by a simple model of three basic subtypes of cells, already apparent in the first 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. Specifically, 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 filling or sectioned processes. However, as discussed, we did not find 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 filling, such as dye blown out of the cells or tapering and spotty filling at the end of visible processes, it is still possible that these cells were incompletely filled 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 difficult to visualize in its entirety. In any case, given that our data was collected from brain slices, we are the first 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.

Specifically, 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 difficult 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 filling. 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 filling 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 significantly different across the three groups, because they reflect 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 significant 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 filled. 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 significantly 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 significantly 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 significantly 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 filling, distinct from the sparse axon branching of groups 2 and 3 (again, not significantly 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, filled 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 filling.

Further evidence that group 2 and 3 are real comes from their electrophysiological properties. Specifically, 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 filled 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 (Wang et al., 2004). The AP half-width of our cells is also in line with room temperature recordings in juvenile rat slices from bitufted interneurons, some identified 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 identified in mouse cortex (Wang et al., 2004; 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 identified 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 identified 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 difficult 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 identified 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 identified in juveniles likely are representative of adult cell types. This study found that there was no significant difference in axonal arbor size between the juvenile and adult. Correlations between firing properties and synaptic properties were also conserved from juvenile to adult. The only significant difference identified 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 firing properties and synaptic properties is not dependent on age, findings 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 difficult to pinpoint exactly the precise spatio-temporal profile 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 field 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 identified as Martinotti cells. Therefore it is possible that the three groups identified 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 reflects their basic heterogeneity, so their classification 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 field 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 (Winfield 1983; Livingstone, 1998). These cells have a direction selective receptive field 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 field 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 fields (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 classification based on qualitative descriptors has led to a multitude of classifications. Unsupervised classification 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 classification 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 defining cell types, it is important to note that it should used in such a way as to optimize its power of objective and quantitative classification, but while also controlling for its weaknesses. One problem of cluster analysis is that it always returns a classification of the input data. Thus it is important that the results of cluster analysis are thoroughly validated for statistical significance and evidence of natural groups, rather than artificial or arbitrary divisions. For example, the quality of the clusters should also be examined with silhouette analysis as well as agreement between classifications 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 classification. 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 reflect 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 define an endpoint that is both meaningful and useful in addition to being statistically significant. 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 definition 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 define functionally distinct groups, most often identified through a combination of characteristic morphology, electrophysiology or molecular profile. The continuation of quantitative classification of interneuron subtypes, coupled with new findings 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 final goal of mapping out the cortical circuitry.

Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Supplementary Material

The Supplementary Material for this article can be found online at  http://www.frontiersin.org/neuralcircuits/paper/10.3389/fncir.2010.00012/