March 24, 2014 15:4 1450009
International Journal of Neural Systems, Vol. 24, No. 4 (2014) 1450009 (10 pages) c World Scientific Publishing Company DOI: 10.1142/S0129065714500099
NEV2lkit: A NEW OPEN SOURCE TOOL FOR HANDLING NEURONAL EVENT FILES FROM MULTI-ELECTRODE RECORDINGS MARKUS BONGARD∗,‡ and DANIEL MICOL†,‡ Institute of Bioengineering, Universidad Miguel Hern´ andez Avda. Universidad s/n, Elche, 03202 Alicante, Spain ∗ [email protected] † [email protected] Int. J. Neur. Syst. 2014.24. Downloaded from www.worldscientific.com by MONASH UNIVERSITY on 02/01/15. For personal use only.
§ ´ EDUARDO FERNANDEZ Institute of Bioengineering Universidad Miguel Hern´ andez and CIBER BBN Avda. Universidad s/n, Elche, 03202 Alicante, Spain [email protected]
Accepted 6 November 2013 Published Online 24 December 2013 The analysis and discrimination of action potentials, or “spikes”, is a central issue to systems neuroscience research. Here we introduce a free open source software for the analysis and discrimination of neural spikes based on principal component analysis and different clustering algorithms. The main objective is to supply a friendly user interface that links the experimental data to a basic set of routines for analysis, visualization and classification of spikes in a consistent framework. The tool has been tested on artificial data sets, on multi-electrode extracellular recordings from ganglion cell populations in isolated superfused mouse, rabbit and turtle retinas, and on electrophysiological recordings from mouse visual cortex. Our results show that NEV2lkit is very reliable and able to satisfy the experimental demands in terms of accuracy, efficiency and consistency across experiments. It performs fast unit sorting in single or multiple experiments and allows the extraction of spikes from over large time intervals in continuously recorded data streams. The tool is implemented in C++ and runs cross-platform on Linux, OS X and Windows systems. To facilitate the adaptation and extension as well as the addition of new routines, tools and algorithms for data analysis, the source code, binary distributions for different operating systems and documentation are all freely available at http://nev2lkit.sourceforge.net. Keywords: Extracellular recording; multi-unit recording; spike detection; spike sorting; unsupervised classification; electrode array; neural prosthesis.
1. Introduction Extracellular recordings of neural activity are the most practical choice in experiments that intend to investigate coding properties in neuronal populations1–7 or to study neural ensemble responses in humans and behaving animals.8–12 These techniques
‡ §
have allowed to get insights into how information is processed in large populations of neurons, and at the same time have provided new challenges in the analysis of neuronal data, particularly when several hundred neurons are recorded in a single experiment.13 However, while the use of extracellular electrode
These authors contributed equally to this work. Corresponding author. 1450009-1
March 24, 2014 15:4 1450009
Int. J. Neur. Syst. 2014.24. Downloaded from www.worldscientific.com by MONASH UNIVERSITY on 02/01/15. For personal use only.
M. Bongard, D. Micol & E. Fern´ andez
arrays allows imaging of multi-neuronal responses with unprecedented spatial and temporal resolution, the development of tools used to analyze this multi-neuronal activity is generally lagging behind the development of the tools used to acquire this data. Furthermore, although it is widely accepted that our understanding of neural coding, information transmission and brain processes would be aided if such data could be shared among individual laboratories,14 the data exchange between research groups using different multi-electrode acquisition systems is hindered by commercial constraints such as exclusive file structures, high priced licenses and hard policies on intellectual rights.15,16 As a result, most commercial software tend to have only limited possibilities to incorporate new tools or to modify existing ones and share them with colleagues. Multi-electrode recordings register the activity of many adjacent single neurons which have to be resolved into independent action potentials. This is possible because voltage variations due to a neuron depend on the inverse of its distance to the recording site.17–22 This implies that action potentials generated by a neuron located near the recording site are represented by higher voltage changes than spikes from neurons farther away. Furthermore, it is assumed that all waveforms from a single neuron are remarkably similar in overall shape, permitting them to be distinguished from noise. Likewise a classical approach to spike sorting includes at least three phases. First, the putative action potentials (spikes) must be detected and recorded using appropriate hardware and software.23 The second step is the isolation of the spikes and the final step involves the classification of the recorded spikes according to a set of features extracted from each individual signal.24–28 In spite of the complexities of acquiring and storing single unit responses from large numbers of neurons, as ensemble recordings gain in popularity, spike sorting becomes a limiting step in the analysis. Therefore we introduce a new free open source software, NEV2lkit (pronounced: nevtoolkit), to help researchers in this process. The main objective is to supply a friendly user interface that links the experimental data to a basic set of routines for analysis, visualization and classification of spikes in a consistent framework.
Most of the current literature on spike sorting considers automatic and non-parametric methods29–31 but sometimes this approach introduces errors in the classification. Thus the selection of criteria for spike sorting using NEV2lkit can be performed automatically or combining semiautomatic and manual approaches. Besides this semiautomatic spike sorting, the clustering algorithms are also aided by a real-time visual representation of the produced results (see Fig. 2). NEV2lkit is implemented in C++ and runs crossplatform on Linux, OS X and Windows systems. It works with single files or with folders containing multi-electrode data files in a variety of formats, such as ASCII-based text, LabView-formats, and Blackrock Microsystem’s Neural Event (NEV) files. Furthermore, it is able to extract spike events with a time base of up to 4 ms from continuous data recordings (i.e. from Multichannel Systems), and sort such data sets. Sorted spikes can be exported in several ASCII based or NEV formats for further usage with other software. To facilitate the adaptation and extension as well as the addition of new routines, tools and algorithms for data analysis, the source code, binary distributions for different operating systems and documentation are all freely available at http://nev2lkit.sourceforge.net. 2.
Methods
2.1. Experimental data Extracellular recordings were obtained from ganglion cell populations in isolated superfused rabbit, mouse and turtle retinas using an array of 100, 1.5 mm long electrodes as reported previously.2,5 Furthermore we also used multi-electrode recordings that were obtained from mouse visual cortex (see Fig. 1). Depending on the experiment, the inter-electrode distance of the multi-electrode array used was either 200 or 400 µm. The signals from the electrodes were amplified, passed through a bandwidth filter (300–2000 Hz) and digitized with a resolution of 16 bits at a sampling rate of 32 or 40 kHz. Neural spike events were detected by comparing the instantaneous electrode signal to level thresholds set for each data channel using standard procedures described elsewhere.16,30,31 Thus, when a supra-threshold event occurs, the signal window surrounding the event is
1450009-2
March 24, 2014 15:4 1450009
Int. J. Neur. Syst. 2014.24. Downloaded from www.worldscientific.com by MONASH UNIVERSITY on 02/01/15. For personal use only.
NEV2lkit: A New Open Source Tool for Handling Neuronal Event Files
Fig. 1. Examples of artificial data and raw recordings. A trace was plotted each time the voltage crossed the lower and/or the upper thresholds. For the artificial data (upper panel), additional data sets with 40% (a) and 80% (b) added Gaussian noise are also showed (for better visual discrimination of the signals, just a single unit is plotted in the figure). Only for the data sets obtained from brain slice recording the spikes were peak aligned.
time-stamped and stored for later, offline analysis. All the selected channels of data as well as the state of the stimulus were digitized with a commercial multiplexed A/D board data acquisition system (Bionic Technologies, Inc) and stored digitally. In the case of the data files recorded from mouse visual cortex, spikes were aligned with respect to the peak above threshold before the unit sorting processing. Additionally we also used artificial data sets that were recorded from a commercially available 100 channel spike simulator (Bionics Technologies Inc) using the same acquisition hardware. These artificial data (see Fig. 1) consisted of sequentially generated spikes with five different amplitudes on each of the 100 recording channels, followed by a burst sequence on all channels. Furthermore, by adding 40% and 80% of Gaussian noise using the “addnoise” program from the TISEAN package,32 we generated two additional data sets that were used in the analysis (Fig. 1, artificial A, B).
All data sets were chosen because they contained a large number of spikes on many electrodes and represented typical experimental results.
2.1.1. Operational flow of NEV2lkit NEV2lkit is structured in two main modules, each designed to perform specific tasks of the neuronal event extraction and unit sorting process. The Input/Output module offers the opening, extraction and storage of spike event data. The spike sorting module incorporates two main steps. The first step is the conversion of the simultaneously recorded spikes into a clustering space constructed with the features of the signal shape extracted from applying the PCA. The second step attempts to separate this data into functionally meaningful “unit” clusters. After the program is launched, the graphical user interface of NEV2lkit provides easy access to all functions (see Fig. 2). The output routines allow
1450009-3
March 24, 2014 15:4 1450009
Int. J. Neur. Syst. 2014.24. Downloaded from www.worldscientific.com by MONASH UNIVERSITY on 02/01/15. For personal use only.
M. Bongard, D. Micol & E. Fern´ andez
(a)
(b)
(c)
Fig. 2. (Color online) Example of clustering results for a representative electrode in a multielectrode recording of rabbit ganglion cells to a moving bar. (a) Raw data before spike sorting. A total of 4622 events were recorded. The time scale is 1.6 ms. (b) Application of PCA to the data (spike vector) using the first three principal components for spike representation and a penalty function of 1, that penalizes overfitting by discouraging models with a large number of parameters. The sorted waveforms are shown in different colors. (c) The spike shapes corresponding to four putative neurons (same colors than b: red, green, blue and cyan). Numbers in parentheses show the total number of spikes. More distant neurons contribute to the background noise (gray).
to save data in NEV format or to export the time stamps of the spike event data as ASCII files with or without the voltage data of each individual spike. The sorting module operates on single or multiple files, allowing to automatically assign putative units over several experiments. Both modules use a predefined default or user defined parameters in the NEV2lkit GUI for their processing but run otherwise
in unsupervised mode. An exhaustive log file is saved to disk and provides the researcher all the information needed to reproduce the spike sorting process. A key feature of NEV2lkit is that it is able to process many files simultaneously. Therefore the assignments of units or clusters are done in many experiments at the same time. As a result the same
1450009-4
March 24, 2014 15:4 1450009
NEV2lkit: A New Open Source Tool for Handling Neuronal Event Files
unit on a given electrode has exactly the same unit label in all selected files. The number of experiments which can be processed simultaneously with NEV2lkit is only limited by the available RAM of the computer on which the program runs. In most cases, the software allows to process experiments with hundreds or even thousands of channels without any problem, therefore in order to speed up the sorting process there is no need to buffer to disk. Nevertheless, if needed, the source code can be easily modified to use the hard disk for temporary storage.
Int. J. Neur. Syst. 2014.24. Downloaded from www.worldscientific.com by MONASH UNIVERSITY on 02/01/15. For personal use only.
2.2. Spike sorting Spike sorting is an important step in the analysis of multi-electrode data. The main purpose is to distinguish between multiple cells recorded together.33–35 In this context, it should be taken into account that the extent to which the responses of a single neuron can be characterized depends on the ability to identify the action potentials originating from a single cell, and to discriminate these action potentials from other sources of electrical activity.36 This process, and therefore all further classification and analysis, depends on the noise present, the action potential amplitude, the conductive properties of the extracellular medium, the micro-electrode configuration, the localization with respect to the cell, and the type and quality of the micro-electrodes.37,38
2.2.1. Principal component analysis We apply a Principal Component Analysis (PCA) to reduce the dimensionality representation of the recorded waveforms in a p-dimensional space. Whereas detailed discussion of PCA is available elsewhere,29,39–41 here we briefly describe the basis of the method that we implemented for those unfamiliar with the technique. In the learning phase a set of spikes are accumulated, and from these the correlation matrix and its eigenvectors are computed. Then, a multidimensional matched filter using principal component vectors is defined and used in a real-time digital filter. The outputs of these filters can be used to classify the recorded spikes. Additionally, we obtain an ordered set of orthogonal basis vectors by calculating the a
eigenvectors of the variance-covariance matrix of the loaded data sets. Each principal component vector represents the most salient features of a given waveform and the eigenvectors with the largest eigenvalues correspond to the dimensions that have the strongest correlation in the dataset. Thus, the first eigenvector represents the largest correlation across the waveforms, the second eigenvector captures as much as possible of the waveforms differences not displayed by the first eigenvector, and so on. As a result, each spike is represented as a weighted sum of eigenvectors, corresponding to the largest eigenvalues of the data sets matrix. Since it is known that the signal-to-noise ratio can be optimized using the first M components of the eigenvectors in which the signal energy exceeds the noisy energy, we calculated the components that were significantly above background noise. We found that the first two components were often sufficient and that using additional components beyond the third one would not improve classification accuracy. These results agreed with previous studies reporting that one or two of the principal component vectors are sufficient to represent 80–90% of the spike waveform.42,43 Therefore, we only used the first three principal component vectors, which in our cases represent approximately 78% of the variation in the analyzed data, reducing the p-dimensional space to three dimensions. We performed the PCA on all multi-electrode channels and therefore the directions of maximum variation on all individual channels can be described simultaneously. Furthermore, we used different matrices for calculating the PCA vectors: the eigenvectors of the correlation one, and of the variance-covariance. 2.2.2. Clustering After the PCA has been calculated, each registered spike is represented as a point in a three-dimensional system spanned by the first three principal components. We used the Klustakwika implementation of the CEM algorithm of Celeux and Govaert to calculate the average separation of the spikes in the p-space. This algorithm computes all different spike groups in the hyperplane, and also takes into account
K.D. Harris, http://klusta-team.github.io/klustakwik/. 1450009-5
March 24, 2014 15:4 1450009
M. Bongard, D. Micol & E. Fern´ andez
the spreading of the spikes within the spike groups. As a result, the algorithm delivers the maximal separation between object groups (in our case different waveforms) combined with minimal spreading within the object groups (the same spike waveform). Klustakwik also allows to use Bayesian information in the clustering process.44 This option can be used to quantify the goodness of the resulting classification which is a useful aid to investigate the isolations of spikes in different clusters.45
Int. J. Neur. Syst. 2014.24. Downloaded from www.worldscientific.com by MONASH UNIVERSITY on 02/01/15. For personal use only.
2.3. Spike detection and extraction We performed spike detection in both supervised and unsupervised modes. In the supervised mode a lower and a higher detection threshold can be individually adjusted for each open data stream. For the extraction, a window time base for the detected events can be defined by the user (up to 4 ms). Extracted events are then displayed in this defined extraction window within the program. In automatic or unsupervised mode, NEV2lkit opens a suitable file, calculates the square root of the time-averaged (mean), squared value of the signal amplitude — the so-called root mean square or RMSE — of the individual signal streams, sets the upper and lower thresholds as RMSE×3.0 and extracts and displays the spikes over the defined time window. 3.
therefore be exactly evaluated. Each recording lasted 60 s and consisted of a sequence of single artificial spikes, each characterized by a different amplitude and time course of the signal, with an inter-spike interval of one second. The sequence of single spikes was followed by a burst of 10 identical spikes. For the analysis we selected randomly 10 electrodes from 20 different recordings. Thresholds on the recording channels were individually set allowing the extraction of different numbers of units on every single electrode. For further performance evaluation additional test sets were generated by adding 40% and 80% of Gaussian noise to this data (Fig. 1, upper panel). Results of the performance of the algorithms implemented in NEV2lkit for sorting the artificial data are summarized in Table 1, where the percentage of spikes used in the classification is expressed as the percentage of the total number of spikes. Erroneously classified spikes are given as a percentage of the total data set. The clustering algorithms were able to automatically remove outliers and artifacts from the data set and approximately 97% of all units were assigned to clusters. The best performance was achieved by choosing the maximization of the correlations of the basis vectors for the PCA. When the maximization of the variance-covariance was chosen to reduce the dimensionality representation of the recorded waveforms, more spikes were assigned to wrong units.
Results
The basic problems in spike sorting are illustrated in Fig. 2, which presents some waveforms recorded by a electrode in a typical multi-electrode experiment. A close inspection shows that the recorded spike waveforms usually contain a mixture of spikes from one or more cells and a significant amount of background activity from unknown sources. In a first step the quality of isolation can be tested by visually inspecting the superimposed spikes on each electrode, but even for relatively well-isolated spikes the recordings include some overlapping of spike events. Hence, the key problem is to correctly classify all the discriminable neurons.
Table 1. Summary of the clustering results for the artificial data. Method for PCA VarianceSum of Correlation covariance squares maximization maximization minimization
% of spikes used % error % of spikes used % error
96.7 0.28 96.7 0.28
3.1. Artificial data The advantage of simulated data is that all signal characteristics, such as the number of units within each recorded file, is known in advance and can
% of spikes used % error
1450009-6
96.7 0.56
0% noise 96.7 5.91 40% noise 96.7 5.91 80% noise 96.7 10.42
96.2 6.93 96.2 9.37 95.1 10.48
March 24, 2014 15:4 1450009
NEV2lkit: A New Open Source Tool for Handling Neuronal Event Files
These results were very reliable, even in the presence of 80% Gaussian noise. Consequently it seems that choosing the “right” calculation method improves significantly the accuracy of the clustering process.
Table 2. Clustering results for the experimental extracellular recordings. Method for PCA VarianceSum of Correlation covariance squares maximization maximization minimization
Int. J. Neur. Syst. 2014.24. Downloaded from www.worldscientific.com by MONASH UNIVERSITY on 02/01/15. For personal use only.
3.2. Extracellular recorded spikes While for experimental recordings the “real” number of unit clusters is not known, the quality of unit isolation can be judged, on a first approach, by looking at the superimposed spikes on every cluster. This delivers a first approximation of how many clusters exist within the data set recorded by an electrode. Figure 2(a) shows a typical data set from rabbit retina and the resulting classification into isolated units by using NEV2lkit (Fig. 2(b)). As we can see, even singular outliers are correctly sorted into an individual cluster (light-gray trace in Fig. 2(c)). The inter-spike interval histograms (ISI) and autocorrelation function histograms further document the high quality of the unit sorting carried out by NEV2lkit and show that there is no event assigned to clusters with an inter-event distance less than the refractory period. The overall qualification of the unit sorting performance for all multi-electrode recording data sets using different methods for the PCA is summarized in Table 2. In an attempt to quantify the goodness of the results, we calculated the ratio between the number of expected classes (from visual inspection of the recordings) and the classes isolated by NEV2lkit. Values higher than 1 mean that the unit sorting process is assigning less classes than the number suggested by visual inspection; 1 means a number of classes equal to the one assumed by visual inspection; and values lower than 1 represent a higher number of classes than expected. As it can be seen the results depend strongly on the PCA method. In order get more insight about the usefulness of every procedure, we also quantify the percentage of spikes classified incorrectly as noise. Due to the experimental setups used, noise from external sources were a rare event in our recordings and its percentage was less than 0.4% of the total signals. Therefore, the increase in the number of spikes assigned to noisy units is somehow related with wrongly assigned spikes. For each data set the noise units were determined and pooled, as well as the
Unit ratio % noisy units
Mouse brain 1 1 1.76 0.44
2 1.76
Unit ratio % noisy units
Rabbit retina 1.5 1 2.54 3.04
1 0.10
Unit ratio % noisy units
Turtle retina 1 0.75 3.61 0.45
1 0.02
Unit ratio % noisy units
Mouse retina 3 3 1.12 3.16
1 1.12
classification error expressed as the percentage of spikes incorrectly assigned to noise units. Taken together these results suggest that it is possible to find an optimal method for the classification of every data set. For example, the artificial data was best sorted into units using the maximization of the correlations as basis vectors for the PCA analysis, whereas for retinal ganglion cells recordings the best performance was achieved by choosing the maximization of the variance-covariances. Because all data were recorded with the same type of multielectrode array and electrophysiological equipment, one can assume that differences in the signal characteristics are indeed depending on the signal origin. 4.
Discussion
One of the most promising advances in neuroscience research is the development of multi-electrode arrays that allow multiple site recordings in various neural systems.7,8,46,47 A major challenge in this context is to acquire meaningful data from a large number of channels but a central problem is how many neurons are present in a given recording. Therefore, finding the spiking activity in electrophysiological recordings is seen as a first step in the decoding of neural activity. However there are usually a wide variability in the results of manual and
1450009-7
March 24, 2014 15:4 1450009
Int. J. Neur. Syst. 2014.24. Downloaded from www.worldscientific.com by MONASH UNIVERSITY on 02/01/15. For personal use only.
M. Bongard, D. Micol & E. Fern´ andez
even automatic spike sorting across a given data set among different subjects and laboratories which suggest the importance of good, widely available, automated spike sorting methods and tools. To facilitate this task we have developed an open source software which has proved to be useful for the classification and management of neural ensemble data. The main idea is to supply a graphical user interface that links the experimental data to a basic set of routines whose purpose is to distinguish between multiple cells recorded together and to assign detected spiking activity to specific neurons. An important issue in this field is to distinguish between single and multi-unit activity. Our results show that the procedures implemented in NEV2lkit using a combination of adjustable PCA- and CEMbased clustering can precisely detect single neurons even in the presence of large amounts of noise. Hence NEV2lkit offers a fast, reliable and accurate implementation for spike sorting. While it is known that template matching can deliver a more accurate classification of the signals,24,25 our results show that the PCA- and CEM-based approach is sufficiently accurate. Furthermore, one advantage of NEV2lkit lies in the fact that it operates fast and in an unsupervised mode over a large number of neuronal recording files so that the user only has to define the minimal and maximal class assignment limits. Thus, by using NEV2lkit it is possible, depending on the overall data size, to identify units across multiple experiments within seconds while the program runs with minimal required human intervention. Our results reveal that CEM-based clustering of the principal components calculated using the maximization of the correlations as basis vectors for the PCA works very well for synthetic data. Interestingly, it turns out that for the real extracellular recordings used in this analysis, the basis vector construction method was not delivering the best performance in the unit sorting process. As a result, the researcher has to define and choose the best method depending on the data and its particular features. NEV2lkit addresses this requirement by offering a selection of different basis vector methods to be used as base for the sorting process. This allows handling unit sorting problems which are presumably partly
depending on the signal origin. Since the physical multi-electrode array characteristics were similar in all the recordings used, our results suggest that there is an optimal PCA approach for spike sorting and unit assignment depending on every data set. Restrictions for the application of NEV2lkit are related to the assumption that extracellular spike shapes are always stationary. This expectation could be realized in our analyzed data sets, because our experimental setup delivers almost stationary recorded spike waveforms. However, there could be additional problems that limit the robustness of the unit sorting carried out by these procedures. For example, in the analyzed experiments the level of background noise remained relatively constant and electrode drift was not observed or was minimal. If these conditions are given, NEV2lkit suits apparently well for fast off-line unit sorting in multi-electrode recordings from different kinds of sources. We are aware that some functions currently in use, as well as advanced functions used in various labs working with multi-electrode data, are not included. We hope that the facilities to extend the code and add more specific algorithms and routines can compensate the former and encourages users to contribute providing ready-made functions useful to the scientific community. Thus, the modular program architecture allows an easy extension of its functionality, such as extending the number of supported data file formats or the integration of external libraries for spike detection. In this context there are already available versions of NEV2lkit that provide extended functionality to the feature extraction procedureb such as variable length data window, variability of the number of principal components and calculation and display of variance.48 Finally, NEV2lkit is free software distributed under the GNU General Public License.c This license grants it users legal permission to copy, distribute and/or modify the software. To this end, the source code is distributed alongside of executable binaries for Linux, OS X and Windows. The only restriction imposed on the user is that should they decide to redistribute the software in its original or modified form, they must in turn grant their users the same rights. We hope to attract an active community of
b http://neurobot.bio.auth.gr/nev2lkit. c
http://www.gnu.org/licenses/gpl.html. 1450009-8
March 24, 2014 15:4 1450009
NEV2lkit: A New Open Source Tool for Handling Neuronal Event Files
users who could benefit from and contribute to this project. Acknowledgments
10.
Int. J. Neur. Syst. 2014.24. Downloaded from www.worldscientific.com by MONASH UNIVERSITY on 02/01/15. For personal use only.
We would like to thank J. V. Segura for his valuable comments. In addition, we wish to thank M. Gerschner, L. van Ahrens and K. Gansel for providing some of the data sets used in this work, and V. Medran and S. Andreu for technical assistance. This work has been partially supported by Spanish Grants MAT2012-39290-C02-01 and IPT-20120574-300000, by ONCE (National Organization of the Spanish Blind) and by the Research Chair in Retinitis Pigmentosa Bidons Egara.
11.
12.
13.
References 1. M. Greschner, M. Bongard, P. Rujan and J. Ammermuller, Retinal ganglion cell synchronization by fixational eye movements improves feature estimation, Nat. Neurosci. 5 (2002) 341–347. 2. E. Fernandez, J. M. Ferrandez, J. Ammermller and R. A. Normann, Population coding in spike trains of simultaneously recorded retinal ganglion cellsx, Brain Res. 887 (2002) 222–229. 3. M. J. Schnitzer and M. Meister, Multineuronal firing patterns in the signal from eye to brain, Neuron 37 (2003) 499–511. 4. S. Takahashi, Y. Anzai and Y. Sakurai, Automatic sorting for multi-neuronal activity recorded with tetrodes in the presence of overlapping spikes, J. Neurophysiol. 89 (2003) 2245–2258. 5. G. J. Ortega, M. Bongard, E. Louis and E. Fernandez, Conditioned spikes: A simple and fast method to represent rates and temporal patterns in multielectrode recordings, J. Neurosci. Methods 133 (2004) 135–141. 6. N. V. Manyakov and M. M. Van Hulle, Decoding grating orientation from microelectrode array recordings in monkey cortical area V4, Int. J. Neural Syst. 20(2) (2010) 95108. 7. F. Franke, D. J¨ ackel, J. Dragas, J. M¨ uller, M. Radivojevic, D. Bakkum and A. Hierlemann, Highdensity microelectrode array recordings and realtime spike sorting for closed-loop experiments: An emerging technology to study neural plasticity, Front. Neural Circuits 6(105) (2012) doi: 10.3389/fncir.2012.00105. 8. G. Van Dijck, K. Seidl, O. Paul, P. Ruther, M. M. Van Hulle and R. Maex, Enhancing the yield of high-density electrode arrays through automated electrode selection, Int. J. Neural Syst. 22(1) (2012) 119. 9. M. A. Lebedev, J. M. Carmena, J. E. O’Doherty, M. Zacksenhouse, C. S. Henriquez, J. C. Principe
14. 15.
16.
17.
18.
19.
20. 21.
22.
1450009-9
and M. A. Nicolelis, Cortical ensemble adaptation to represent velocity of an artificial actuator controlled by a brain-machine interface, J. Neurosci. 25 (2005) 4681–4693. R. Gutierrez, J. M. Carmena, M. A. Nicolelis and S. A. Simon, Orbitofrontal ensemble activity monitors licking and distinguishes among natural rewards, J. Neurophysiol. 95 (2006) 119–133. L. R. Hochberg, M. D. Serruya, G. M. Friehs, J. A. Mukand, M. Saleh, A. H. Caplan, A. Branner, D. Chen, R. D. Penn and J. P. Donoghue, Neuronal ensemble control of prosthetic devices by a human with tetraplegia, Nature 442 (2006) 164–171. J. P. Donoghue, A. Nurmikko, M. Black and L. R. Hochberg, Assistive technology and robotic control using motor cortex ensemble-based neural interface systems in humans with tetraplegia, J. Physiol. 579 (2007) 603–611. N. G. Hatsopoulose and J. P. Donoghue, The science of neural interface systems, Annu. Rev. Neurosci. 32 (2007) 249–266. D. Gardner, Neurodatabase.org: Networking the microelectrode, Nat. Neurosci. 7 (2004) 486–487. U. Egert, T. Knott, C. Schwarz, M. Nawrot, A. Brandt, S. Rotter and M. Diesmann, MEA-Tools: An open source toolbox for the analysis of multielectrode data with MATLAB, J. Neurosci. Methods 117 (2002) 33–42. M. P. Bonomini, J. M. Ferrandez, J. A. Bolea and E. Fernandez, DATA-MEAns: An open source tool for the classification and management of neural ensemble recordings, J. Neurosci. Methods 148 (2005) 137–146. N. R. Luque, J. A. Garrido, J. Ralli, J. J. Laredo and E. Ros, From sensors to spikes: Evolving receptive fields to enhance sensori-motor information in a robot-arm, Int. J. Neural Syst. 22(4) (2012) 1250013. K. Ramanathan, N. Ning, D. Dhanasekar, L. Guoqi, S. Luping and P. Vadakkepat, Presynaptic learning and memory with a persistent firing neuron and a habituating synapse: A model of short term persistent habituation, Int. J. Neural Syst. 22(4) (2012) 1250015. W. K. Wong, Z. Wang and B. Zhen, Relationship between applicability of current-based synapses and uniformity of firing patterns, Int. J. Neural Syst. 22(4) (2012) 1250017. S. Ghosh-Dastidar and H. Adeli, Spiking neural networks, Int. J. Neural Syst. 19(4) (2009) 295–308. S. Ghosh-Dastidar and H. Adeli, A new supervised learning algorithm for multiple spiking neural networks with application in epilepsy and seizure detection, Neural Netw. 22(10) (2009) 1419–1431. S. Ghosh-Dastidar and H. Adeli, Improved spiking neural networks for EEG classification and epilepsy and seizure detection, Integr. Comput-Aided. Eng. 14(3) (2007) 187–212.
March 24, 2014 15:4 1450009
Int. J. Neur. Syst. 2014.24. Downloaded from www.worldscientific.com by MONASH UNIVERSITY on 02/01/15. For personal use only.
M. Bongard, D. Micol & E. Fern´ andez
23. K. Sameshima and M. Bach, Trends in multichannel neural ensemble recording instrumentation, in Methods for Neural Ensemble Recordings, ed. M. A. Nicolelis (CRC Press, Boca Raton, Florida, 1999), pp. 47–60. 24. R. A. Normann, D. J. Warren, J. Ammermuller, E. Fernandez and S. Guillory, High-resolution spatiotemporal mapping of visual pathways using multielectrode arrays, Vision Res. 41 (2001) 1261–1275. 25. S. Shoham, M. R. Fellows and R. A. Normann, Robust, automatic spike sorting using mixtures of multivariate t-distributions, J. Neurosci. Methods 127 (2003) 111–122. 26. U. Rutishauser, E. M. Schuman and A. N. Mamelak, Online detection and sorting of extracellularly recorded action potentials in human medial temporal lobe recordings, in vivo, J. Neurosci. Methods 154 (2006) 204–224. 27. S. Ghosh-Dastidar and H. Adeli, Wavelet-clusteringneural network model for freeway incident detection, Comput. Aided Civ. Inf. 18(5) (2003) 325–338. 28. M. Ahmadlou and H. Adeli, Enhanced probabilistic neural network with local decision circles: A robust classifier, Integr. Comput. Aid. E 17(3) (2010) 197–210. 29. E. Chan, V. Hok, A. Della-Chiesa, J. J. H. Miller, S. M. O’Mara and R. B. Reilly, Automated spike sorting algorithm based on Laplacian eigenmaps and k-means clustering, J. Neural Eng. 8(1) (2011) 016006. 30. Z. Tiganj and M. Mboup, A non-parametric method for automatic neural spikes clustering based on the non-uniform distribution of the data, J. Neural Eng. 8(6) (2011) 066014. 31. Y. Yuan, C. Yang and J. Si, The M-Sorter: An automatic and robust spike detection and classification system, J. Neurosci. Meth. 210 (2012) 281–290. 32. R. Hegger, H. Kantz and T. Schreiber, Practical implementation of nonlinear time series methods: The TISEAN package, Chaos 9 (1999) 413–435. 33. T. Takekawa, Y. Isoura and T. Fukai, Accurate spike sorting for multi-unit recordings, Eur. J. Neurosci. 31 (2010) 263272. 34. I. Cajigas, W. Q. Malik and E. N. Brown, nSTAT: Open-source neural spike train analysis toolbox for Matlab, J. Neurosci. 211 (2012) 45264. 35. C. Vargas-Irwin and J. P. Donoghue, Automated spike sorting using density grid contour clustering and subtractive waveform decomposition, J. Neurosci. Methods 164 (2007) 1–18. 36. N. Schmitzer-Torbert, J. Jackson, D. Henze, K. Harris and A. D. Redish, Quantitative measures of cluster quality for use in extracellular recordings, Neuroscience 131 (2005) 1–11.
37. H. Kaneko, H. Tamura and S. S. Suzuki, Tracking spike-amplitude changes to improve the quality of multineuronal data analysis, IEEE Trans. Biomed. Eng. 54 (2007) 262–272. 38. T. Takekawa, Y. Isoura and T. Fukai, Spike sorting of heterogeneous neuron types by multimodalityweighted PCA and explicit robust variational Bayes, Front. Neuroinform. 6(5) (2012) doi: 10.3389/fninf.2012.00005. 39. S. Ghosh-Dastidar, H. Adeli and N. Dadmehr, Principal component analysis-enhanced cosine radial basis function neural network for robust epilepsy and seizure detection, IEEE Trans. Biomed. Eng. 55(2) (2008) 512–518. 40. M. Al-Naser and U. Soderstrom, Reconstruction of occluded facial images using asymmetrical principal component analysis, Integr. Comput. Aided Eng. 19(3) (2012) 273–283. 41. A. Meraoumia, S. Chitroub and A. Bouridane, 2D and 3D palmprint information, PCA and HMM for an improved person recognition performance, Integr. Comput. Aided Eng. 20(3) (2013) 303–319. 42. R. Vollgraf, M. Munk and K. Obermayer, Optimal filtering for spike sorting of multi-site electrode recordings, Network 16 (2005) 85–113. 43. P. M. Horton, A. U. Nicol, K. M. Kendrick and J. F. Feng, Spike sorting based upon machine learning algorithms (SOMA), J. Neurosci. Methods 160 (2007) 52–68. 44. K. D. Harris, D. A. Henze, J. Csicsvari, H. Hirase and G. Buzsaki, Accuracy of tetrode spike separation as determined by simultaneous intracellular and extracellular measurements, J. Neurophysiol. 84 (2000) 401–414. 45. J. Wilda, Z. Prekopcsakc, T. Siegera, D. Novaka and R. Jech, Accuracy of tetrode spike separation as determined by simultaneous intracellular and extracellular measurements, J. Neurosci. Methods 203 (2012) 369–376. 46. M. R. Cohen and A. Kohn, Measuring and interpreting neuronal correlations, Nature 14(7) (2012) 811–819. 47. H. A. Wark, R. Sharma, K. S. Mathews, E. Fernandez, J. Yoo, B. Christensen, P. Tresco, L. Rieth, F. Solzbacher, R. A. Normann and P. Tathireddy, A new high-density (25 electodres/mm3) penetrating microelectrode array for recording and stimulating sub-millimeter neuroanatomical structures, J. Neural Eng. 10 (2013) 045003. 48. D. A. Adamos, E. K. Kosmidis and G. Theophildis, Performance evaluation of PCA-based spike sorting algorithms, Comput. Methods Programs Biomed. 91 (2008) 232–244.
1450009-10
International Journal of Neural Systems, Vol. 24, No. 4 (2014) 1450009 (10 pages) c World Scientific Publishing Company DOI: 10.1142/S0129065714500099
NEV2lkit: A NEW OPEN SOURCE TOOL FOR HANDLING NEURONAL EVENT FILES FROM MULTI-ELECTRODE RECORDINGS MARKUS BONGARD∗,‡ and DANIEL MICOL†,‡ Institute of Bioengineering, Universidad Miguel Hern´ andez Avda. Universidad s/n, Elche, 03202 Alicante, Spain ∗ [email protected] † [email protected] Int. J. Neur. Syst. 2014.24. Downloaded from www.worldscientific.com by MONASH UNIVERSITY on 02/01/15. For personal use only.
§ ´ EDUARDO FERNANDEZ Institute of Bioengineering Universidad Miguel Hern´ andez and CIBER BBN Avda. Universidad s/n, Elche, 03202 Alicante, Spain [email protected]
Accepted 6 November 2013 Published Online 24 December 2013 The analysis and discrimination of action potentials, or “spikes”, is a central issue to systems neuroscience research. Here we introduce a free open source software for the analysis and discrimination of neural spikes based on principal component analysis and different clustering algorithms. The main objective is to supply a friendly user interface that links the experimental data to a basic set of routines for analysis, visualization and classification of spikes in a consistent framework. The tool has been tested on artificial data sets, on multi-electrode extracellular recordings from ganglion cell populations in isolated superfused mouse, rabbit and turtle retinas, and on electrophysiological recordings from mouse visual cortex. Our results show that NEV2lkit is very reliable and able to satisfy the experimental demands in terms of accuracy, efficiency and consistency across experiments. It performs fast unit sorting in single or multiple experiments and allows the extraction of spikes from over large time intervals in continuously recorded data streams. The tool is implemented in C++ and runs cross-platform on Linux, OS X and Windows systems. To facilitate the adaptation and extension as well as the addition of new routines, tools and algorithms for data analysis, the source code, binary distributions for different operating systems and documentation are all freely available at http://nev2lkit.sourceforge.net. Keywords: Extracellular recording; multi-unit recording; spike detection; spike sorting; unsupervised classification; electrode array; neural prosthesis.
1. Introduction Extracellular recordings of neural activity are the most practical choice in experiments that intend to investigate coding properties in neuronal populations1–7 or to study neural ensemble responses in humans and behaving animals.8–12 These techniques
‡ §
have allowed to get insights into how information is processed in large populations of neurons, and at the same time have provided new challenges in the analysis of neuronal data, particularly when several hundred neurons are recorded in a single experiment.13 However, while the use of extracellular electrode
These authors contributed equally to this work. Corresponding author. 1450009-1
March 24, 2014 15:4 1450009
Int. J. Neur. Syst. 2014.24. Downloaded from www.worldscientific.com by MONASH UNIVERSITY on 02/01/15. For personal use only.
M. Bongard, D. Micol & E. Fern´ andez
arrays allows imaging of multi-neuronal responses with unprecedented spatial and temporal resolution, the development of tools used to analyze this multi-neuronal activity is generally lagging behind the development of the tools used to acquire this data. Furthermore, although it is widely accepted that our understanding of neural coding, information transmission and brain processes would be aided if such data could be shared among individual laboratories,14 the data exchange between research groups using different multi-electrode acquisition systems is hindered by commercial constraints such as exclusive file structures, high priced licenses and hard policies on intellectual rights.15,16 As a result, most commercial software tend to have only limited possibilities to incorporate new tools or to modify existing ones and share them with colleagues. Multi-electrode recordings register the activity of many adjacent single neurons which have to be resolved into independent action potentials. This is possible because voltage variations due to a neuron depend on the inverse of its distance to the recording site.17–22 This implies that action potentials generated by a neuron located near the recording site are represented by higher voltage changes than spikes from neurons farther away. Furthermore, it is assumed that all waveforms from a single neuron are remarkably similar in overall shape, permitting them to be distinguished from noise. Likewise a classical approach to spike sorting includes at least three phases. First, the putative action potentials (spikes) must be detected and recorded using appropriate hardware and software.23 The second step is the isolation of the spikes and the final step involves the classification of the recorded spikes according to a set of features extracted from each individual signal.24–28 In spite of the complexities of acquiring and storing single unit responses from large numbers of neurons, as ensemble recordings gain in popularity, spike sorting becomes a limiting step in the analysis. Therefore we introduce a new free open source software, NEV2lkit (pronounced: nevtoolkit), to help researchers in this process. The main objective is to supply a friendly user interface that links the experimental data to a basic set of routines for analysis, visualization and classification of spikes in a consistent framework.
Most of the current literature on spike sorting considers automatic and non-parametric methods29–31 but sometimes this approach introduces errors in the classification. Thus the selection of criteria for spike sorting using NEV2lkit can be performed automatically or combining semiautomatic and manual approaches. Besides this semiautomatic spike sorting, the clustering algorithms are also aided by a real-time visual representation of the produced results (see Fig. 2). NEV2lkit is implemented in C++ and runs crossplatform on Linux, OS X and Windows systems. It works with single files or with folders containing multi-electrode data files in a variety of formats, such as ASCII-based text, LabView-formats, and Blackrock Microsystem’s Neural Event (NEV) files. Furthermore, it is able to extract spike events with a time base of up to 4 ms from continuous data recordings (i.e. from Multichannel Systems), and sort such data sets. Sorted spikes can be exported in several ASCII based or NEV formats for further usage with other software. To facilitate the adaptation and extension as well as the addition of new routines, tools and algorithms for data analysis, the source code, binary distributions for different operating systems and documentation are all freely available at http://nev2lkit.sourceforge.net. 2.
Methods
2.1. Experimental data Extracellular recordings were obtained from ganglion cell populations in isolated superfused rabbit, mouse and turtle retinas using an array of 100, 1.5 mm long electrodes as reported previously.2,5 Furthermore we also used multi-electrode recordings that were obtained from mouse visual cortex (see Fig. 1). Depending on the experiment, the inter-electrode distance of the multi-electrode array used was either 200 or 400 µm. The signals from the electrodes were amplified, passed through a bandwidth filter (300–2000 Hz) and digitized with a resolution of 16 bits at a sampling rate of 32 or 40 kHz. Neural spike events were detected by comparing the instantaneous electrode signal to level thresholds set for each data channel using standard procedures described elsewhere.16,30,31 Thus, when a supra-threshold event occurs, the signal window surrounding the event is
1450009-2
March 24, 2014 15:4 1450009
Int. J. Neur. Syst. 2014.24. Downloaded from www.worldscientific.com by MONASH UNIVERSITY on 02/01/15. For personal use only.
NEV2lkit: A New Open Source Tool for Handling Neuronal Event Files
Fig. 1. Examples of artificial data and raw recordings. A trace was plotted each time the voltage crossed the lower and/or the upper thresholds. For the artificial data (upper panel), additional data sets with 40% (a) and 80% (b) added Gaussian noise are also showed (for better visual discrimination of the signals, just a single unit is plotted in the figure). Only for the data sets obtained from brain slice recording the spikes were peak aligned.
time-stamped and stored for later, offline analysis. All the selected channels of data as well as the state of the stimulus were digitized with a commercial multiplexed A/D board data acquisition system (Bionic Technologies, Inc) and stored digitally. In the case of the data files recorded from mouse visual cortex, spikes were aligned with respect to the peak above threshold before the unit sorting processing. Additionally we also used artificial data sets that were recorded from a commercially available 100 channel spike simulator (Bionics Technologies Inc) using the same acquisition hardware. These artificial data (see Fig. 1) consisted of sequentially generated spikes with five different amplitudes on each of the 100 recording channels, followed by a burst sequence on all channels. Furthermore, by adding 40% and 80% of Gaussian noise using the “addnoise” program from the TISEAN package,32 we generated two additional data sets that were used in the analysis (Fig. 1, artificial A, B).
All data sets were chosen because they contained a large number of spikes on many electrodes and represented typical experimental results.
2.1.1. Operational flow of NEV2lkit NEV2lkit is structured in two main modules, each designed to perform specific tasks of the neuronal event extraction and unit sorting process. The Input/Output module offers the opening, extraction and storage of spike event data. The spike sorting module incorporates two main steps. The first step is the conversion of the simultaneously recorded spikes into a clustering space constructed with the features of the signal shape extracted from applying the PCA. The second step attempts to separate this data into functionally meaningful “unit” clusters. After the program is launched, the graphical user interface of NEV2lkit provides easy access to all functions (see Fig. 2). The output routines allow
1450009-3
March 24, 2014 15:4 1450009
Int. J. Neur. Syst. 2014.24. Downloaded from www.worldscientific.com by MONASH UNIVERSITY on 02/01/15. For personal use only.
M. Bongard, D. Micol & E. Fern´ andez
(a)
(b)
(c)
Fig. 2. (Color online) Example of clustering results for a representative electrode in a multielectrode recording of rabbit ganglion cells to a moving bar. (a) Raw data before spike sorting. A total of 4622 events were recorded. The time scale is 1.6 ms. (b) Application of PCA to the data (spike vector) using the first three principal components for spike representation and a penalty function of 1, that penalizes overfitting by discouraging models with a large number of parameters. The sorted waveforms are shown in different colors. (c) The spike shapes corresponding to four putative neurons (same colors than b: red, green, blue and cyan). Numbers in parentheses show the total number of spikes. More distant neurons contribute to the background noise (gray).
to save data in NEV format or to export the time stamps of the spike event data as ASCII files with or without the voltage data of each individual spike. The sorting module operates on single or multiple files, allowing to automatically assign putative units over several experiments. Both modules use a predefined default or user defined parameters in the NEV2lkit GUI for their processing but run otherwise
in unsupervised mode. An exhaustive log file is saved to disk and provides the researcher all the information needed to reproduce the spike sorting process. A key feature of NEV2lkit is that it is able to process many files simultaneously. Therefore the assignments of units or clusters are done in many experiments at the same time. As a result the same
1450009-4
March 24, 2014 15:4 1450009
NEV2lkit: A New Open Source Tool for Handling Neuronal Event Files
unit on a given electrode has exactly the same unit label in all selected files. The number of experiments which can be processed simultaneously with NEV2lkit is only limited by the available RAM of the computer on which the program runs. In most cases, the software allows to process experiments with hundreds or even thousands of channels without any problem, therefore in order to speed up the sorting process there is no need to buffer to disk. Nevertheless, if needed, the source code can be easily modified to use the hard disk for temporary storage.
Int. J. Neur. Syst. 2014.24. Downloaded from www.worldscientific.com by MONASH UNIVERSITY on 02/01/15. For personal use only.
2.2. Spike sorting Spike sorting is an important step in the analysis of multi-electrode data. The main purpose is to distinguish between multiple cells recorded together.33–35 In this context, it should be taken into account that the extent to which the responses of a single neuron can be characterized depends on the ability to identify the action potentials originating from a single cell, and to discriminate these action potentials from other sources of electrical activity.36 This process, and therefore all further classification and analysis, depends on the noise present, the action potential amplitude, the conductive properties of the extracellular medium, the micro-electrode configuration, the localization with respect to the cell, and the type and quality of the micro-electrodes.37,38
2.2.1. Principal component analysis We apply a Principal Component Analysis (PCA) to reduce the dimensionality representation of the recorded waveforms in a p-dimensional space. Whereas detailed discussion of PCA is available elsewhere,29,39–41 here we briefly describe the basis of the method that we implemented for those unfamiliar with the technique. In the learning phase a set of spikes are accumulated, and from these the correlation matrix and its eigenvectors are computed. Then, a multidimensional matched filter using principal component vectors is defined and used in a real-time digital filter. The outputs of these filters can be used to classify the recorded spikes. Additionally, we obtain an ordered set of orthogonal basis vectors by calculating the a
eigenvectors of the variance-covariance matrix of the loaded data sets. Each principal component vector represents the most salient features of a given waveform and the eigenvectors with the largest eigenvalues correspond to the dimensions that have the strongest correlation in the dataset. Thus, the first eigenvector represents the largest correlation across the waveforms, the second eigenvector captures as much as possible of the waveforms differences not displayed by the first eigenvector, and so on. As a result, each spike is represented as a weighted sum of eigenvectors, corresponding to the largest eigenvalues of the data sets matrix. Since it is known that the signal-to-noise ratio can be optimized using the first M components of the eigenvectors in which the signal energy exceeds the noisy energy, we calculated the components that were significantly above background noise. We found that the first two components were often sufficient and that using additional components beyond the third one would not improve classification accuracy. These results agreed with previous studies reporting that one or two of the principal component vectors are sufficient to represent 80–90% of the spike waveform.42,43 Therefore, we only used the first three principal component vectors, which in our cases represent approximately 78% of the variation in the analyzed data, reducing the p-dimensional space to three dimensions. We performed the PCA on all multi-electrode channels and therefore the directions of maximum variation on all individual channels can be described simultaneously. Furthermore, we used different matrices for calculating the PCA vectors: the eigenvectors of the correlation one, and of the variance-covariance. 2.2.2. Clustering After the PCA has been calculated, each registered spike is represented as a point in a three-dimensional system spanned by the first three principal components. We used the Klustakwika implementation of the CEM algorithm of Celeux and Govaert to calculate the average separation of the spikes in the p-space. This algorithm computes all different spike groups in the hyperplane, and also takes into account
K.D. Harris, http://klusta-team.github.io/klustakwik/. 1450009-5
March 24, 2014 15:4 1450009
M. Bongard, D. Micol & E. Fern´ andez
the spreading of the spikes within the spike groups. As a result, the algorithm delivers the maximal separation between object groups (in our case different waveforms) combined with minimal spreading within the object groups (the same spike waveform). Klustakwik also allows to use Bayesian information in the clustering process.44 This option can be used to quantify the goodness of the resulting classification which is a useful aid to investigate the isolations of spikes in different clusters.45
Int. J. Neur. Syst. 2014.24. Downloaded from www.worldscientific.com by MONASH UNIVERSITY on 02/01/15. For personal use only.
2.3. Spike detection and extraction We performed spike detection in both supervised and unsupervised modes. In the supervised mode a lower and a higher detection threshold can be individually adjusted for each open data stream. For the extraction, a window time base for the detected events can be defined by the user (up to 4 ms). Extracted events are then displayed in this defined extraction window within the program. In automatic or unsupervised mode, NEV2lkit opens a suitable file, calculates the square root of the time-averaged (mean), squared value of the signal amplitude — the so-called root mean square or RMSE — of the individual signal streams, sets the upper and lower thresholds as RMSE×3.0 and extracts and displays the spikes over the defined time window. 3.
therefore be exactly evaluated. Each recording lasted 60 s and consisted of a sequence of single artificial spikes, each characterized by a different amplitude and time course of the signal, with an inter-spike interval of one second. The sequence of single spikes was followed by a burst of 10 identical spikes. For the analysis we selected randomly 10 electrodes from 20 different recordings. Thresholds on the recording channels were individually set allowing the extraction of different numbers of units on every single electrode. For further performance evaluation additional test sets were generated by adding 40% and 80% of Gaussian noise to this data (Fig. 1, upper panel). Results of the performance of the algorithms implemented in NEV2lkit for sorting the artificial data are summarized in Table 1, where the percentage of spikes used in the classification is expressed as the percentage of the total number of spikes. Erroneously classified spikes are given as a percentage of the total data set. The clustering algorithms were able to automatically remove outliers and artifacts from the data set and approximately 97% of all units were assigned to clusters. The best performance was achieved by choosing the maximization of the correlations of the basis vectors for the PCA. When the maximization of the variance-covariance was chosen to reduce the dimensionality representation of the recorded waveforms, more spikes were assigned to wrong units.
Results
The basic problems in spike sorting are illustrated in Fig. 2, which presents some waveforms recorded by a electrode in a typical multi-electrode experiment. A close inspection shows that the recorded spike waveforms usually contain a mixture of spikes from one or more cells and a significant amount of background activity from unknown sources. In a first step the quality of isolation can be tested by visually inspecting the superimposed spikes on each electrode, but even for relatively well-isolated spikes the recordings include some overlapping of spike events. Hence, the key problem is to correctly classify all the discriminable neurons.
Table 1. Summary of the clustering results for the artificial data. Method for PCA VarianceSum of Correlation covariance squares maximization maximization minimization
% of spikes used % error % of spikes used % error
96.7 0.28 96.7 0.28
3.1. Artificial data The advantage of simulated data is that all signal characteristics, such as the number of units within each recorded file, is known in advance and can
% of spikes used % error
1450009-6
96.7 0.56
0% noise 96.7 5.91 40% noise 96.7 5.91 80% noise 96.7 10.42
96.2 6.93 96.2 9.37 95.1 10.48
March 24, 2014 15:4 1450009
NEV2lkit: A New Open Source Tool for Handling Neuronal Event Files
These results were very reliable, even in the presence of 80% Gaussian noise. Consequently it seems that choosing the “right” calculation method improves significantly the accuracy of the clustering process.
Table 2. Clustering results for the experimental extracellular recordings. Method for PCA VarianceSum of Correlation covariance squares maximization maximization minimization
Int. J. Neur. Syst. 2014.24. Downloaded from www.worldscientific.com by MONASH UNIVERSITY on 02/01/15. For personal use only.
3.2. Extracellular recorded spikes While for experimental recordings the “real” number of unit clusters is not known, the quality of unit isolation can be judged, on a first approach, by looking at the superimposed spikes on every cluster. This delivers a first approximation of how many clusters exist within the data set recorded by an electrode. Figure 2(a) shows a typical data set from rabbit retina and the resulting classification into isolated units by using NEV2lkit (Fig. 2(b)). As we can see, even singular outliers are correctly sorted into an individual cluster (light-gray trace in Fig. 2(c)). The inter-spike interval histograms (ISI) and autocorrelation function histograms further document the high quality of the unit sorting carried out by NEV2lkit and show that there is no event assigned to clusters with an inter-event distance less than the refractory period. The overall qualification of the unit sorting performance for all multi-electrode recording data sets using different methods for the PCA is summarized in Table 2. In an attempt to quantify the goodness of the results, we calculated the ratio between the number of expected classes (from visual inspection of the recordings) and the classes isolated by NEV2lkit. Values higher than 1 mean that the unit sorting process is assigning less classes than the number suggested by visual inspection; 1 means a number of classes equal to the one assumed by visual inspection; and values lower than 1 represent a higher number of classes than expected. As it can be seen the results depend strongly on the PCA method. In order get more insight about the usefulness of every procedure, we also quantify the percentage of spikes classified incorrectly as noise. Due to the experimental setups used, noise from external sources were a rare event in our recordings and its percentage was less than 0.4% of the total signals. Therefore, the increase in the number of spikes assigned to noisy units is somehow related with wrongly assigned spikes. For each data set the noise units were determined and pooled, as well as the
Unit ratio % noisy units
Mouse brain 1 1 1.76 0.44
2 1.76
Unit ratio % noisy units
Rabbit retina 1.5 1 2.54 3.04
1 0.10
Unit ratio % noisy units
Turtle retina 1 0.75 3.61 0.45
1 0.02
Unit ratio % noisy units
Mouse retina 3 3 1.12 3.16
1 1.12
classification error expressed as the percentage of spikes incorrectly assigned to noise units. Taken together these results suggest that it is possible to find an optimal method for the classification of every data set. For example, the artificial data was best sorted into units using the maximization of the correlations as basis vectors for the PCA analysis, whereas for retinal ganglion cells recordings the best performance was achieved by choosing the maximization of the variance-covariances. Because all data were recorded with the same type of multielectrode array and electrophysiological equipment, one can assume that differences in the signal characteristics are indeed depending on the signal origin. 4.
Discussion
One of the most promising advances in neuroscience research is the development of multi-electrode arrays that allow multiple site recordings in various neural systems.7,8,46,47 A major challenge in this context is to acquire meaningful data from a large number of channels but a central problem is how many neurons are present in a given recording. Therefore, finding the spiking activity in electrophysiological recordings is seen as a first step in the decoding of neural activity. However there are usually a wide variability in the results of manual and
1450009-7
March 24, 2014 15:4 1450009
Int. J. Neur. Syst. 2014.24. Downloaded from www.worldscientific.com by MONASH UNIVERSITY on 02/01/15. For personal use only.
M. Bongard, D. Micol & E. Fern´ andez
even automatic spike sorting across a given data set among different subjects and laboratories which suggest the importance of good, widely available, automated spike sorting methods and tools. To facilitate this task we have developed an open source software which has proved to be useful for the classification and management of neural ensemble data. The main idea is to supply a graphical user interface that links the experimental data to a basic set of routines whose purpose is to distinguish between multiple cells recorded together and to assign detected spiking activity to specific neurons. An important issue in this field is to distinguish between single and multi-unit activity. Our results show that the procedures implemented in NEV2lkit using a combination of adjustable PCA- and CEMbased clustering can precisely detect single neurons even in the presence of large amounts of noise. Hence NEV2lkit offers a fast, reliable and accurate implementation for spike sorting. While it is known that template matching can deliver a more accurate classification of the signals,24,25 our results show that the PCA- and CEM-based approach is sufficiently accurate. Furthermore, one advantage of NEV2lkit lies in the fact that it operates fast and in an unsupervised mode over a large number of neuronal recording files so that the user only has to define the minimal and maximal class assignment limits. Thus, by using NEV2lkit it is possible, depending on the overall data size, to identify units across multiple experiments within seconds while the program runs with minimal required human intervention. Our results reveal that CEM-based clustering of the principal components calculated using the maximization of the correlations as basis vectors for the PCA works very well for synthetic data. Interestingly, it turns out that for the real extracellular recordings used in this analysis, the basis vector construction method was not delivering the best performance in the unit sorting process. As a result, the researcher has to define and choose the best method depending on the data and its particular features. NEV2lkit addresses this requirement by offering a selection of different basis vector methods to be used as base for the sorting process. This allows handling unit sorting problems which are presumably partly
depending on the signal origin. Since the physical multi-electrode array characteristics were similar in all the recordings used, our results suggest that there is an optimal PCA approach for spike sorting and unit assignment depending on every data set. Restrictions for the application of NEV2lkit are related to the assumption that extracellular spike shapes are always stationary. This expectation could be realized in our analyzed data sets, because our experimental setup delivers almost stationary recorded spike waveforms. However, there could be additional problems that limit the robustness of the unit sorting carried out by these procedures. For example, in the analyzed experiments the level of background noise remained relatively constant and electrode drift was not observed or was minimal. If these conditions are given, NEV2lkit suits apparently well for fast off-line unit sorting in multi-electrode recordings from different kinds of sources. We are aware that some functions currently in use, as well as advanced functions used in various labs working with multi-electrode data, are not included. We hope that the facilities to extend the code and add more specific algorithms and routines can compensate the former and encourages users to contribute providing ready-made functions useful to the scientific community. Thus, the modular program architecture allows an easy extension of its functionality, such as extending the number of supported data file formats or the integration of external libraries for spike detection. In this context there are already available versions of NEV2lkit that provide extended functionality to the feature extraction procedureb such as variable length data window, variability of the number of principal components and calculation and display of variance.48 Finally, NEV2lkit is free software distributed under the GNU General Public License.c This license grants it users legal permission to copy, distribute and/or modify the software. To this end, the source code is distributed alongside of executable binaries for Linux, OS X and Windows. The only restriction imposed on the user is that should they decide to redistribute the software in its original or modified form, they must in turn grant their users the same rights. We hope to attract an active community of
b http://neurobot.bio.auth.gr/nev2lkit. c
http://www.gnu.org/licenses/gpl.html. 1450009-8
March 24, 2014 15:4 1450009
NEV2lkit: A New Open Source Tool for Handling Neuronal Event Files
users who could benefit from and contribute to this project. Acknowledgments
10.
Int. J. Neur. Syst. 2014.24. Downloaded from www.worldscientific.com by MONASH UNIVERSITY on 02/01/15. For personal use only.
We would like to thank J. V. Segura for his valuable comments. In addition, we wish to thank M. Gerschner, L. van Ahrens and K. Gansel for providing some of the data sets used in this work, and V. Medran and S. Andreu for technical assistance. This work has been partially supported by Spanish Grants MAT2012-39290-C02-01 and IPT-20120574-300000, by ONCE (National Organization of the Spanish Blind) and by the Research Chair in Retinitis Pigmentosa Bidons Egara.
11.
12.
13.
References 1. M. Greschner, M. Bongard, P. Rujan and J. Ammermuller, Retinal ganglion cell synchronization by fixational eye movements improves feature estimation, Nat. Neurosci. 5 (2002) 341–347. 2. E. Fernandez, J. M. Ferrandez, J. Ammermller and R. A. Normann, Population coding in spike trains of simultaneously recorded retinal ganglion cellsx, Brain Res. 887 (2002) 222–229. 3. M. J. Schnitzer and M. Meister, Multineuronal firing patterns in the signal from eye to brain, Neuron 37 (2003) 499–511. 4. S. Takahashi, Y. Anzai and Y. Sakurai, Automatic sorting for multi-neuronal activity recorded with tetrodes in the presence of overlapping spikes, J. Neurophysiol. 89 (2003) 2245–2258. 5. G. J. Ortega, M. Bongard, E. Louis and E. Fernandez, Conditioned spikes: A simple and fast method to represent rates and temporal patterns in multielectrode recordings, J. Neurosci. Methods 133 (2004) 135–141. 6. N. V. Manyakov and M. M. Van Hulle, Decoding grating orientation from microelectrode array recordings in monkey cortical area V4, Int. J. Neural Syst. 20(2) (2010) 95108. 7. F. Franke, D. J¨ ackel, J. Dragas, J. M¨ uller, M. Radivojevic, D. Bakkum and A. Hierlemann, Highdensity microelectrode array recordings and realtime spike sorting for closed-loop experiments: An emerging technology to study neural plasticity, Front. Neural Circuits 6(105) (2012) doi: 10.3389/fncir.2012.00105. 8. G. Van Dijck, K. Seidl, O. Paul, P. Ruther, M. M. Van Hulle and R. Maex, Enhancing the yield of high-density electrode arrays through automated electrode selection, Int. J. Neural Syst. 22(1) (2012) 119. 9. M. A. Lebedev, J. M. Carmena, J. E. O’Doherty, M. Zacksenhouse, C. S. Henriquez, J. C. Principe
14. 15.
16.
17.
18.
19.
20. 21.
22.
1450009-9
and M. A. Nicolelis, Cortical ensemble adaptation to represent velocity of an artificial actuator controlled by a brain-machine interface, J. Neurosci. 25 (2005) 4681–4693. R. Gutierrez, J. M. Carmena, M. A. Nicolelis and S. A. Simon, Orbitofrontal ensemble activity monitors licking and distinguishes among natural rewards, J. Neurophysiol. 95 (2006) 119–133. L. R. Hochberg, M. D. Serruya, G. M. Friehs, J. A. Mukand, M. Saleh, A. H. Caplan, A. Branner, D. Chen, R. D. Penn and J. P. Donoghue, Neuronal ensemble control of prosthetic devices by a human with tetraplegia, Nature 442 (2006) 164–171. J. P. Donoghue, A. Nurmikko, M. Black and L. R. Hochberg, Assistive technology and robotic control using motor cortex ensemble-based neural interface systems in humans with tetraplegia, J. Physiol. 579 (2007) 603–611. N. G. Hatsopoulose and J. P. Donoghue, The science of neural interface systems, Annu. Rev. Neurosci. 32 (2007) 249–266. D. Gardner, Neurodatabase.org: Networking the microelectrode, Nat. Neurosci. 7 (2004) 486–487. U. Egert, T. Knott, C. Schwarz, M. Nawrot, A. Brandt, S. Rotter and M. Diesmann, MEA-Tools: An open source toolbox for the analysis of multielectrode data with MATLAB, J. Neurosci. Methods 117 (2002) 33–42. M. P. Bonomini, J. M. Ferrandez, J. A. Bolea and E. Fernandez, DATA-MEAns: An open source tool for the classification and management of neural ensemble recordings, J. Neurosci. Methods 148 (2005) 137–146. N. R. Luque, J. A. Garrido, J. Ralli, J. J. Laredo and E. Ros, From sensors to spikes: Evolving receptive fields to enhance sensori-motor information in a robot-arm, Int. J. Neural Syst. 22(4) (2012) 1250013. K. Ramanathan, N. Ning, D. Dhanasekar, L. Guoqi, S. Luping and P. Vadakkepat, Presynaptic learning and memory with a persistent firing neuron and a habituating synapse: A model of short term persistent habituation, Int. J. Neural Syst. 22(4) (2012) 1250015. W. K. Wong, Z. Wang and B. Zhen, Relationship between applicability of current-based synapses and uniformity of firing patterns, Int. J. Neural Syst. 22(4) (2012) 1250017. S. Ghosh-Dastidar and H. Adeli, Spiking neural networks, Int. J. Neural Syst. 19(4) (2009) 295–308. S. Ghosh-Dastidar and H. Adeli, A new supervised learning algorithm for multiple spiking neural networks with application in epilepsy and seizure detection, Neural Netw. 22(10) (2009) 1419–1431. S. Ghosh-Dastidar and H. Adeli, Improved spiking neural networks for EEG classification and epilepsy and seizure detection, Integr. Comput-Aided. Eng. 14(3) (2007) 187–212.
March 24, 2014 15:4 1450009
Int. J. Neur. Syst. 2014.24. Downloaded from www.worldscientific.com by MONASH UNIVERSITY on 02/01/15. For personal use only.
M. Bongard, D. Micol & E. Fern´ andez
23. K. Sameshima and M. Bach, Trends in multichannel neural ensemble recording instrumentation, in Methods for Neural Ensemble Recordings, ed. M. A. Nicolelis (CRC Press, Boca Raton, Florida, 1999), pp. 47–60. 24. R. A. Normann, D. J. Warren, J. Ammermuller, E. Fernandez and S. Guillory, High-resolution spatiotemporal mapping of visual pathways using multielectrode arrays, Vision Res. 41 (2001) 1261–1275. 25. S. Shoham, M. R. Fellows and R. A. Normann, Robust, automatic spike sorting using mixtures of multivariate t-distributions, J. Neurosci. Methods 127 (2003) 111–122. 26. U. Rutishauser, E. M. Schuman and A. N. Mamelak, Online detection and sorting of extracellularly recorded action potentials in human medial temporal lobe recordings, in vivo, J. Neurosci. Methods 154 (2006) 204–224. 27. S. Ghosh-Dastidar and H. Adeli, Wavelet-clusteringneural network model for freeway incident detection, Comput. Aided Civ. Inf. 18(5) (2003) 325–338. 28. M. Ahmadlou and H. Adeli, Enhanced probabilistic neural network with local decision circles: A robust classifier, Integr. Comput. Aid. E 17(3) (2010) 197–210. 29. E. Chan, V. Hok, A. Della-Chiesa, J. J. H. Miller, S. M. O’Mara and R. B. Reilly, Automated spike sorting algorithm based on Laplacian eigenmaps and k-means clustering, J. Neural Eng. 8(1) (2011) 016006. 30. Z. Tiganj and M. Mboup, A non-parametric method for automatic neural spikes clustering based on the non-uniform distribution of the data, J. Neural Eng. 8(6) (2011) 066014. 31. Y. Yuan, C. Yang and J. Si, The M-Sorter: An automatic and robust spike detection and classification system, J. Neurosci. Meth. 210 (2012) 281–290. 32. R. Hegger, H. Kantz and T. Schreiber, Practical implementation of nonlinear time series methods: The TISEAN package, Chaos 9 (1999) 413–435. 33. T. Takekawa, Y. Isoura and T. Fukai, Accurate spike sorting for multi-unit recordings, Eur. J. Neurosci. 31 (2010) 263272. 34. I. Cajigas, W. Q. Malik and E. N. Brown, nSTAT: Open-source neural spike train analysis toolbox for Matlab, J. Neurosci. 211 (2012) 45264. 35. C. Vargas-Irwin and J. P. Donoghue, Automated spike sorting using density grid contour clustering and subtractive waveform decomposition, J. Neurosci. Methods 164 (2007) 1–18. 36. N. Schmitzer-Torbert, J. Jackson, D. Henze, K. Harris and A. D. Redish, Quantitative measures of cluster quality for use in extracellular recordings, Neuroscience 131 (2005) 1–11.
37. H. Kaneko, H. Tamura and S. S. Suzuki, Tracking spike-amplitude changes to improve the quality of multineuronal data analysis, IEEE Trans. Biomed. Eng. 54 (2007) 262–272. 38. T. Takekawa, Y. Isoura and T. Fukai, Spike sorting of heterogeneous neuron types by multimodalityweighted PCA and explicit robust variational Bayes, Front. Neuroinform. 6(5) (2012) doi: 10.3389/fninf.2012.00005. 39. S. Ghosh-Dastidar, H. Adeli and N. Dadmehr, Principal component analysis-enhanced cosine radial basis function neural network for robust epilepsy and seizure detection, IEEE Trans. Biomed. Eng. 55(2) (2008) 512–518. 40. M. Al-Naser and U. Soderstrom, Reconstruction of occluded facial images using asymmetrical principal component analysis, Integr. Comput. Aided Eng. 19(3) (2012) 273–283. 41. A. Meraoumia, S. Chitroub and A. Bouridane, 2D and 3D palmprint information, PCA and HMM for an improved person recognition performance, Integr. Comput. Aided Eng. 20(3) (2013) 303–319. 42. R. Vollgraf, M. Munk and K. Obermayer, Optimal filtering for spike sorting of multi-site electrode recordings, Network 16 (2005) 85–113. 43. P. M. Horton, A. U. Nicol, K. M. Kendrick and J. F. Feng, Spike sorting based upon machine learning algorithms (SOMA), J. Neurosci. Methods 160 (2007) 52–68. 44. K. D. Harris, D. A. Henze, J. Csicsvari, H. Hirase and G. Buzsaki, Accuracy of tetrode spike separation as determined by simultaneous intracellular and extracellular measurements, J. Neurophysiol. 84 (2000) 401–414. 45. J. Wilda, Z. Prekopcsakc, T. Siegera, D. Novaka and R. Jech, Accuracy of tetrode spike separation as determined by simultaneous intracellular and extracellular measurements, J. Neurosci. Methods 203 (2012) 369–376. 46. M. R. Cohen and A. Kohn, Measuring and interpreting neuronal correlations, Nature 14(7) (2012) 811–819. 47. H. A. Wark, R. Sharma, K. S. Mathews, E. Fernandez, J. Yoo, B. Christensen, P. Tresco, L. Rieth, F. Solzbacher, R. A. Normann and P. Tathireddy, A new high-density (25 electodres/mm3) penetrating microelectrode array for recording and stimulating sub-millimeter neuroanatomical structures, J. Neural Eng. 10 (2013) 045003. 48. D. A. Adamos, E. K. Kosmidis and G. Theophildis, Performance evaluation of PCA-based spike sorting algorithms, Comput. Methods Programs Biomed. 91 (2008) 232–244.
1450009-10
