Robust detection of neural spikes using sparse coding based features

1.   Introduction

Neurons are the elementary structure and functional units of nervous system. Monitoring and studying the neuronal activity can provide a useful mean for understanding the interactions between neurons and higher brain functions ^[1,2,3]. Early researches on neuroelectrophysiology have discovered that neurons represent and transmit information by firing electric impulses, called action potentials or, briefly, spikes ^[4,5,6]. According to this significant discovery, electrode intracranial recording technique has been widely used to collect the extracellular activity of neurons for decades ^[7,8,9].

The analysis and processing of extracellular recording signals are generally summarized as three main steps (see Figure 1). The first and crucial step is the detection of spikes from the recorded signals, which promises to be able to extract necessary spike data for all subsequent analyses ^[10]. Simply speaking, the aim of spike detection is to accurately find and localize the occurrence of individual spikes. To this end, in the past two decades, various algorithms have been developed with their own advantages and drawbacks.

The early works about spike detection mainly focus on the traditional simple characteristics of spikes. An example is the most widely used amplitude thresholding (AT). This technique only takes signal amplitude as the feature for the identification of spikes ^[11,12]. Nonlinear energy operator (NEO) is another popular detection method, which takes into account the instantaneous energy of the signal ^[13,14]. The major attraction of these methods is their simplicity and low computational complexity, and hence, they are easy to implement in hardware. However, it is proven that these traditional simple characteristics based methods are sensitive to background noise.

To address this problem, researchers pay more attention on the complex structural characteristics of spikes. Template matching (TM) ^[15,16] and wavelet transform (WT) ^[17,18] are the most representative techniques. TM searches for spike events in the signal that "closely resemble" the template. In other words, the similarity between template and the signal determines whether a spike occurs. It has been demonstrated that TM performs much better than the simple threshold algorithms and can identify spikes with different shapes ^[19]. However, it does not work well with the presence of high-level background noise. The core idea of wavelet-based technique is to separate spike signals from background noise based on the wavelet coefficients. Although this technique performs well even in the high noise level, its performance is sensitive to the selected wavelet basis that is not optimal for the spikes from different neurons.

As mentioned above, while the existing methods for spike detection have made great progress, it is desirable to develop a more reliable spike detection method that can satisfy the following demanding requirements at the same time. Firstly, it should be robust to the background noise from various sources, such as ambient environment, recording hardware and the activities of distant neurons. Secondly, as the spikes from different neurons have different shapes, the detector should have strong flexibility in the spike shape. Considering the above requirements, in this work, we propose a novel spike detection algorithm based on the promising sparse representation theory. Combining the internal generation mechanism of neural spikes and the characteristics of extracellular neural recordings, we first propose a target-driven sparse representation framework for extracellular recording signals. Under this framework, the neural spike signals can be effectively separated from background noise. In addition, to improve the flexibility of the proposed method, a universal dictionary is learned to give a sparse representation of various spike signals. Finally, the information (location and number) of spikes in the recorded signal are achieved by comprehensively analyzing the sparse features.

The rest of this paper is organized as follows. Section 2 gives a brief review of relevant background knowledge. Then, the proposed method is detailedly analyzed in Section 3. In Section 4, the proposed method is evaluated on both synthesized and real data. Finally, the conclusion is drawn to summarize the work of the paper in Section 5.

2.   Background

In this section, we briefly introduce the theory of sparse representation to provide the relevant background for later section.

Sparse representation has received increasing interests in recent years. As a powerful signal processing technique, it has been widely studied and successfully used in various applications, such as recognition ^[20], denoising ^[21], reconstruction^[22], segmentation ^[23] and so on. The major attraction of sparse representation is its ability to approximately represent the target signal using a few basic elements^[24]. Mathematically, suppose that $\mathrm{y}\in \mathrm{R}^L$ is the signal to be processed, and $\mathrm{D} = [\mathrm{d}_1, \mathrm{d}_2, \ldots, \mathrm{d}_N]\in \mathrm{R}^{L\times N}$ is a given dictionary with $N$ atoms. The sparse representation of $\mathrm{y}$ with respect to $\mathrm{D}$ is to search for a sparse vector $\mathrm{\alpha} = [\alpha_1, \alpha_2, \ldots, \alpha_N]^T\in \mathrm{R}^N$ such that

Technically, the sparsity of $\mathrm{\alpha}$ can be measured by the $l_0$-norm of it, which is equivalent to the number of nonzero components in $\mathrm{\alpha}$. In view of this fact, the problem of sparse representation can be cast as the following $l_0$-minimization problem,

Considering the fact that the received signal is inevitably contaminated with noise, the linear constraint $\mathrm{y} = \mathrm{D\alpha}$ is often relaxed to the term $\|\mathrm{y}-\mathrm{D\alpha}\|^2_2\leq \epsilon$, where $\epsilon$ is the error tolerance. Then, the problem (2.2) can be rewritten as,

Alternatively, $\mathrm{\alpha}$ can also be obtained by

where $\lambda$ is a nonzero constant that balances the dual objectives of minimizing both accuracy and sparsity.

To address the above minimization problem, researchers have made great efforts and various algorithms are proposed. According to the working principle of the algorithm, the existing algorithms can be roughly classified into two categories. The first category includes algorithms that are developed based on greedy principle and can be directly used to solve problem (2.4), such as matching pursuit (MP), orthogonal matching pursuit (OMP) etc ^[25,26,27]. The algorithms fall into the second category mainly take advantage of convex relaxation to replace problem (2.4) with a convex optimization problem, including basis pursuit (BP), interior-point methods, gradient methods and so on ^[28,29,30].

3.   Proposed spike detection method

In this section, we present a novel method based on the theory of sparse representation for the task of spike detection. Firstly, we build a target-driven sparse representation framework to separate the target spike signal from background noise. Moreover, to give an effective representation of various spikes, we further construct a universal dictionary. Finally, the information of spikes in the recorded signal can be determined by analyzing the sparse coefficients. Figure 2 gives an intuitive illustration of our method. Details of the proposed method are presented as follows.

3.1.   Sparse model of extracellular recording signals

In general, the characteristics of extracellular recording signals can be summarized in two aspects. Firstly, the single extracellular recording signal usually contains the spikes with different shapes and sizes from a number of neurons near electrode tip. Secondly, the recorded signal is inevitably contaminated by the interference noise. In view of these facts, the recorded signal $\mathrm{y}\in \mathrm{R}^L$ can be represented as a linear combination of target spike signal and background noise,

where $\mathrm{y}^s$ is the clean spike signal, and $\mathrm{y}^n$ represents the background noise.

It should be noted that the recorded signal contains a large number of background noise from various sources, such as ambient environment, recording hardware and the activities of distant neurons, which make the problem of spike detection difficult. Inspired by the theory of sparse representation, we then take advantage of this famous technique to reduce the influence of background noise. That is

Under this framework, we intend to extract the sparse coefficients as features for the identification of spikes. Note that its inherent assumption is that the spike signal is sparse over the dictionary, and the noise is not sparse over the dictionary. However, we find that this assumption is difficult to be full satisfied because of the diversity of background noise. In particular, part of noise (such as neural noise) is potentially spike-like and may be generate a sparse projection with the dictionary. As a result, it will result in performance degradation to a certain degree. This fact will be given a detailed explanation in the Section 3.3.

To address this problem, we here construct two different dictionaries to express $\mathrm{y}^s$ and $\mathrm{y}^n$, respectively, that is

where $\mathrm{\alpha}^{s}$ is a sparse vector specific to spike signal $\mathrm{y}^{s}$ and $\mathrm{\alpha}^{n}$ is a sparse vector specific to noise signal $\mathrm{y}^{n}$. Moreover, $\mathrm{D}^s$ and $\mathrm{D}^n$ are dictionaries associated with sparse vectors $\mathrm{\alpha}^{s}$ and $\mathrm{\alpha}^{n}$, respectively. Therefore, (3.1) can be further written as

Note that with the help of coupled dictionaries, the noise is prone to be absorbed by the noise dictionary. Therefore, the sparse vectors $\mathrm{\alpha}^{s}$ and $\mathrm{\alpha}^{n}$ can capture the difference between the spike signal and background noise. When a spike occurs in the recorded signal $\mathrm{y}$, it will lead to a large coefficient in $\mathrm{\alpha}^{s}$. On the contrary, if there is no spike in $\mathrm{y}$, the nonzero coefficients are mainly scattered in $\mathrm{\alpha}^{n}$. In this way, we can determine whether there is a spike in the recorded signal.

3.2.   Dictionary construction

As stated in ^[31,32], the dictionary learned from data always results in better performance than the off-the-shelf dictionaries. In view of this fact, in this work, we use the real data as template to construct the necessary dictionaries $\mathrm{D}^s$ and $\mathrm{D}^n$.

Since the subject of this work is to detect spikes in the recorded signal, we try to find a signal dependent dictionary to represent and discriminate spike signals. With this consideration, we make use of the similarity of the spike signals, and propose to directly use spike waveforms to form the dictionary $\mathrm{D}^s$. In this work, we utilize the classical $K$-means algorithm ^[33] to learn the dictionary $\mathrm{D}^s$ from the precollected spike data, that is $\mathrm{D}^s = [\mathrm{c}_1, \mathrm{c}_2, \ldots, \mathrm{c}_K]$, where $\{ {\mathrm{c}_i}\}_{i = 1}^{K}$ denotes the $K$ clustering centers. In summary, there are two mainly advantages with the application of $K$-means algorithm. Firstly, it will produce a compact dictionary only with a few atoms, which results in the decrease of computational burden in the subsequent sparse optimization. Secondly, it maintains the diversity of atoms, and ensures that the suitable atom can be properly found and assigned to each spike to be distinguished.

Different from the construction of $\mathrm{D}^s$, we use the structural features of the signal as atoms to construct the dictionary $\mathrm{D}^n$. Similarly, we first extract large numbers of noise signal segments from the recorded signals to construct the necessary noise data. Note that the length of noise signal segments is set to the size of neural spikes. Then, the $K$-SVD algorithm ^[34] is introduced to extract the structural features of the collected noise samples, and denoted by $\mathrm{u}_1, \mathrm{u}_2, \dots, \mathrm{u}_{N_n}$. With these structural features, dictionary $\mathrm{D}^n$ can be obtained and formulated as $\mathrm{D}^n = [\mathrm{u}_1, \mathrm{u}_2, \dots, \mathrm{u}_{N_n}]$.

3.3.   Sparse feature extraction

In this subsection, we will show how to extract sparse coefficients of the recorded signal with respect to the learned dictionaries. Due to the limitation of atom size, the recorded signal is piecewise processed with a sliding window. In addition, the length of the sliding window should correspond to the size of atoms.

Let $\mathrm{y}_k$ be the $k$th segment of the recorded signal $\mathrm{y}$. Based on the equation (3.4), $\mathrm{y}_k$ can be represented as

where $\mathrm{D}$ is the combination of dictionaries $\mathrm{D}^s$ and $\mathrm{D}^n$, and $\mathrm{\alpha}_k$ is the column vector consisted of $\mathrm{\alpha}^{s}_{k}$ and $\mathrm{\alpha}^{n}_{k}$. Meanwhile, $\mathrm{\alpha}_k$ can be can be estimated by

To address above optimization problem, we use MP algorithm as the solver for its simplicity and easy implementation. Figure 3 gives an intuitive illustration of sparse coefficient distribution. From Figure 3(c) and (e), it can be seen that the spike in the signal will induce a large coefficient in $\mathrm{\alpha}_{k}^{s}$. By contrast, if no spike occurs in the signal, the numerical values of coefficients will be small and mainly locate in $\mathrm{\alpha}_{k}^{n}$. Therefore, we take the maximum coefficient of $\mathrm{\alpha}_{k}^{s}$ as the feature to distinguish the spike from background noise. In order to facilitate subsequent analysis, we save the maximum coefficient $\alpha_{kMax}$ into a new vector $\mathrm{Z}$, i.e., $\mathrm{Z}(k) = \alpha_{kMax}$. After applying the same operation to all segments, the coefficient signal $\mathrm{Z}$ is finally obtained.

Figure 3 also gives an intuitive explanation about the effectiveness of noise dictionary $\mathrm{D}^n$. From Figure 3(b), it can be seen that if we only use spike dictionary $\mathrm{D}^s$ to build the sparse model, the background noise may still yields a sparse projection with respect to $\mathrm{D}^s$. However, the results in Figure 3(c) demonstrate that by introducing the noise dictionary $\mathrm{D}^n$, the coefficients of noise mainly locate in the projection of $\mathrm{D}^n$. Moreover, from Figure 3(d)-(e), it can also be seen that the noise is prone to be absorbed by the noise dictionary. Therefore, we can conclude that the combination of spike dictionary $\mathrm{D}^s$ and noise dictionary $\mathrm{D}^n$ contributes to the improvement of detection performance.

3.4.   Spike discrimination

In order to detect the occurrence of spikes, an amplitude threshold is necessary to be set. If the coefficient of the resulting signal $\mathrm{Z}$ exceeds the threshold, a neural spike is identified. In this work, the threshold is determined based on the statistical characteristics of the resulting signal $\mathrm{Z}$, and is set as follows:

where $\bar{\mathrm{Z}}$ is the mean of output signal $\mathrm{Z}$, and $b$ is a constant ($3 \leq b \leq 5$) that scale the threshold.

4.   Experimental results

In this section, we report our experimental results of applying the proposed method into spike detection. The proposed method is firstly tested on synthesized data to verify the flexibility in the spike shape and the superior detection ability at different noise levels. And then we apply the propose method to real extracellular recording signals.

Experimental data: Both synthesized and real data are used to verify the performance of the proposed spike detection method.

(1) Synthesized extracellular recording signals: The synthesized signal consists of two types of signals: target spike trains and background noise. In each spike train, spike events are randomly placed based on Poisson distribution with a certain refractory period. Background noise is extracted from real recorded signals. In order to make a more intuitive explanation, we show an example of synthesized extracellular signal generation in Figure 4.

(2) Real data: It is collected from the prefrontal cortex of mice when performing a working memory task ^[35]. All acquisition and processing procedures were accomplished by the institute of Neuroscience, Chinese Academy of Sciences. The recorded signals are high-pass filtered at $250$ Hz by a digital Butterworth filter and sampled at $40$ kHz.

Evaluation criteria: Hit rate (HR) and false rate (FR) are used as the criteria to measure the performance of the algorithms. HR measures the detection accuracy and FR represents the error between ground truth and detection result. The higher HR and lower FR indicates better performance of the spike detection method. These two criteria are defined as follows ^[36]:

where $N_{c}$ is the number of correctly detected spike events (true positive), $N_{t}$ is the total number of true spike events recorded in the signal, $N_{f}$ is the number of falsely detected spike events (false positive) and $N_{d}$ is the total number of spike events detected by the method.

4.1.   Experiments on synthesized data

4.1.1.   Flexibility Analysis

In this subsection, an experiment is first carried out to verify the flexibility of our method. To this end, we use about $100$ neural spikes (not included in spike data for dictionary construction) to synthesize a set of test signals, in which each test signal only contains one neuron's spike train. In addition, we add different levels of background noise to the synthesized signals, and the noise level ranges from $0.10$ to $0.45$.

Figure 5 shows the average performance of the proposed method, including average HR, average FR and the corresponding deviations (Dev). As observed from Figure 5, the proposed method performs with higher HR and lower FR under all different scenarios. Furthermore, it can be seen that the standard deviation values obtained by the proposed method for all the cases are always small, which means that the performance of the proposed method is steady and believable. These facts indicate that the proposed method has good generalization performance and is able to detect distinct spike trains effectively.

4.1.2.   Performance Comparison

To verify the superior detection ability of the proposed method, we compare it with three classical algorithms (AT, NEO and TM mentioned in Section 1) and two latest ones (SWT ^[37] and NLWD ^[18]). In addition, we also make a comparison with a widely used sparse denoising method (SOLD) ^[38]. In order to make a fair comparison, the parameters of each method are properly tuned so as to achieve the best possible result. Considering the fact that the accuracy of spike detection will directly impact all subsequent analyses, we here define the "best possible result" as high HR and least possible FR.

Since the microelectrode generally picks up the electrical activities from more than one neuron, we analysis the performance of the algorithm with more realistic synthesized signals. In this experiment, we use five most representative spikes (see Figure 6) to generate the target spike trains in the test signals. For each noise level, the average results over $50$ trials are obtained for all methods.

The performance comparisons of seven methods are shown in Figure 7. From Figure 7(b), it can be seen that the FR of the proposed method is less than those of other methods, which implies that the proposed method is more accurate. On the other hand, the results in Figure 7(a) demonstrate that the proposed method achieves higher HR compared with the other six methods, especially with the presence of high-level noise. It should be mentioned that our method still maintains over $63 \%$ HR even if the noise level increases to $0.45$, whereas the HR of other six methods is less than $55 \%$. These findings indicate that the proposed method is able to detect more spikes from the recorded signal, which contributes to extracting more neurons' information for future researches.

4.2.   Experiments on real data

As we know, real neural signals always differ from synthesized signals due to many uncontrollable internal or external factors. Therefore, in this subsection, we use real data to test the performance of the proposed method. The detailed description of the real signals can be seen Experimental data.

In fact, it is hard to evaluate the performance of a spike detection method with real recorded signals because the ground truth information (such as the number and locations of spike events) is unknown to us. In order to test the performance of the proposed method on the real signals, the first task is to obtain the ground truth information for error measurement. To this end, two experienced EEGers are invited to locate the spike events in real signals. It should be noted that only the signal waveforms which are labeled by both EEGers are recognized as spike events.

In order to make an intuitive analysis, all spike detection methods are first applied to a slice of real recorded signal and the outputs are shown in Figure 8. For better observation, the true spikes in the original signal are labeled with blue asterisks, and the spikes detected by the method are labeled with red circles. From Figure 8, we can see that the proposed method can effectively identify spike events from background noise. In addition, we further compare the proposed method with other established ones on several real recorded signals and present the results in Table 1. It can be seen that the proposed method performs much better than other six ones. Therefore, we can conclude that the performance of the proposed method is superior to other established ones.

5.   Conclusions

In this work, a novel method based on sparse representation was introduced for high accuracy spike detection. The proposed method mainly focused on how to improve the robustness to noise and the flexibility in the spike shape. To this end, we learned the dictionaries from the precollected datasets to give a sparse and comprehensive representation of the extracellular recording signals. Under the constructed sparse representation framework, sparse coefficients can be used as the features for the discrimination of spikes. Experimental results on both synthesized extracellular neural recordings and real data demonstrated that the proposed method outperforms the existing methods in terms of both robustness and flexibility.

Acknowledgments

This work was supported by the Fund of Science and Technology Department of Guizhou Province (No. [2020]1Y253) and the Fund of Guizhou University of Finance and Economics (No. 2019XYB02).

Conflict of interest

The authors declared that they have no conflicts of interest to this work.