|
|
||||||||
INNOVATIVE METHODOLOGY
1Max Planck Institute for Biological Cybernetics, Tübingen, Germany, 2Department of Neuroscience, Baylor College of Medicine, Houston, Texas; and 3Division of Biology, Division of Engineering and Applied Science, California Institute of Technology, Pasadena, California
Submitted 7 March 2007; accepted in final form 15 October 2007
|
|
ABSTRACT |
|---|
|
|
|
INTRODUCTION |
|---|
|
1) A prevailing hypothesis about the formation of memories is that they are gradually established across distributed neocortical circuits (Squire et al. 2004
). Therefore, having the ability to track the responses of the very same cells over days and weeks across different parts of the brain is of paramount importance in studying the mechanisms of learning. In addition, the ability to monitor the properties of individual neurons will also enable us to determine the type of neurons that are more prone to change their properties during learning. For instance, neurons that are broadly tuned may be more plastic than those that are already sharply tuned.
2) A fundamental component of learning is the ability to generalize and make useful predictions under novel situations not encountered before. This aspect of learning distinguishes it from mere memory encoding and lies at the core of intelligence (Poggio and Bizzi 2004
). Being able to understand the reorganization principles with which learning is accomplished across neuronal networks may also help us gain insights into how generalization to novel instances is accomplished and the reasons for limitations in generalization (Fahle 2005
; Fahle and Edelman 1993
). Yet, to study the principles of generalization of learning, one needs to extensively study the properties of neuronal circuits under a wide variety of different conditions. Given the limitations of how many trials a monkey can execute each day such detailed studies are difficult to achieve during a single recording session. Chronic stable recordings can thus also be of great practical importance.
3) The information contained in the activity of a neuronal ensemble critically depends on the interactions between its constituent neurons in addition to the response tuning functions of each of these individual neurons. Accordingly, changes in the interactions among neurons may be a principal mechanism for increasing the information content of a neuronal ensemble (Stopfer and Laurent 1999
; Stopfer et al. 1997
). Multiple chronically implanted electrode arrays can be used to characterize how the interactions in a neural circuit change during learning and memory consolidation.
Previous studies have reported recording stability across days based on similarity between action potential waveforms (Greenberg and Wilson 2004
; Rousche and Normann 1998
; Schmitzer-Torbert and Redish 2004
; Taylor et al. 2002
; Wilson et al. 2003
) or waveform features (Nicolelis et al. 2003
). Unfortunately, most of these studies are based on the premise that observing similar waveforms implies recording from the same neurons. However, it is not clear a priori how similar the waveforms of different adjacent neurons recorded from the same electrode are. This needs to be estimated empirically and be the foundation of any assessment of recording stability. To this end, we developed a statistical framework to quantify the accuracy with which stable neurons can be identified across successive recording days and weeks. Using this framework we demonstrate the feasibility of recording from the same neurons across days. We recorded data using the technique of multiple, chronically implanted tetrodes (Gray et al. 1995
; Wilson and McNaughton 1993
), which we adapted and further developed for research in awake, behaving macaques.
|
|
METHODS |
|---|
|
Recording chambers were positioned stereotaxically based on high-resolution magnetic resonance anatomical images. These images were collected in a vertical 4.7 T scanner with a 40-cm-diameter bore (Biospec 47/40c; Bruker Medical, Ettlingen, Germany). The system had a 50 mT/m (180 µs rise time) actively shielded gradient coil (B-GA 26, Bruker Medical) of 26 cm diameter. We used a custom chair and custom system for positioning the monkeys in the magnet. Anatomical data were collected with T1-weighted high-resolution (256 x 256 x 160 real data points at 0.5 mm isotropic linear resolution) images using three-dimensional (3D) MDEFT (modified driven equilibrium Fourier transform) pulse sequences, with an echo time (TE) of 4 ms, repetition time (TR) of 22 ms, flip angle (FA) of 20°, and four segments. These anatomical scans were done while the animals were anesthetized. Morphological methods were used to extract the skull parameters (Paravision, Bruker Medical) and create a 3D rendered surface (Analyze; Mayo Foundation, Rochester, MN) to be used for designing the cranial headpost and the recording chamber to mirror-image the skull surface. These form-specific implants were built using a five-axis CNC machine (Willemin-Macodel W428, Bassecourt, CH Energy Group), which resulted in an excellent fit between the implants and the skull surface. The animal received titanium implants. Before implantation the implants were coated with hydroxyapatite (HA). In addition to the titanium chamber a form-specific acrylic-free histocompatible material with the highest available mechanical strength was custom designed, constructed, and used to replace the bone removed inside the chamber. Similar to the 3D custom design of the recording chamber, high-resolution MRI data were used to construct the bottom surface of a bone replacement plate to mirror-image the surface of the dura. Given the excellent fit between the titanium implants and skull, and the bone replacement plate and dura, the implants were secured only with screws without the use of any dental acrylic or bone cement. Titanium HA-coated screws were used to secure the implants on the skull.
Surgical methods
The cranial headpost, scleral eye coil, recording chambers, and chronic tetrode array and its related parts were implanted under general anesthesia using aseptic conditions. After the subcutaneous injection of Robinul (0.01 mg/kg) and Ketavet (15 mg/kg), we prepared the monkeys for intubation by intravenous injections of the analgesic fentanyl (0.003 mg/kg), the barbiturate anesthetic trapanal (5 mg/kg) and the paralytic lysthenon (3 mg/kg). Throughout the duration of the surgery the animal was maintained under balanced anesthesia consisting of 1.3% isoflurane.
The recording chamber was implanted at a location determined by stereotactic coordinates over the operculum in area V1. A craniotomy was performed inside the chamber followed by a duratomy in the center of the craniotomy coinciding with the desired recording location. A custom-designed and -built dura ring was inserted inside the duratomy. The main components of the dura ring are two thin plates inserted above and below the intact dura surrounding the duratomy. Thereafter the bone replacement plate was positioned inside the chamber and around the dura ring. Custom surgical tools were designed and built to guide the placement of the bone replacement plate while the dura ring was held into position. Finally, the tetrode array was positioned. Each tetrode could now be independently adjusted. Other than the tetrodes themselves nothing else (such as metallic guide tubes) was inserted into the cortex, thus avoiding unnecessary damage to the brain.
Chronic electrophysiological and eye movement recordings
Neural activity was recorded using custom-built tetrodes whose tips were electroplated with gold. No preselection functional criteria were applied for the neurons recorded. Neural activity was sampled at 32 kHz, digitized, and stored using the Cheetah data acquisition system (Neuralynx, Tucson, AZ). Whenever the recorded voltage exceeded a predefined threshold (typically 25 µV) on any tetrode channel, a segment of 32 samples was extracted across all four channels and stored for later use. We refer to these putative spikes as "events." After an event was detected, the triggering mechanism was disabled for a period of 0.5650 ms (18 samples) for the tetrode on which the event occurred to avoid double-triggering the same spike.
The animals were implanted with a scleral search coil (Judge et al. 1980
) and their eye movements were monitored on-line. Data were also collected for off-line analysis using both a QNX-based data acquisition system at 200 Hz and the Cheetah data acquisition system at 2,000 Hz.
Behavioral paradigm and visual stimulation
Visual stimuli were displayed using a dedicated graphics workstation (TDZ 2000; Intergraph Systems, Huntsville, AL) with a resolution of 1,280 x 1,024 and an 85 Hz refresh rate, running an OpenGL-based stimulation program. The behavioral aspects of the experiment were controlled using the QNX real-time operating system (QSSL, Ontario, Canada). After the monkey acquired fixation on a colored square target (0.2°) for 300 ms, a sine wave grating stimulus was presented (typical grating parameters were: size of grating, 5° in diameter; spatial frequency, four cycles per degree; 100% contrast in case of moving gratings; 16 different directions of motion equally distributed; and 0.5 cycles/s). The grating stimuli were displayed for 500 ms and the animal was required to maintain fixation for another 500 ms. At the end of each successful trial a drop of apple juice was used for reward. The fixation window was ±1°.
Clustering methods for single-unit isolation
OVERVIEW.
Our spike-sorting method relies on fitting a mixture of Gaussians (MoG) model to the data. After spike alignment, we perform PCA to reduce data dimensionality and fit the mixture model using the split and merge expectation maximization (SMEM) algorithm (Ueda et al. 2000
). To assess clustering quality we use a contamination measure based on the assignment certainty to the individual mixture components.
SPIKE ALIGNMENT.
The waveforms, consisting of 32 samples on each channel, are realigned to reduce jitter induced by the discretization. After up-sampling by a factor of 10 using cubic spline interpolation, every waveform is aligned to its center of mass, calculated from a continuous region that surpasses 50% of the peak value, and resampled accordingly (Sahani 1999
). Events with center of mass not lying within a window of ±2 points (±0.06 ms) around the 8th data point are excluded from spike sorting because they represent overlapping spikes and could hinder reliable assignment (typically between 5 and 10%). Therefore to account for the maximum allowable shift of ±2 points, the number of samples representing one event has to be reduced to 28.
FEATURE EXTRACTION.
We selected the first three principal components (PCs) on each channel, which typically account for around 80% of the total variance and 73% of the total energy on the average, as features. This projects the 112-dimensional data (28 samples on four channels) into a 12-dimensional subspace (three PCs on four channels). Although it was previously shown (Görür et al. 2004
) that dimensionality reduction and clustering can in fact be performed concurrently in a mixture of factor analyzers model, we found this approach computationally intractable for large data sets.
MODEL DESCRIPTION.
The MoG method is used to model the waveform clusters of each neuron. Waveforms are assumed to be constant over time but superimposed by normally distributed noise. For m components, the probability density function evaluated at point x is given by
![]() | (1) |
The mixing proportions
i represent the a priori probability of any point belonging to a certain cluster and are normalized to
i=1m
i = 1. Intuitively,
i is the fraction of events belonging to component i. Each mixture component is modeled by a multivariate normal distribution with mean µi and covariance matrix Qi. Its probability density function is given by
![]() | (2) |
For the purpose of clustering in general and spike sorting in particular, one seeks a unique assignment of data points to clusters. This is achieved by choosing the component that maximizes the posterior probability
i
(x; µi, Qi) of having generated the data point.
A widely made assumption is that most of the variability in recorded spike waveforms is due to low-amplitude background signals from distant neurons (Pouzat et al. 2002
; Sahani 1999
). Under this assumption, the covariance matrices of all mixture components are expected to be dominated by the noise covariance and thus highly similar. Furthermore, this common covariance matrix could be easily estimated from recorded noise. In practice, however, we found that using the same covariance matrix for all components (homoscedastic model) only poorly describes waveform variability—a fact that can be attributed to a number of reasons: First, spike waveforms produced by a certain neuron are not constant, and therefore the covariance matrices of different neurons may be significantly different, especially in the presence of a high signal-to-noise ratio. Second, the activity of individual neurons will have a different temporal relationship to local network activity. Accordingly, the background noise resulting from network activity will not be the same for all recorded neurons. Third, although jitter effects due to the discrete triggering process are reduced by the spike alignment we perform, it still affects the covariance matrix of a cluster in a way that depends on the waveform shape. For these reasons, we decided to allow different covariance matrices for all components (heteroscedastic model) to cluster the data.
One common interpretation of a mixture model enforces a one-to-one correspondence of mixture components and neurons. Two observations led us to digress from this restricted perspective. First, a notable fraction of triggered events is made up of background noise, which includes action potentials of distant neurons. In general, such events feature low amplitudes and evade satisfactory classification by any spike-sorting method. Thus it is sensible to introduce the notion of a noise cluster responsible for all spikes of unknown origin and groups of spikes that cannot be isolated. The characteristics of this collection of noise events render it intangible to be modeled by a single Gaussian, an observation easily attributed to the varying event sources. Consequently, we propose a mixture of normally distributed components as an appropriate model for the noise cluster.
The second issue is related to complex-spiking neurons (Gray et al. 1995
), which exhibit a decrease in spike amplitude from spike to spike within a burst. Small, overlapping clusters in peak-to-peak representation as well as in PC space are the hallmark of complex-spiking neurons alongside their distinct cross-correlograms. To accommodate the high variability of spike waveforms, we resort to modeling complex-spiking neurons by a mixture of normal components as well. Because of the difficulties associated with automatically identifying complex-spiking neurons, user intervention is necessary for this step to group mixture components together based on the preceding indicators. In our data set, bursting neurons of this sort were extraordinarily rare.
MODEL FITTING.
To fit a mixture model to the data, we use the SMEM algorithm (Ueda et al. 2000
). SMEM is less prone to converge to local maxima than the expectation-maximization (EM) algorithm because it repeatedly tries to split and merge components. In each of these steps, two mixture components Ci, Cj are merged into one and a different component Ck is split into two. When the normal EM steps have converged, a list of all potential triplets (i, j, k), which represent the merge of clusters Ci and Cj and the split of cluster Ck, is compiled and sorted by the potential gain in likelihood of the model. For detailed explanations of the selection process and the ensuing modifications during a split and merge step, we refer to Ueda et al. (2000)
.
As previously indicated, initialization is crucial for good performance of the EM algorithm (McLachlan and Peel 2000
). Instead of starting from many random initial values, we perform k-means clustering to provide an initialization of high quality, which avoids the necessity of repeated runs of the SMEM algorithm.
The second issue in the application of EM to fit mixture models is the identification of the correct number of components. It is possible to determine this by fitting the model on a training subset of the data and evaluate its performance in a cross-validation scheme based on the model likelihood on the left-out validation subset. Further methods rely on information criteria like the Bayesian Information Criterion (BIC) or the Akaike Information Criterion (AIC). In all these approaches, the model likelihood is penalized by the number of free parameters in some form. The BIC is defined as
![]() | (3) |
Although determining the model complexity and fitting the model are done fully automatically, because most of the triggered events are in the highly non-Gaussian noise cluster, the algorithm sometimes slightly overestimates the model complexity, which may result in multiple components being fit to a single neuron. In this generally rare case, we merge these components manually and rerun EM to obtain the final result. Note that in this refitting step no splitting and merging are performed to prevent overclustering artifacts from begin reintroduced into the model. Additionally, as described earlier, mixture components of complex-spiking neurons are grouped and clusters that do not show a pronounced furrow in their cross-correlogram (typical for well-isolated units because of the refractory period) are discarded into the noise cluster. For these manual steps we rely on a modified version of the Klusters software (L. Hazan, Buzsáki lab, Rutgers, Newark, NJ). Note that these manual steps differ fundamentally from manual cluster cutting in that their purpose is not to manually define cluster boundaries but rather to help the algorithm to determine the model complexity and to validate and consolidate its output.
Quantitative evaluation of single-unit isolation quality
To characterize the functional changes associated with learning at the level of single cells, it is imperative to have excellent single-unit isolation. The quality of single-unit isolation is also critical for studying the principles of neural coding such as identifying the existence of precise temporal spike patterns in the firing of single neurons (e.g., bursts) and more complex spatiotemporal spike patterns across neuronal assemblies. Spike clustering from tetrode recordings relies on using four recording channels, previously shown to provide superior single-unit isolation over using information from single wires (Gray et al. 1995
; Harris et al. 2000
).
We quantified the quality of isolation of each identified cluster in the following way. The clustering process yields a model with m components, characterized by their means µi, covariance matrices Qi, and priors
i for i = 1, ..., m. It can be used as a generative model to create synthetic data that follow the distribution in Eq. 1 and should be virtually indistinguishable from the recorded data set. Using the same classification function
![]() | (4) |

iN
samples from
(µi, Qi)]. Consequently, the index g(x)
{1, ..., m} of the component that generated the sample x is known for all x
S. Then the rate of false positives eiFP and false negatives ("misses") eiFN for component i is given by
![]() | (5) |
![]() | (6) |
We calculated cont(i) for all single units in our data sets and included only neurons that satisfy cont(i) <0.05.
We found our contamination measure to be in good agreement with the visualization of cluster separation in four-dimensional space (both peak-to-peak and principal components, shown in Figs. 1 and 3). However, there are some instances where the four-dimensional space is not sufficient to observe cluster separation. For example, we note that in the second row of Fig. 3 the yellow cluster does not exhibit a clear boundary to other clusters and background activity, although it is in fact very well isolated as reflected by our single-unit isolation criterion. This apparent discrepancy is explained by the fact that these cluster plots show projections only onto the first principal component axis of each channel (or peak-to-peak amplitude). However, clustering is performed in a 12-dimensional space derived from the first three principal component axes of each channel (see METHODS). In this case, the yellow cluster is separable from noise when spike waveform information, represented by the other principal component axes, is taken into account.
|
|
We evaluated single-unit stability by computing the distance between average waveforms of pairs of single units recorded on consecutive days. We will refer to this simply as the distance between two average waveforms. The distance is computed in the 28 x 4-dimensional space of aligned waveforms (see Clustering methods for single-unit isolation), where 28 is the number of samples for each average waveform on four channels.
The set of average waveforms of neurons recorded from a given tetrode on a given day is denoted by Cn = {X
28x4:X recorded on day n}. We denote by xi the average waveform on channel i (i.e., the ith column of the matrix X).
For each channel, for a pair of average waveforms x and y, we first scale x by
to minimize the sum of squared differences between x and y
![]() | (7) |
(x, y). We then compute two different distance measures d1 and d2. d1 is a normalized Euclidean distance between the scaled waveforms
![]() | (8) |
![]() | (9) |
![]() | (10) |
![]() | (11) |
![]() | (12) |
![]() | (13) |
![]() | (14) |
of the cumulative density function F can be obtained from the observed pairwise distances in the early unstable data. By setting the threshold at t =
–1(
) the probability for a type I error for each pair of waveform clusters with distance less than t is less than
. Because, for each cluster, the distances to all k clusters on the following day are computed, the threshold has to be adjusted for multiple comparisons. Thus for each neuron, the threshold is given by
![]() | (15) |
|
Orientation tuning functions were fitted using a least-squares algorithm and modified von Mises circular distribution functions given by
![]() | (16) |
is the neuron's preferred orientation, b determines the width of the orientation tuning function (the larger, the more sharply tuned; flat if zero), and f0 and a determine the peak firing rate and the offset above zero. Normalized Euclidean distance
A normalized Euclidean distance was used to compute the similarity between the orientation (or direction of motion) tuning functions of neurons. The normalized Euclidean distance between x and y is defined as
![]() | (17) |
|
|
RESULTS |
|---|
|
We recorded the spiking activity of single neurons from the primary visual cortex (V1) of the awake, behaving macaque using multiple chronically implanted tetrodes. To evaluate single-unit stability we did not adjust the location of the tetrodes for many months. We analyzed data recorded during this time where we ran the same experiments to map orientation tuning functions in all recording sessions. Multiple single neurons could be isolated from five of seven tetrodes in most recording sessions. Excluding two tetrodes that yielded only eight single units in 22 sessions, the mean number of well-isolated single units per session was 2.49/tetrode (range = 1.32–4.91; SD = 1.46). Our method for sorting single-unit spiking activity is based on fitting a MoG model in a 12-dimensional feature space using the SMEM algorithm (Ueda et al. 2000
) (see METHODS for details). We demonstrate the quality of single-unit isolation using examples from recording sessions during which more than four units were isolated from a tetrode. The scatterplots in Fig. 1A show peak-to-peak amplitudes of waveforms from all pairs of channels. Multiple clusters are clearly evident representing the spikes of putative different single units. The nine single units isolated by our clustering algorithm are plotted in different colors. In Fig. 1B we plot the waveforms projected onto the first principal component axis of each channel. In Fig. 1C we plot individual waveforms and the average waveform of the identified single units for each tetrode channel. The autocorrelograms of the units are shown on the right panel of Fig. 1C. The quality of single-unit isolation for all nine identified neurons recorded from this tetrode is demonstrated in Fig. 2. In Fig. 2A we demonstrate the single-unit isolation quality of the nine identified clusters from all other events not belonging to each one of these clusters. For each cluster, we plot the cumulative density of distances to the cluster center. The separation between pairs of clusters is illustrated by projecting onto the axis maximizing discriminability (see Fig. 2B captions for details). Out of these nine single units identified, seven exceeded our single-unit isolation criterion and are marked with red frames (Figs. 1C and 2A). Only neurons satisfying our single-unit isolation criterion are used for all analyses in this study (see METHODS for details). Additional examples of such well-isolated single-unit clusters are plotted in color in Fig. 3.
|
During periods when we did not adjust the location of the tetrodes, we observed that the spike waveforms and orientation tuning functions of some neurons recorded from the same tetrode were very similar across consecutive days (Fig. 4). This observation could be the result of recording stability and indicative of the fact that we were recording from the same neurons across days. To quantitatively determine single-unit stability across multiple days, the distance (or similarity) between all pairs of waveforms for two consecutive days were computed for each tetrode (distribution of pairwise waveform distances). Accordingly, two waveform clusters recorded across two days generated by the same neuron should have a small distance. However, a small distance does not necessarily imply that the two waveform clusters represent the same neuron. First, it is not a priori clear how similar different neurons in the same area are. This also strongly depends on the area from which one is recording. Second, although the amplitude ratio across the four tetrode channels uniquely depends on the location of the spike-generating neuron, this location is always relative to the tetrode. Thus on tetrode movement, in principle, two different neurons could generate highly similar waveforms on the same tetrode on consecutive days. It is important to point out that, for this reason, no measure of similarity—neither a distance metric nor a statistical test—can answer the question of recording stability solely based on two sets of waveforms recorded on consecutive days. Therefore first one has to empirically determine a null distribution of pairwise distances using data from a period where one is confident that there was no recording stability (e.g., because the tetrodes were adjusted). When tetrodes are not adjusted across days, then a significant deviation from this null distribution toward small distances would be an indication of recording stability.
|
The average waveforms of two neurons were compared using two different distance measures. For each channel of a tetrode separately, the average waveforms of the two neurons to be compared were scaled such that the sum of squared differences between them was minimized. The four scaling factors derived from this process were combined to form a distance parameter plotted on the ordinate of Fig. 5A. The sum of the normalized Euclidean distances between the scaled average waveforms for the four channels was plotted on the abscissa (see METHODS for details).
Figure 5B shows the distribution of pairwise distances between average waveforms recorded from the same tetrode across consecutive recording sessions when the tetrodes were not adjusted (experimental distribution). The experimental distribution is different from the null distribution in that an additional cluster considerably shifted toward the origin is clearly evident in the experimental distribution, indicating recording stability. To create the experimental distribution, data from a total of 22 recording sessions in two periods were used. Period A consisted of 13 recording sessions over the course of 23 days starting about 15 mo after implantation. The minimum temporal distance between two consecutive sessions was 1 day, the maximum 6 days. During this period, 13.2 ± 2.7 (mean ± SD) well-isolated single units could be recorded per session (172 in total). Period B consisted of 9 recording sessions on consecutive days starting about 6 mo after implantation. During this period, 12.2 ± 3.3 (mean ± SD) well-isolated single units could be recorded per session (110 in total).
Determining stable neurons can now be seen as a form of statistical hypothesis testing, where the null distribution defines the probability that two different neurons have a small waveform distance by chance. Because the null hypothesis of no stability should be rejected only for values close to the origin, defining a threshold in the two-dimensional distribution is nontrivial. Therefore we project both the null distribution and the experimental distribution into the one-dimensional subspace that optimally discriminates between the two modes in the experimental distribution (Fig. 5C). This subspace is found by fitting a mixture of two Gaussians model to the experimental distribution and computing the Fisher linear discriminant between the two Gaussians (see METHODS for details). It corresponds to taking a weighted average of the two distance measures. By computing the optimal subspace we empirically determine which weighting optimally separates stable neurons from different neurons. Note that the optimality criterion does not bias the statistical test in a way that could cause more false positives because testing is done using the marginalized null distribution, which has not been used to determine this subspace.
Using the preceding statistical framework (
= 0.05) we calculated for how many days the recordings remained sufficiently stable to reliably track the same neurons. This can be quantified in two slightly different ways as shown in Fig. 6. The first way we considered is to compute the average fraction of neurons recorded on a given day that can be recorded during the next n days. Given the average yield of single units per day, one can compute from these numbers how many neurons can be simultaneously recorded over a certain period of time. Although for period A we do not have data from every day, we consider it reasonable to assume that neurons that could be tracked across consecutive sessions spaced by >1 day would also have been tracked in between these two sessions if recordings had been made. Under this assumption, on average 61% (8.0 of 13.2) of the neurons recorded on any given day during period A and 41% (5.0 of 12.2) during period B could be tracked for
2 days; 52% (6.9 of 13.2; period A) and 26% (3.1 of 12.2; period B) could be tracked for
3 days; and 36% (4.7 of 13.2; period A) and 8% (1.0 of 12.2; period B) could be tracked for
7 days. The complete distributions for periods A and B are shown in Fig. 6A.
|
2 days, from 38 for
3 days, and from 10 for
1 wk. Validation of single-unit recording stability algorithm using orientation tuning functions
The data we used to compute the experimental distribution were collected during experiments where the orientation tuning functions of the neurons were measured. This allowed us to determine the stability of orientation tuning functions of neurons that we found to be stable across days. For each neuron, the orientation tuning functions were measured with either static (8 orientations) or moving gratings (16 directions of motion; see METHODS for details). For each visually responsive and orientation-tuned neuron recorded across two consecutive recording sessions the similarity between its orientation (or direction of motion) tuning functions was quantified using a normalized Euclidean distance that can assume values between 0 and 1 (a value of 0 indicates identical tuning functions). Not surprisingly, we find that neurons that were stable according to their waveform characteristics also had similar tuning functions across days (Fig. 7). The distribution of tuning function similarity of the same neurons was significantly different from the distribution of tuning function similarity between different neurons recorded on the same tetrode (Fig. 7, Wilcoxon rank-sum test, P < 10–36, n1 = 131, n2 = 632, rank sum = 20,029). Although the two distributions are very different, they do overlap to some extent; however, this is expected for several reasons.
|
2) If a tetrode moves a little too much, a neuron might still be recorded but, due to different channel-to-channel amplitude ratios, it might not be detected as stable. Given that this will be the same neuron, it will have the same tuning function.
3) As with any statistical method, there are false positives and misses, which to some extent create an overlap.
Taken together, it is not surprising that there are pairs of neurons classified as different that have small tuning differences. Finally, Fig. 8 shows additional examples of stable recordings across different numbers of days.
|
|
|
DISCUSSION |
|---|
|
Single-unit stability across days
We used a linear combination of two distance measures in spike waveform space to quantitatively evaluate single-unit stability across days. A null distribution of waveform distances for pairs of clusters that were generated by different neurons was constructed. This distribution was used in a statistical hypothesis testing procedure to determine which neurons were stable across days. Analysis of orientation tuning function differences revealed that neurons identified as stable have highly consistent orientation tuning functions across multiple days. This is to be expected because the animal was performing a passive fixation task, during which no perceptual learning is likely to occur.
Previously, recording stability across days has been claimed based on similarity between action potential waveforms (Greenberg and Wilson 2004
; Rousche and Normann 1998
; Schmitzer-Torbert and Redish 2004
; Taylor et al. 2002
; Wilson et al. 2003
) or waveform features (Nicolelis et al. 2003
). However, action potentials generated by different neurons might have highly similar waveforms, in particular in the case of single electrode recordings. The use of tetrodes partly resolves this problem by adding spatial information about the location of a neuron that is contained in the amplitude ratios across channels. This information, however, is always relative to the coordinate system of the tetrode. Thus at least in theory, two different neurons could produce exactly the same waveform on a tetrode on two different days. Therefore it is necessary to examine the distribution of pairwise waveform distances between different neurons recorded on the same tetrode. To the best of our knowledge, to date all approaches to evaluate recording stability across days share a common problem originating from the uncertainty described earlier—they all fail to give a reliable statistical confidence measure of single-unit stability. Other studies use functional stability of recorded neurons to verify recording stability (Greenberg and Wilson 2004
). Although this can be used to merely demonstrate recording stability, it cannot be used to study neurons whose tuning functions change during learning.
Recent work, somewhat similar in spirit to our method, tried to resolve these problems (Liu et al. 2006
). The distribution of feature differences for clusters generated by different neurons (equivalent to our null distribution) is approximated by comparing clusters across different electrodes. On the one hand, this guarantees that different neurons are being used. On the other hand, an estimate of the distribution might be systematically biased toward larger feature differences due to the different characteristics of different electrodes. In contrast, in our method clusters recorded from the same tetrode are used, avoiding this problem. Liu and colleagues (2006)
derived a threshold for feature differences based on the probability that two clusters were classified as one neuron, although they were generated by two different neurons. Unfortunately, this probability depends on the true percentage of neurons that are stable across days, a quantity that is unknown (in fact, it is part of what one tries to compute). To overcome this problem, these authors assumed this percentage to be
50%. Therefore the resulting error probability holds only if this assumption is true. It is important to realize that the validity of the assumption cannot be verified by the percentage of neurons that were found to be stable because this result was derived using the assumption and this would be a circular reasoning. Additionally, because the null distribution tends to overestimate distances between different neurons, the error probabilities need not hold, even if the assumed degree of a priori stability is true.
Although our approach to determine a null distribution based on early data right after implantation when the drive was adjusted on a regular basis does not have the disadvantages of comparing waveforms across different tetrodes (Liu et al. 2006
), it might not yield an unbiased estimate either. One could argue that tetrode properties could change with time, affecting the null distribution. Furthermore, it is possible that during the adjustment process different types of neurons are sampled from day to day, biasing the null distribution toward larger waveform differences. However, this does not seem to be the case for the results reported in this study because the large mode in the experimental distribution (which includes data recorded up to >1 year after implantation) looks highly similar to the null distribution. This is strong evidence that our estimate of the null distribution is relatively unbiased.
Future directions
One major improvement could be made to our method if the exact tetrode geometry inside the brain was known. If the locations of all four tetrode tips relative to each other are known, this can be used as a coordinate system and, using amplitude information, the location of each source can be triangulated. A method to measure tetrode geometry using a microscope and a camera as well as the derivation of the source localization using amplitude information is described in Chelaru and Jog (2005)
and Jog et al. (2002)
. Because tetrode movement has a much greater effect on neurons close to the tetrode than on distant ones, incorporating the spatial distance of two clusters across days could result in much better separability of clusters generated by the same versus different neurons.
|
|
GRANTS |
|---|
|
|
|
ACKNOWLEDGMENTS |
|---|
|
|
|
FOOTNOTES |
|---|
The costs of publication of this article were defrayed in part by the payment of page charges. The article must therefore be hereby marked "advertisement" in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.
Address for reprint requests and other correspondence: A. S. Tolias, Department of Neuroscience, Baylor College of Medicine, One Baylor Plaza, Suite S553, Houston, TX 77030 (E-mail: atolias{at}cns.bcm.edu)
|
|
REFERENCES |
|---|
|
Fahle M. Perceptual learning: specificity versus generalization. Curr Opin Neurobiol 15: 154–160, 2005.[CrossRef][Web of Science][Medline]
Fahle M, Edelman S. Long-term learning in vernier acuity: effects of stimulus orientation, range and of feedback. Vision Res 33: 397–412, 1993.[CrossRef][Web of Science][Medline]
Figueiredo M, Jain AK. Unsupervised learning of finite mixture models. IEEE Trans Pattern Anal Machine Intell 24: 381–396, 2002.[CrossRef]
Gilbert CD, Sigman M, Crist RE. The neural basis of perceptual learning. Neuron 31: 681–697, 2001.[CrossRef][Web of Science][Medline]
Görür D, Rasmussen CE, Tolias AS, Sinz F, Logothetis NK. Editors. Modelling Spikes with Mixtures of Factor Analysers. Berlin: Springer-Verlag, 2004, p. 391–398.
Gray CM, Maldonado PE, Wilson M, McNaughton B. Tetrodes markedly improve the reliability and yield of multiple single-unit isolation from multi-unit recordings in cat striate cortex. J Neurosci Methods 63: 43–54, 1995.[CrossRef][Web of Science][Medline]
Greenberg PA, Wilson FA. Functional stability of dorsolateral prefrontal neurons. J Neurophysiol 92: 1042–1055, 2004.
Harris KD, Henze DA, Csicsvari J, Hirase H, Buzsáki G. Accuracy of tetrode spike separation as determined by simultaneous intracellular and extracellular measurements. J Neurophysiol 84: 401–414, 2000.
Jog MS, Connolly CI, Kubota Y, Iyengar DR, Garrido L, Harlan R, Graybiel AM. Tetrode technology: advances in implantable hardware, neuroimaging, and data analysis techniques. J Neurosci Methods 117: 141–152, 2002.[CrossRef][Web of Science][Medline]
Judge SJ, Wurtz RH, Richmond BJ. Vision during saccadic eye movements. I. Visual interactions in striate cortex. J Neurophysiol 43: 1133–1155, 1980.
Liu X, McCreery DB, Bullara LA, Agnew WF. Evaluation of the stability of intracortical microelectrode arrays. IEEE Trans Neural Syst Rehabil Eng 14: 91–100, 2006.[CrossRef][Web of Science][Medline]
McLachlan J, Peel D. Finite Mixture Models. New York: Wiley, 2000.
Miyashita Y. Inferior temporal cortex: where visual perception meets memory. Annu Rev Neurosci 16: 245–263, 1993.[Web of Science][Medline]
Nicolelis MA, Dimitrov D, Carmena JM, Crist R, Lehew G, Kralik JD, Wise SP. Chronic, multisite, multielectrode recordings in macaque monkeys. Proc Natl Acad Sci USA 100: 11041–11046, 2003.
Poggio T, Bizzi E. Generalization in vision and motor control. Nature 431: 768–774, 2004.[CrossRef][Medline]
Pouzat C, Mazor O, Laurent G. Using noise signature to optimize spike-sorting and to assess neuronal classification quality. J Neurosci Methods 122: 43–57, 2002.[CrossRef][Web of Science][Medline]
Rousche PJ, Normann RA. Chronic recording capability of the Utah Intracortical Electrode Array in cat sensory cortex. J Neurosci Methods 82: 1–15, 1998.[CrossRef][Web of Science][Medline]
Sahani M. Latent Variable Models for Neural Data Analysis. Pasadena, CA: California Institute of Technology, 1999.
Schmitzer-Torbert N, Redish AD. Neuronal activity in the rodent dorsal striatum in sequential navigation: separation of spatial and reward responses on the multiple T task. J Neurophysiol 91: 2259–2272, 2004.
Squire LR, Stark CE, Clark RE. The medial temporal lobe. Annu Rev Neurosci 27: 279–306, 2004.[CrossRef][Web of Science][Medline]
Stopfer M, Bhagavan S, Smith BH, Laurent G. Impaired odour discrimination on desynchronization of odour-encoding neural assemblies. Nature 390: 70–74, 1997.[CrossRef][Medline]
Stopfer M, Laurent G. Short-term memory in olfactory network dynamics. Nature 402: 664–668, 1999.[CrossRef][Medline]
Taylor DM, Tillery SI, Schwartz AB. Direct cortical control of 3D neuroprosthetic devices. Science 296: 1829–1832, 2002.
Ueda N, Nakano R, Ghahramani Z, Hinton GE. SMEM algorithm for mixture models. Neural Comput 12: 2109–2128, 2000.[CrossRef][Web of Science][Medline]
Wilson FA, Ma YY, Greenberg PA, Ryou JW, Kim BH. A microelectrode drive for long term recording of neurons in freely moving and chaired monkeys. J Neurosci Methods 127: 49–61, 2003.[CrossRef][Web of Science][Medline]
Wilson MA, McNaughton BL. Dynamics of the hippocampal ensemble code for space. Science 261: 1055–1058, 1993.
This article has been cited by other articles:
![]() |
A. S. Dickey, A. Suminski, Y. Amit, and N. G. Hatsopoulos Single-Unit Stability Using Chronically Implanted Multielectrode Arrays J Neurophysiol, August 1, 2009; 102(2): 1331 - 1339. [Abstract] [Full Text] [PDF] |
||||
![]() |
A. Hernandez, V. Nacher, R. Luna, M. Alvarez, A. Zainos, S. Cordero, L. Camarillo, Y. Vazquez, L. Lemus, and R. Romo Procedure for recording the simultaneous activity of single neurons distributed across cortical areas during sensory discrimination PNAS, October 28, 2008; 105(43): 16785 - 16790. [Abstract] [Full Text] [PDF] |
||||
![]() |
M. J. Rasch, A. Gretton, Y. Murayama, W. Maass, and N. K. Logothetis Inferring Spike Trains From Local Field Potentials J Neurophysiol, March 1, 2008; 99(3): 1461 - 1476. [Abstract] [Full Text] [PDF] |
||||
| ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
| HOME | HELP | FEEDBACK | SUBSCRIPTIONS | ARCHIVE | SEARCH | TABLE OF CONTENTS |
| Visit Other APS Journals Online |