PMCCPMCCPMCC

Search tips
Search criteria 

Advanced

 
Logo of nihpaAbout Author manuscriptsSubmit a manuscriptHHS Public Access; Author Manuscript; Accepted for publication in peer reviewed journal;
 
Neural Comput. Author manuscript; available in PMC 2010 July 8.
Published in final edited form as:
PMCID: PMC2900194
NIHMSID: NIHMS214086

A Classification Method to Distinguish Cell-Specific Responses Elicited by Current Pulses in Hippocampal CA1 Pyramidal Cells

Abstract

The suprathreshold electrophysiological responses of pyramidal cells have been grouped into large classes such as bursting and spiking. However, it is not known whether, within a class, response variability ranges uniformly across all cells or whether each cell has a unique and consistent profile that can be differentiated. A major difficulty when comparing suprathreshold responses is that slight variations in spike timing in otherwise very similar traces render traditional metrics ineffective.

To address these issues, we developed a novel distance measure based on fiducial points to quantify the similarity among traces with trains of action potentials and applied it together with classification techniques to a set of in vitro patch clamp recordings from CA1 pyramidal cells. We tested if responses to repeated current stimulation of a given cell would cluster together yet remain distinct from those of other cells.

We found that depolarizing and hyperpolarizing current pulses elicited responses in each cell that clustered and were systematically distinguishable from responses in other cells. The fiducial-point distance measure was more effective than other methods based on spike times and voltage-gradient phase planes. Depolarizing traces were more reliably differentiated than hyperpolarizing traces, and combining both scores was even more effective.

These results suggest that each CA1 pyramidal cell has unique properties that can be detected and quantified with methods discussed here. This uniqueness may be due to slight variations in morphology or membrane channel densities and kinetics, or to large, coordinated variations in these elements. Ascertaining the actual sources and their degree of variability is important when constructing network models of neural function to ensure that key mechanisms are robust in the face of variations within these ranges. The analytical tools presented here can assist in constructing detailed cell models to match experimental records to elucidate the sources of electrophysiological variability in neurons.

1 Introduction

Repetitive firing in pyramidal neurons has been described in numerous studies attempting to elucidate the ionic currents that regulate its initiation, timing, termination, repolarization, afterhyperpolarization, and other qualities (Madison & Nicoll, 1984; Storm, 1987, 1990; Sah, Gibb, & Gage, 1988). Over time a picture has emerged of a complex interplay among a multitude of channels having different activation and deactivation voltage and time dependencies, with densities that vary over the cell’s morphological axis (Magee & Johnston, 1995; Hoffman, Magee, Colbert, & Johnston, 1997; Magee, 1998).

The complexity of the processes involved raises the question of whether successive responses are similar or whether elaborate interactions cause dissimilar responses, perhaps via chaotic-like processes. And differences among neurons in terms of cell-specific factors such as morphology and membrane channels raise the question of how different the responses are across cells. Visual inspection of published reports suggests that trains of action potentials may have consistent features when elicited repeatedly (see, e.g., Figures 1 and 4 of Madison & Nicoll, 1984). At the same time, although traces from different cells exhibit some general common features, trace dissimilarities may be sufficient to systematically differentiate among cells.

To date, there are no reports that systematically quantify the similarities or dissimilarities among traces with trains of action potentials, though spike timing has received some attention (Mainen & Sejnowski, 1995). Waveforms are frequently compared using the sum of squared differences over the length of their traces. However, action potentials in a train are narrow with slight variations in timing that make direct application of this measure problematic, perhaps explaining the paucity of studies quantifying their similarities. Techniques to compare traces containing action potentials have received more attention in modeling studies (Vanier & Bower, 1999; Keren, Peled, Korngreen, 2005; Achard & De Schutter, 2006; Weaver & Wearne, 2006), though their application has mostly remained confined to simulations where they are typically used as match or objective functions coupled to automated parameter search methods.

Between-cell variability has received attention recently with the realization that each neuron is constantly rebuilding itself from its constituent proteins, replacing ion channel proteins that turn over in the membrane with half-lives of minutes, hours or longer, and by consideration of how this process may be regulated (Marder & Goaillard, 2006). Results from experimental and modeling studies have begun to challenge beliefs that held that similar cell behaviors imply similarities in constituent elements. Experimental estimates of channel densities in cells that have similar electrophysiological profiles have shown substantially different ratios of inward Na+ and Ca2+ currents in Purkinje cells (Swensen & Bean, 2005) or among various K+ currents in stomatogastric ganglion neurons (Schulz, Goaillard, & Marder, 2006). Modeling studies have produced analogous results (Prinz, Bucher, & Marder, 2004; Achard & De Schutter, 2006).

In this study, we asked whether hippocampal CA1 pyramidal neurons exhibit cell-specific profiles in their electrophysiological responses. We collected a number of responses from several neurons to test if, for a given cell, responses group together yet remain distinct from those of other cells. To do this, we developed a novel method based on fiducial points to measure the similarity of traces with trains of action potentials and compared its performance to other published methods. Some of the results presented here have appeared in abstract form (Ambros-Ingerson, Grover, & Holmes, 2005).

2 Methods

2.1 Experimental Procedures

The data set of in vitro patch clamp recording used in the study was obtained as described below.

2.1.1 Slice Preparation and Maintenance

Sprague-Dawley rats (30–60 days old) were sedated by inhalation of a CO2/air mixture and decapitated. The brain was removed and placed into chilled artificial cerebrospinal fluid (ACSF, equilibrated with 95% O2, 5% CO2), then trimmed with a razor blade, glued to the stage of a vibratome, and sectioned into 400 μm thick slices. Following sectioning, brain slices were dissected further to free the hippocampus from surrounding structures. To prevent spontaneous bursting in the presence of the GABAA receptor antagonist bicuculline (see below), area CA3 was removed by a small cut.

Hippocampal slices were stored in a holding chamber at the interface of ACSF and humidified 95% O2, 5% CO2, at room temperature. Individual slices were transferred, as needed, to a small volume (~200 μl) interface-recording chamber (35°C), where they were perfused at 1 to 1.5 ml/min with ACSF. The upper surfaces of the slices were exposed to warmed, humidified atmosphere consisting of 95% O2, 5% CO2.

2.1.2 Electrophysiology

Intracellular recordings were made using the “blind” whole cell recording technique (Blanton, Lo Turco, & Kriegstein, 1989; Grover, 1998; Grover & Chen, 1999). The pipette solution was composed of 140 mM potassium gluconate, 10 mM Na-HEPES, 3 mM MgCl2, pH 7.2, 280–290 mOsm. Potentials were recorded with an Axoclamp 2B (Axon Instruments). Membrane potentials for the F-5 thru F-11 cells were amplified 10 times, digitized at 100 kHz, stored on a PC, digitally low-pass-filtered at 10 kHz, and digitally downsampled to 20 kHz. The rest were low-pass-filtered at 3 kHz, amplified 10 times, digitized at 20 kHz, and stored on a personal computer. Access resistance (compensated using the Axoclamp bridge circuit) and membrane resistance were monitored throughout the recordings. Recordings were terminated if large or abrupt changes in access or membrane resistance occurred.

Recordings of voltage at the soma were made in response to 10 types of stimuli delivered at 5 s intervals: short (2 ms) pulses of +750 pA and +1500 pA, medium (20 ms) pulses of −50 pA and −200 pA, and long (500 ms) pulses of −200 pA, −50 pA, +100 pA, +200 pA, +400 pA, and +800 pA. The complete set was delivered repeatedly at least five times. Only the 500 ms +800 pA, −50 pA, and −200 pA traces were used in the analysis presented in this report.

Recordings were made in the presence of three different media; 11 cells in ACSF (F group), 3 cells with the addition of 10 μM bicuculline, and 1μM CPG-55485 to block GABAA and GABAB currents, respectively (G group) and 4 cells with both GABA inhibitors plus 50 μM D-AP5 to block NMDA currents (A group). Test stimulation was delivered until a minimum 5 minute period of stable recording was achieved. All membrane potentials were corrected for the calculated liquid junction potential (−10 mV).

2.1.3 Cell Selection Criteria

Cells were accepted for analysis if their resting membrane potentials were at least −60 mV (after correction for the liquid junction potential) and their input resistance estimated from the voltage change at the end of the 200 ms current pulses was at least 30 MΩ.

2.2 Distance Functions

Comparison of traces was carried out according to the following definitions. The waveform distance Dw between traces Ta (t) and Tb (t) was given by

Dw(Ta,Tb)=1te[0teTa(t)Tb(t)pdt]1/p,
(2.1)

where te corresponds to the time of the last point under consideration (assumed the same for both traces) and p ≥ 1 is a parameter typically set at 1 or 2. When p = 2 the waveform distance corresponds to the commonly used sum of squared differences normalized to trace length.

The fiducial-point distance Dfp was given by

Dfp(Ta,Tb)=1te[i=0Ns0l¯iTa(miat+qia)Tb(mibt+qib)pdt]1/p,
(2.2)

where Ns = min(Na, Nb), Nk is the number of spikes in trace k [set membership]{a, b}, qik are the fiducial points for trace k with qik=tik(0<i<Ns),q0k=0,qNsk=max(tNsa,tNsb),qNs+1k=te,tik is the time of the ith spike of trace Tk, mik=lik/l¯i is the stretch factor for interspike interval i of trace k, lik=qi+1kqik is the duration of the ith interspike interval of trace k, and l¯i=(lia+lib)/2 is the average of the ith interspike interval durations of traces a and b. Note that Dfp reduces to Dw if all fiducial points agree.

More informally, in calculating the fiducial point distance, we used the time points of action potentials as fiducial (i.e., standard reference) points to break each trace into segments, which were then matched up, one by one, with those of the other trace for comparison (see Figure 1). Each pair of segments was then linearly stretched in time to the average of their durations, and their differences computed and aggregated over all matched pairs and normalized to trace length.

Figure 1
Schematic of the details for computing the fiducial point distance between traces Ta (t) and Tb (t). See equation 2.2 for labeling definitions.

Note that the time point of the next-to-last fiducial point qNsk=max(tNsa,tNsb) is the same for both traces; it corresponds to a spike peak time in one trace that in general is not a spike peak time of the other trace. Thus, the last interval is of equal length on both traces, so that for this interval, the formula reduces to computing the conventional Dw distance. Note that when traces have different numbers of spikes (NaNb), this last interval will have |NaNb | spikes in one trace and none in the other, so that the contribution of this interval to the Dfp distance will increase with increased discrepancies in the number of spikes between traces. This definition for qNsk was found to be important; alternative definitions had inferior multilevel nearest-neighbor performance (see below) and sometimes would produce results that violated the triangle inequality (see Table 2 and section 4).

Table 2
Scorecard of the Metric Conditions Satisfied by Each of the Distance Functions Used in This Study.

The interspike interval distance Dis was based on the differences between the interspike intervals and was given by

Dis(Ta,Tb)=1Ns[i=0Nslialibp]1/p,
(2.3)

and the spike time distance Dst was based on the differences between spike times and was given by

Dst(Ta,Tb)=1Ns[i=0Nstiatibp]1/p.
(2.4)

The phase plane distance Dph was given by (LeMasson & Maex, 2001)

Dph(Ta,Tb)=[ij|Ma(i,j)NspaMb(i,j)Nspb|p]1/p,
(2.5)

where the density matrix Mk (i, j) counts the number of times phase-plane points (v, v) of trace Tk (t) sampled at frequency Δt fall within the boundaries of a box of size Δv × Δv centered at (iΔv, jΔv), Nspk=[tek/Δt]+1 is the number of sampled points for trace k, and v = dv/dt is the first time derivative of v. Δv and Δv are empirically determined parameters; those used in the tables that follow were obtained by beginning with values similar to those reported in Achard and De Schutter (2006) and optimizing the multilevel nearest-neighbor classification results by twofold multiplications or divisions. The optimization process resulted in parameter values very different from those previously reported for Δv (0.05–0.1 versus 40 mV/ms); most of the time we report the optimized values but occasionally show the results for the values used in the above report (see Figure 2 and Table 1). Closer inspection of the Δv discrepancy with previously used values suggests that improved classification with optimized parameters is due to fine resolution in the region of spike initiation and regeneration, where the enhanced Δv resolution makes a difference. However, the fine Δv grid results in very little overlap between traces in other regions of the Mk density matrix. A variable resolution Δv may be able to reconcile the demand for fine resolution in the region of spike initiation with a coarser grid elsewhere.

Figure 2
Sample of current-clamp recording sets. All panels show sets of traces of current clamp responses to 500 ms pulses of +800 pA, −200 pA, and −50 pA. (A) Set of five current clamp responses for the cell with smallest (A1; ...
Table 1
Classification Errors for Traces Evoked by 500 ms Pulses of +800 pA.

More informally, to compute the phase-plane distance, a plot of the trace’s magnitude versus its first time derivative (known as a phase-plane in physics) is obtained for each trace, and the density of overlap between the plots is obtained. The trace overlap density is calculated by placing a grid over the plot and comparing the numbers of sampled points in each grid box. Spike trains in the phase-plane typically appear as repeating closed loops, one for each spike.

The spike time alignment distance Dspike was given by (Victor & Purpura, 1997),

Dspike[q](Ta,Tb)=minS0,S1,SrK(Sj,Sj1),
(2.6)

where S0 is the sequence t1a,t2a,,tNaa of spike times in trace Ta, Sr is the corresponding sequence of spike times in trace Tb, and K (Sj, Sj−1) is the cost of an elementary transformation of sequence Sj−1 to Sj. The allowable elementary transformations are adding or deleting a single spike with a cost of one, and adding Δt to a spike time with a cost of qt|, where q ≥ 0 is a parameter with units 1/s.

The interval time alignment distance Dinterval is similar to that of Dspike (Victor & Purpura, 1997), except that the sequences consist of interspike intervals; S0 is the sequence of interspike intervals d1a,d2a,,dNa+1a in trace Ta, where d1a=t1a,dia=tiati1a(1<iNa),dNa+1a=tetNaa, and Sr is the corresponding sequence of interspike intervals in trace Tb. The results reported below are for the best result obtained for q in 0, 0.01, 0.02, 0.05, 0.1, 0.2, 0.5, 1, 2, 5, 10, 50, and 100 1/s.

More informally, the interval alignment distances are based on the minimum “cost” of transforming the spike train sequence (spikes times or interspike intervals) of one trace to that of the other. The parameter q represents the cost of moving a spike (or extending an interval) t seconds in relation to the cost of adding or deleting a spike altogether. When q = 0 both alignment distances reduce to the difference in the number of spikes. (See Victor, 2005, for examples and applications of these metrics.)

Note that the units of Dw and Dfp are mV, those of Dis and Dst are ms, and Dph, Dspike, and Dinterval are dimensionless.

2.3 Clustering and Classification

Agglomerative hierarchical clustering techniques are useful to explore for the presence of clusters in the data set, though their quantitative assessment of cluster validity is limited. Hierarchical dendrograms were obtained using agglomerative hierarchical clustering with Ward’s method as described in Everitt (1974) using Lance and Williams’s flexible method (1966). Agglomerative methods begin by computing a distance matrix between all pairs in the set and then fuse individuals, or groups of individuals, using a similarity criterion. Ward’s criterion minimizes the loss of information that results from the grouping of individuals into clusters as given by the total sum of squared deviations of every point from the mean of the cluster to which it belongs.

To quantitatively assess the presence of trace cell clusters, we estimated the proximity and separation of traces from each given cell by counting how often their nearest neighbors were recordings from the same or different cells. More formally, multilevel nearest-neighbor classification was performed by computing the distance of test traces to all other traces in the set and sorting the results. A correct classification at level i (i < number of traces per cell) was counted if the provenance of the first through ith nearest neighbors all agreed with the provenance of the test trace.

Good scores on the multilevel nearest-neighbor classification test indicate a significant degree of clustering in the data set by providing an estimate of classifier accuracy over all possible training and testing data set partitions. The accuracy of classifiers is frequently measured by partitioning the available data into two sets; one is used for training the classifier, while the classifier’s performance on the other is reported as an estimate of classification accuracy. To improve reliability, many studies report aggregate results for a selection of partitions. Ideally, the performance on all possible partitions is sought, though in practice this can be computationally infeasible as the number of possible partitions is generally very large. The accuracy of one of the simplest classifiers—the nearest-neighbor classifier—on a data set (n samples per class) with a perfect score on the first level of the described multilevel nearest-neighbor classification test would be perfect on all possible partitions of the data set into n − 1 training samples per class and one testing sample per class. Similarly, a data set with a perfect score at level n − 1 on the multilevel nearest-neighbor classification test translates into perfect nearest-neighbor classifier performance on all possible partitions with one training sample per class and n − 1 testing samples per class.

3 Results

3.1 Characteristics of Cells Studied

Figure 2 shows voltage recordings from a selection of CA1 pyramidal cells used in the study. Responses were evoked by 500 ms current pulses of −200, −50, and +800 pA in sequence; five sets were obtained approximately 60 s apart. Hyperpolarizing currents evoked a characteristic “sag” originating from cationic h currents; we estimated the input resistance at the peak of the response (RNpk) and at the end of the pulse after reaching a steady-state amplitude (RNe). Depolarizing currents evoked a sustained train of action potentials that, as time progressed, widened, decreased in amplitude, and had longer interspike intervals (often referred to as accommodation).

Figure 2A shows all five traces from the cell that had the smallest (left) and largest (right) average fiducial point distance (Dfp, p = 1) between traces evoked by depolarizing pulses. Note that even in the best of cases, the agreement in spike shape and timing deteriorates as the train progresses.

The following panels illustrate the variability of traces evoked from different cells. Figure 2B shows the traces of cell pairs that had the smallest (left) and largest (right) depolarizing Dfp scores. Note the remarkable similarity between the depolarizing traces on the left, where not only spike timing but also spike shape are virtually identical, though at close inspection, spike initiation waveforms past about 150 ms differ slightly. Figure 2C shows traces that had the smallest Dph score for two sets of parameters of this measure. Note that the selection of parameters has a significant impact on Dph scores and on which traces are found more similar.

As noted in section 2, the cells in the study were recorded in three different media; 11 in ACSF (F group), 3 with blockers for GABAA and GABAB receptors (G group), and 4 with both GABA inhibitors plus an NMDA channel blocker (A group). One-factor ANOVAs for 20 measurements, including the resting potential, the peak and steady-state input resistance for the large and small pulse hyperpolarizing traces, and a number of features of the spike trains of depolarizing traces, including the amplitude, width, and interspike intervals for the first, second, and third spike, were computed to determine if the measures were systematically affected by the medium’s composition. The results were not significant in all cases, though there was a tendency for the A group to have hyperpolarized resting potentials and reduced spike widths when compared to the F group (F = 3.441, P = 0.06). As there were no statistical differences between groups, and given that these blockers are considered to affect primarily synaptic transmission and intercellular activity is not thought to be involved in the origin of the traces under study (which were evoked by intracellular current injection), we did not differentiate further between the three groups in the rest of our analysis.

Figure 3 shows box-and-whisker plots for the resting potential, the input resistance of hyperpolarizing traces, and a number of features of the spike trains of depolarizing traces. The median resting potential was −71 mV, ranging from −65 to −82 mV. “Voltage sag” caused input resistances to decrease at hyperpolarized voltages, indicating the presence of cationic currents. Input resistance estimates at larger currents (RNpk= 59 MΩ, RNe= 48 MΩ, −200 pA) were invariably smaller than at smaller currents (RNpk= 69 MΩ, RNe= 56 MΩ, −50 pA).

Figure 3
Box-and-whisker plots of selected measurements from the set of recordings evoked by current pulses in 18 CA1 pyramidal cells. The upper row shows the resting potential before the pulse (left) and the input resistance for 500 ms, −200 pA, and −50 ...

For depolarizing traces, the median number of spikes during the 500 ms pulse was 28, ranging from 11 to 33. Median peak amplitude for the first spike was 40 mV, which was reduced for subsequent spikes to 20 mV by the 20th spike. In contrast, the median spike width (at −20 mV) increased from 0.9 ms for the first spike to 1.5 ms by the 20th spike. Similarly, median interspike intervals increased from 5.1 ms for the first interval to 20 ms by the 20th interval. The measure that showed the most increase in variability as the train progressed was the interspike interval. The variability in spike width and amplitude remained roughly the same after the third spike. For all the measurements in Figure 3, the variability for recordings across cells was always larger than the variability of repeated recordings from the same cell (not shown).

3.2 Classification of Depolarizing Traces

To test if responses from a given cell are similar to each other yet distinct from those of other cells, we first used agglomerative hierarchical clustering to assess the existence of cell clusters and then proceeded with a quantitative analysis based on nearest-neighbor classification.

To find the hierarchical relationships, we computed the distances for all possible pairs of traces associated with a given stimulus protocol and constructed dendrograms using agglomerative techniques with Ward’s method (see section 2). The result of applying this procedure to the set of 90 traces—five +800 pA traces from each of 18 neurons—when using the fiducial point distance Dfp is shown in Figure 4 (top). Note that all five traces from each cell join at low height levels before joining traces from other cells. This suggests that responses from individual cells do cluster and that fiducial point distances can be used to systematically distinguish responses from one neuron from those of other neurons. The asterisks correspond to the depolarizing traces shown in Figure 2; note that traces from cell A1 join at very low levels, while those from F1 join at slightly higher levels, reflecting the distinctions highlighted in Figure 2A.

Figure 4
Dendrograms obtained by hierarchical clustering. Agglomerative hierarchical clustering using Ward’s method was applied to the distance matrix that resulted from computing fiducial point distances (p = 1, top), or phase-plane distances (p = 1, ...

The results when using other distance functions did not produce dendrograms where responses from the same cell would separate clearly from those of other cells. For example, the bottom of Figure 4 shows the dendrogram obtained when using the phase-plane distance with optimized parameters. Note that even though the traces for the first 5 cells on the left and the last 9 cells on the right, for a total of 14, show separate clusters, there is an overlap of the responses for cells A2 and G1 and for cells A1 and G3. Similarly deficient results were obtained with all other distance functions described in section 2. For example, the phase-plane distance with unoptimized parameters (p = 2, Δv = 2, Δv = 40) showed overlap on 7 cells and thus had only 11 cells that showed separate clusters (not shown).

While the dendrograms suggest that cell clusters exist, they do not provide a quantitative measure of clustering. To quantify the accuracy with which fiducial point and other distance measures can be used to identify the provenance of neuron traces, we used nearest-neighbor classification at multiple levels (see section 2). Distances from each trace to all other traces were computed and sorted. For each tested trace we determined, among the remaining 89 traces, the provenance of the first, second, third, and fourth nearest neighbors. For each level, a correct classification was counted if all trace provenances up to that level agreed with the provenance of the trace being tested. For example, a correct classification at the fourth level would require that the traces of the first through fourth nearest neighbors come from the same cell as the test trace.

Table 1 shows the results of applying this procedure to the set of 90 +800 pA traces using fiducial point distances and a number of distance functions that have been described in the literature. Fiducial point distance classification with p = 1 had the best performance: it was perfect at the first and second levels, had four errors at level three, and 10 errors at the fourth level. Thus, for any trace, 100% of its two nearest neighbors, 96.6% of its three nearest neighbors, and 88.9% of its four nearest neighbors agreed in their provenance, indicating that when using fiducial point distances, traces from each cell are closely grouped together and are distinct from those of other cells. Phase-plane distance performed second best (p = 1). The worst performers were those based on the timing of spikes, though the spike Dspike and interval Dinterval alignment measures did better than the absolute spike times Dst and interspike intervals Dis. All distance functions did better in the p = 1 than the p = 2 variant.

3.3 Distance Functions and Conditions for a Metric Space

A match (or distance) function that assigns a numerical value to pairs of traces needs to satisfy four conditions to give rise to what in mathematics is known as a metric space: (1) nonnegativity, (2) identity, (3) symmetry, and (4) triangle inequality. Table 2 lists the results of evaluating the various distance functions from the perspective of the metric space conditions. The nonnegativity and symmetry conditions were satisfied by all tested functions. The waveform distance is the only distance function that satisfied the identity condition. For the others, when two traces are equal, their distance is zero, but the converse may not hold; it is possible to construct traces with zero distance that are distinct from each other. For example, two traces may agree on their spike times (so that Dst and Dis are zero) yet have very different waveforms between the spikes.

The waveform distances Dw, Dspike, and Dinterval have been shown mathematically to satisfy the triangle inequality (Victor & Purpura, 1997; Rynne & Youngson, 2000). Mathematical proofs for the other methods are beyond the scope of this letter. Instead we tested the triangle inequality for all possible x-y-z combinations (90 · 89 · 88 = 794, 880) in the distance tables constructed for the p = 1 variants in Table 1. Only one small violation was found in the fiducial point case (perhaps due to rounding error), and there were zero violations for the phase-plane case. There were, however, substantial violations for the functions based on spike times: 620 and 2181 for the interspike and spike time distances, respectively. Furthermore, we should note that failure to include the 1/p power in the distance definitions invariably produced large triangle inequality violations for Dw, Dfp, and Dph when p ≠ 1.

3.4 Classifying Depolarizing Traces by Combination of Distance Functions

We explored whether distance function combinations can be constructed to improve on the results so far presented. In particular, we considered the p = 1 variants of Dfpw = Dfp + kDw and Dphw = Dph + kDw where k > 0 is a mix factor (p = 1, Δv = 1, Δv = 0.05). We found no values of k that would make Dfpw improve on the Table 1 results for Dfp. However, k = 0.16 did result in improvements for Dphw over Dph and Dw; the number of multilevel misclassifications for Dphw was 0, 4, 15, 19. We found a slight improvement over the Table 1 results for Dfp for the three-way mixture Dfp + k1 Dph + kDw; when setting k = 1 and k1 = 0.01–0.1 we obtained a 0, 0, 3, 7 score. We should note that with the addition of Dw, when k > 0, Dfpw, Dphw and the three-way mix satisfy the second metric space condition (identity), and therefore may satisfy all four metric space conditions.

3.5 Classifying Hyperpolarizing Traces

The results of applying the multiple-level nearest-neighbor analysis to the set of 90 traces (five from each of 18 neurons) obtained when stimulation amplitude was set at −200 pA for 500 ms is shown in Table 3. Since hyperpolarizing currents did not evoke action potentials, we should note that Dw = Dfp.

Table 3
Multilevel Nearest-Neighbor Classification Errors for 90 Traces Elicited by 500 ms −200 pA Current Steps in CA1 Pyramidal Cells.

The p = 2 variant of Dw had one error at the first level and in general did better than the p = 1 variant. The results for the phase-plane distance Dph were comparable to those of the waveform distance Dw.

Neither the resting potential Vr nor the input resistances RNpk and RNe were good classification indices; their multilevel classification errors were substantial.

To isolate the contribution of trace shape to classification, we considered two cases. In the first, we subtracted the baseline recorded before the current pulse from all traces before computing the Dw distances. In all cases, these operations increased the number of errors. Application of the nearest-neighbor analysis to traces evoked by −50 pA resulted in more classification errors than in the −200 pA case (compare Tables 3 and and4),4), though the phase-plane distance performed better than the waveform distance. Subtraction of trace averages increased the number of errors by a large margin (not shown).

Table 4
Multilevel Nearest-Neighbor Classification Errors for 90 Traces Elicited by 500 ms −50 pA Current Steps in CA1 Pyramidal Cells.

3.6 Classifying Depolarizing and Hyperpolarizing Trace Pairs

To evaluate the degree of redundancy between hyperpolarizing and depolarizing traces, we tested whether the classification results were improved by considering a distance function based on scores from both kinds of traces. Some of the results of this analysis are shown in Table 5 where the distance functions consisted of adding the scores for the traces stimulated by +800 pA pulses to those stimulated by −200 pA or −50 pA.

Table 5
Multilevel Nearest-Neighbor Classification Errors for 90 Pairs of Traces, Elicited by 500 ms Steps of +800 pA and Either −200 pA or −50 pA.

In general, the combination of depolarizing and hyperpolarizing traces resulted in improved scores across the board; addition of −200 pA traces fared better than addition of −50 pA traces, and the p = 1 variants did better than their p = 2 counterparts. Fiducial point distances scored best; they had 0, 0, 1, 4 and 0, 0, 2, 9 scores when combined with the −200 pA and the −50 pA traces (p = 1 variant), respectively. The p = 1 variants of the phase-plane distances were second best. The worst scores were obtained for the p = 2 variants when the +800 pA phase-plane distances were added to waveform distances of hyperpolarizing traces. Combining the −50 pA and −200 pA waveform hyperpolarizing trace distances did not improve results over those obtained when isolated (not shown).

4 Discussion

The nature of the variability inherent in the electrophysiology of pyramidal cells in hippocampus and other brain regions has an impact on scientific issues that range from the construction of realistic neuron models to the design of plausible pharmacological interventions. For example, a model for long-term potentiation in a homogeneous neuronal network may fail in a significant portion of cells in a heterogeneous population where the proposed mechanism is ineffective; only if the proposed mechanism is robust in the face of the range of variability of the heterogeneous population would the candidate mechanism be viable. The same argument applies to the design of pharmacological interventions.

Electrophysiological recordings show substantial cell-to-cell variability; whether it arises from experimental or intrinsic causes has received little attention. In this study, we recorded multiple traces from each of several pyramidal cells in hippocampal CA1 and asked whether the variability we observed ranged uniformly over all cells or whether each cell had a consistent profile that varied among cells. Although all recorded waveforms shared general features, we found that each cell had a consistent response profile that could be systematically differentiated from that of other cells.

Traces arising from depolarizing currents were more readily differentiated than those arising from hyperpolarizing pulses, while traces evoked by large hyperpolarizing currents provided better cell discrimination than those evoked by small ones. Combining the scores from depolarized and hyperpolarized responses provided improved discrimination, indicating that the two stimulation conditions are not redundant and suggesting that cell-specific factors are represented in both positive and negative voltage excursions. The finding that adding synaptic receptor–related blockers to the recording medium did not have significant effects on the results indicates that the responses used in the study rely on intrinsic cell properties alone.

4.1 Measuring Waveform Similarities

Evaluation of match functions for traces with trains of action potentials has typically been based on their impact on some aspect in a parameter search of simulated target traces, such as the speed of convergence or the precision with which parameter values are estimated (Vanier & Bower, 1999; Keren et al., 2005; Achard & De Schutter, 2006; Weaver & Wearne, 2006). This approach fails to separate the effectiveness of the search method from that of the match function, making it difficult to interpret results.

Here we took a different approach. The criterion presented here is based on the notion that the design of the match function should be sensitive to experimentally found regularities. Thus, to the best extent possible, experimental traces from the same cell should be rated closer to each other than to those from other cells. Match functions that perform well under this criterion can be expected to perform well with parameter search methods, as the match function will give shape to a search space where trace variants known to come from a neuron presumably with identical parameter values will form a neighborhood separate from traces known to come from other cells with presumably different parameter values.

Measuring the similarity between recorded waveforms is essential to detecting cell-specific profiles, and it is important to have a method that can do this job well. While the sum of squared differences (Dw above) is adequate when comparing most traces, this method is inadequate when action potentials are present. The narrow shape and subtle variation in peak times of action potentials would result in a poor score for a pair of traces shifted by a spike width but otherwise identical.

To address this issue, we developed a novel match function based on fiducial (i.e., standard-reference) points (see section 2 and Figure 1). Fiducial point matching Dfp outperformed all other methods in identifying cell-specific profiles. Phase-plane matching Dph, where traces are transformed to the voltage-gradient phase plane for comparison, performed second best, even after extensive tuning of its two parameters (Δv, Δv). Measures based on the timing of spikes were the worst performers.

Interestingly, the p = 1 variants outperformed the p = 2 variants most of the time, suggesting that the common practice of squaring differences may require some scrutiny. In general, when constructing a similarity metric analogous to equation 2.2, the best choice of p may be domain specific. Insofar as 1 ≤ p < ∞ is satisfied, there is no mathematical consideration in measure theory to favor any particular value. A basic difference between the two variants is that differences that are smaller than one get much smaller after squaring when compared to taking their absolute value. It may be that these small differences play an important role in the waveforms of this study.

A match (or distance) function that assigns a nonnegative numerical value to pairs of elements in a set gives rise to a mathematical space, the nature of which varies depending on the definition of the match function. For this space to become a metric space, it has to satisfy the four conditions given in Table 2. For example, the x-y plane with D(p,q)=(pxqx)2+(pyqy)2 is a typical metric space. Failure to satisfy any of these four conditions (or combination of them) results in spaces that do not correspond to our intuition and, furthermore, may preclude application of common mathematical tools. Satisfaction of the metric space conditions is considered to be so basic that it is hardly mentioned at all in search method publications. However, most parameter search methods have been developed with the assumption that metric space properties are satisfied, and their rates of convergence, if not the definition of convergence itself, typically rely on having a metric space.

We found that only the waveform distance Dw satisfied all four conditions of a metric space; all the rest failed the identity condition. The triangle inequality was satisfied by the spike alignment measures Dspike and Dinterval, perhaps by the fiducial point and phase-plane distances, but not the spike-based distances. Nevertheless, the identity condition can be satisfied by adding the waveform distance to the fiducial and phase-plane distances. This improved the scores of the phase-plane distance function and did not perturb the already very good performance of the fiducial point distance function. Some published reports have used distance function definitions that violate the metric properties; omission of the 1/p power results in triangle inequality violations (Keren et al., 2005; Achard & De Schutter, 2006), and separate treatment of model and target data violates the symmetry condition (Vanier & Bower, 1999).

4.2 Origin of Cell-Specific Electrophysiological Profiles

Our findings suggest that CA1 pyramidal cells have unique properties that confer them with a degree of individuality in electrophysiological terms. These unique properties can arise from cell-to-cell differences in morphology, passive electrotonic properties, or active properties (densities and kinetics of membrane channels).

The best-documented source of cell-to-cell variability is morphology. Morphological reconstructions of CA1 pyramidal cells (Ishizuka, Cowan, & Amaral, 1995; Carnevale, Tsai, Claiborne, & Brown, 1997; Pyapali, Sik, Penttonen, Buzsaki, & Turner, 1998; Megias, Emri, Freund, & Gulyas, 2001) show considerable variation, though there is reason to believe that as a whole, these reports overestimate intrinsic morphological heterogeneity due to reconstruction issues that have significant impact on response profiles (Scorcioni, Lazarewicz, & Ascoli, 2004; Szilagyi & De Schutter, 2004; Ambros-Ingerson & Holmes, 2005).

Passive electrotonic properties may vary both within a cell and between cells. Direct measurement of Cm in four types of neuron or glia has yielded average values ranging from 0.85 to 1.06 (Gentet, Stuart, & Clements, 2000). Values of Ri are likely to be sensitive to the presence of large proteins in organelles and in the cytoplasmic matrix, and thus may vary within and between cells (Rall et al., 1992). Membrane resistivity (Rm) is a term used to aggregate the effects of all voltage-independent ionic conductances. The variability in resting potentials observed in this and other studies may have its origin in cell-to-cell variability in Rm.

The finding that cells were better differentiated by depolarizing pulses than by hyperpolarizing pulses suggests that active properties are important sources of cell-to-cell variability. The list of voltage-activated channels currently thought to generate the currents that give rise to the recordings in our study include the fast and persistent Na+ currents, the K+ delayed rectifier and A currents, Ca2+ currents, Ca2+-mediated K currents, and the hyperpolarizing-activated cationic h current (Borg-Graham, 1999). Their interaction is complex and believed to be highly nonlinear, and small variations may result in discernible profiles.

Some insight into the plausible cell-to-cell variability of active parameters can be gleaned from the range of values reported for a common kinetic variable such as the V1/2 of activation of the fast Na+ channel in hippocampal pyramidal cells; reports go from −55 mV (Steinhauser, Tennigkeit, Matthies, & Gundel, 1990), through −38 mV (Sah et al., 1988) to −30 mV (Magee & Johnston, 1995). A change of this magnitude would produce radically different results in any computational neuron model. While some of the variability may be attributable to differences in methodology such as rat strain or recording method, the contribution of cell-to-cell variability may be significant.

On the other hand, even large variations in channel densities may result in traces with similar properties; in a remarkable study, different ratios of inward Na+ and Ca2+ currents were found at the same membrane potential in Purkinje neurons with almost identical activity profiles (Swensen & Bean, 2005). These findings support the notion that substantial coordinated variations in cell-specific factors may be a common occurrence in neurons.

4.3 Implications for Hippocampal CA1 Pyramidal Cell Models

The results from this study suggest that the electrophysiological responses from CA1 pyramidal cells, while having many common features, may not be homogeneous, and that each neuron has a specific response profile. The existence of cell-specific response profiles and the ability to identify them may make it possible to discern whether the underlying factors in the electrophysiology of CA1 pyramidal cells are homogeneous or heterogeneous. Homogeneity predicts that underlying factors such as channel densities and kinetic parameters have small variations around mean values; heterogeneity predicts clustered variations and, possibly, large or abrupt coordinated variation among diverse factors.

The available evidence is insufficient to apportion the relative contribution of passive or active mechanisms to cell-specific profiles. Elucidating the specific mechanisms that modulate cell-specific profiles, and their relative importance, does not appear to be straightforward. A possible avenue is to construct detailed computational models to match the waveforms of several cells and compare the parameter values that achieve this. Automated parameter search methods that perform well in a very large space with many local minima, are able to cope with a large number of parameters, and are parallelizable will surely be needed. Our results indicate that matching several traces at once is likely to improve accuracy and that using the sum of fiducial point Dfp and waveform Dw distances is the most promising match function to compute and minimize the distance between simulated and experimental traces.

4.4 Experimental Database Resource

The experimental recording data set from this study will be made available at axon.zool.ohiou.edu##.

Acknowledgments

This work was supported by NIH grant AA-014294 and the NSF-NIH CRCNS program. It was aided in part by an allocation of computing time from the Ohio Supercomputer Center.

Contributor Information

José Ambros-Ingerson, Department of Biological Sciences, Neuroscience Program and Quantitative Biology Institute, Ohio University, Athens, OH 45701, U.S.A.

Lawrence M. Grover, Department of Pharmacology, Physiology and Toxicology, Marshall University School of Medicine, Huntington, WV 25704, U.S.A.

William R. Holmes, Department of Biological Sciences, Neuroscience Program and Quantitative Biology Institute, Ohio University, Athens, OH 45701, U.S.A.

References

  • Achard P, De Schutter E. Complex parameter landscape for a complex neuron model. PLoS Comput Biol. 2006;2:e94. [PubMed]
  • Ambros-Ingerson J, Holmes WR. Analysis and comparison of morphological reconstructions of hippocampal field CA1 pyramidal cells. Hippocampus. 2005;15:302–315. [PubMed]
  • Ambros-Ingerson J, Grover LM, Holmes WR. Abstract viewer itinerary planner. Washington, DC: Society for Neuroscience; 2005. A parameter search method for optimal matching of impulse firing waveforms from hippocampal CA1 pyramidal cells.
  • Blanton MG, Lo Turco JJ, Kriegstein AR. Whole cell recording from neurons in slices of reptilian and mammalian cerebral cortex. J Neurosci Methods. 1989;30:203–210. [PubMed]
  • Borg-Graham LJ. Interpretations of data and mechanisms for hippocampal pyramidal cell models. In: Ulinski PS, Jones EG, Peters A, editors. Cerebral cortex: Models of cortical circuits. New York: Plenum; 1999. pp. 19–138.
  • Carnevale NT, Tsai KY, Claiborne BJ, Brown TH. Comparative electrotonic analysis of three classes of rat hippocampal neurons. J Neurophysiol. 1997;78:703–720. [PubMed]
  • Everitt B. Cluster analysis. New York: Halsted Press; 1974.
  • Gentet LJ, Stuart GJ, Clements JD. Direct measurement of specific membrane capacitance in neurons. Biophys J. 2000;79:314–320. [PubMed]
  • Grover LM. Evidence for postsynaptic induction and expression of NMDA receptor independent LTP. J Neurophysiol. 1998;79:1167–1182. [PubMed]
  • Grover LM, Chen Y. Evidence for involvement of group II/III metabotropic glutamate receptors in NMDA receptor-independent long-term potentiation in area CA1 of rat hippocampus. J Neurophysiol. 1999;82:2956–2969. [PubMed]
  • Hoffman DA, Magee JC, Colbert CM, Johnston D. K+ channel regulation of signal propagation in dendrites of hippocampal pyramidal neurons. Nature. 1997;387:869–875. [PubMed]
  • Ishizuka N, Cowan WM, Amaral DG. A quantitative analysis of the dendritic organization of pyramidal cells in the rat hippocampus. J Comp Neurol. 1995;362:17–45. [PubMed]
  • Keren N, Peled N, Korngreen A. Constraining compartmental models using multiple voltage recordings and genetic algorithms. J Neurophysiol. 2005;94:3730–3742. [PubMed]
  • Lance G, Williams W. Computer programs for hierarchical polythetic classification. Computer Journal. 1966;9:60–64.
  • LeMasson G, Maex R. Introduction to equation solving and parameter fitting. In: De Schutter E, editor. Computational neuroscience: Realistic modelling for experimentalists. London: CRC Press; 2001. pp. 1–25.
  • Madison DV, Nicoll RA. Control of the repetitive discharge of rat CA 1 pyramidal neurones in vitro. J Physiol. 1984;354:319–331. [PubMed]
  • Magee JC. Dendritic hyperpolarization-activated currents modify the integrative properties of hippocampal CA1 pyramidal neurons. J Neurosci. 1998;18:7613–7624. [PubMed]
  • Magee JC, Johnston D. Characterization of single voltage-gated Na+ and Ca2+ channels in apical dendrites of rat CA1 pyramidal neurons. J Physiol. 1995;487(Pt 1):67–90. [PubMed]
  • Mainen ZF, Sejnowski TJ. Reliability of spike timing in neocortical neurons. Science. 1995;268:1503–1506. [PubMed]
  • Marder E, Goaillard JM. Variability, compensation and homeostasis in neuron and network function. Nat Rev Neurosci. 2006;7:563–574. [PubMed]
  • Megias M, Emri Z, Freund TF, Gulyas AI. Total number and distribution of inhibitory and excitatory synapses on hippocampal CA1 pyramidal cells. Neuroscience. 2001;102:527–540. [PubMed]
  • Prinz AA, Bucher D, Marder E. Similar network activity from disparate circuit parameters. Nat Neurosci. 2004;7:1345–1352. [PubMed]
  • Pyapali GK, Sik A, Penttonen M, Buzsaki G, Turner DA. Dendritic properties of hippocampal CA1 pyramidal neurons in the rat: Intracellular staining in vivo and in vitro. J Comp Neurol. 1998;391:335–352. [PubMed]
  • Rall W, Burke RE, Holmes WR, Jack JJ, Redman SJ, Segev I. Matching dendritic neuron models to experimental data. Physiol Rev. 1992;72:S159–186. [PubMed]
  • Rynne BP, Youngson MA. Linear functional analysis. London: Springer; 2000.
  • Sah P, Gibb AJ, Gage PW. The sodium current underlying action potentials in guinea pig hippocampal CA1 neurons. J Gen Physiol. 1988;91:373–398. [PMC free article] [PubMed]
  • Schulz DJ, Goaillard JM, Marder E. Variable channel expression in identified single and electrically coupled neurons in different animals. Nat Neurosci. 2006;9:356–362. [PubMed]
  • Scorcioni R, Lazarewicz MT, Ascoli GA. Quantitative morphometry of hippocampal pyramidal cells: Differences between anatomical classes and reconstructing laboratories. J Comp Neurol. 2004;473:177–193. [PubMed]
  • Steinhauser C, Tennigkeit M, Matthies H, Gundel J. Properties of the fast sodium channels in pyramidal neurones isolated from the CA1 and CA3 areas of the hippocampus of postnatal rats. Pflugers Arch. 1990;415:756–761. [PubMed]
  • Storm JF. Action potential repolarization and a fast after-hyperpolarization in rat hippocampal pyramidal cells. J Physiol. 1987;385:733–759. [PubMed]
  • Storm JF. Potassium currents in hippocampal pyramidal cells. Prog Brain Res. 1990;83:161–187. [PubMed]
  • Swensen AM, Bean BP. Robustness of burst firing in dissociated Purkinje neurons with acute or long-term reductions in sodium conductance. J Neurosci. 2005;25:3509–3520. [PubMed]
  • Szilagyi T, De Schutter E. Effects of variability in anatomical reconstruction techniques on models of synaptic integration by dendrites: A comparison of three Internet archives. Eur J Neurosci. 2004;19:1257–1266. [PubMed]
  • Vanier MC, Bower JM. A comparative survey of automated parameter-search methods for compartmental neural models. J Comput Neurosci. 1999;7:149–171. [PubMed]
  • Victor JD. Spike train metrics. Curr Opin Neurobiol. 2005;15:585–592. [PMC free article] [PubMed]
  • Victor JD, Purpura KP. Metric-space analysis of spike trains: Theory, algorithms and application. Network: Computation in Neural Systems. 1997;8:127–164.
  • Weaver CM, Wearne SL. The role of action potential shape and parameter constraints in optimization of compartmental models. Neurocomputing. 2006;69:1053–1057.