Data generation for simulation

The synthetic spike generation procedure can serve two purposes: firstly, to provide simulation data for examining the effects of the approach proposed in this study, and secondly, to provide training data for learning a classifier to identify overlapping spikes.

To generate the simulation data, we first assumed that raw extracellular signals would be recorded through a single electrode at a sampling rate of 30 kHz, with a measurement duration limited to 60 s. We initially set the number of neurons to three, assuming that when overlapping spikes occur near this electrode, at least two single units would fire spikes simultaneously or within a short delay. Additionally, it was assumed that the amplitude of each neuron’s spike waveform would decrease according to an inverse-square law with increasing distance from the electrode [26, 27]. In other words, we determined the magnitude based on the inverse-square law, which is expressed as:

$$_=\frac_}_\right)}^}$$

(1)

where Ij denotes the magnitude of the j-th neuron’s spike waveform, dj is the distance from the electrode tip to neuron j, and I0,j is the magnitude when d = 0. The magnitude was set to I0 = 120 μV when fully contacting the electrode tip, based on the actual action potential morphology of the extracellular recordings. Initially, the distances between neurons and electrodes were randomly sampled between 20 and 60 µm to align with the number of spike templates. These distances are then randomly shuffled in each simulation and applied element-wise as weights to the spike templates, which have been rescaled between 0 and 1. Here, we ensured that at least two neurons were placed a short distance apart to reproduce overlapping spikes. Furthermore, we assumed that the ISI of the spike train firing from each neuron strictly followed a gamma distribution [28, 29].

Spike templates, representing the action potential of an isolated single unit, are necessary to generate plausible synthetic spike trains (as detailed in “Spike template estimation” section). To support this, we postulated that the spike templates would be the action potentials of single units, detectable around the electrode tip. In this study, initial spike templates were obtained from a real dataset (see “Real data acquisition” section). These spike templates were convolved with the spike timings to generate simulated spike trains, assuming the firing activities of an ideally isolated single unit. Since the ISIs of real neuronal spike timing are known to follow a gamma distribution, we generate the gamma random numbers to represent the ISIs. The shape and scale parameters of the gamma distribution define the spike firing rate, fj, using the following equation:

where αj represents the shape parameter and θj represents the scale parameter of the single unit j. Considering the refractory period for each single unit, we initialized α with a randomly sampled value from a uniform distribution between 1.01 and 2 [1]. The firing rates were set to 60 Hz to ensure that the probability of overlapping spikes occurring exceeded 20%, resulting in each θ being approximately 0.05 (see Figure S1 in Supplementary Material). Note that the firing rate was set uniformly to ensure a consistent proportion of overlapping spikes among all detected spikes (see Figure S2 in Supplementary Material).

We generate random numbers representing ISIs following the gamma distribution based on these parameters. In this process, we use the “gamrnd” function from the Matlab toolbox. Following the ISI generation, samples within the 48-point refractory period are excluded, reflecting the physiological phenomenon where subsequent spikes do not occur within 2.5 ms in the same neuron [1]. The spike trains are then created by taking the cumulative sum of the ISI samples. Spike templates are convolved with the spike trains randomly generated as the number of spike templates individually.

These ideally isolated single unit activities are then summed across multiple spike templates. This allows us to access information regarding the occurrence of overlapping spikes and the unit labels involved. Finally, Gaussian white noise with a signal-to-noise ratio (SNR) of a specific range (from 1 to 3 detailed in “Assessment” section) is applied. The SNR is defined as the minimum peak-to-peak spike waveform scale relative to the root mean square of the spike-free noise segment, as detailed in [30]. Note that in this study, we did not consider background signal drift.

Real data acquisition

A real dataset was obtained by chronically implanting a 96-channel microelectrode array into the primary motor cortex (M1) of one Rhesus macaque. Spontaneous neural activities were recorded while the monkey freely moved its arm without performing any task instruction. In this study, we evaluated the data obtained from only channels in which distinct neural spikes were observed. Neural signals were acquired using a Cerebus system (Blackrock Neurotech, Salt Lake City, UT) for a duration of 180 s at a sampling rate of 30 kHz.

Band-pass filtering and spike detection

The preprocessing procedures for extracellular activity in our proposed approach follow the traditional spike sorting pipeline. We constructed a 6th-order Butterworth filter for band-pass filtering. Specifically, we filtered the raw extracellular signals through a 300–3000 Hz band-pass filter. Subsequently, for spike detection, we employed the method proposed by Quiroga et al. [31]. According to their method, the threshold, thr, is determined based on the standard deviation of the background noise of the filtered signal, x, as defined in Eq. 3.

$$thr=c\cdot \text\left(\frac\right|}\right)$$

(3)

where the denominator, 0.6745, is derived from the inverse of the cumulative distribution function for the Gaussian distribution, and c is the constant reflecting spike detection sensitivity, which can be calculated as std(x)/median(|x|). After the threshold was crossed, we stored the putative spike waveforms, U =  ∈ ℝ44×M, with 44 samples (equal to 1.5 ms) for each detected spike, where M denotes the number of detected spikes. We aligned spike waveforms based on their minimum peaks, with 10 samples preceding and 34 samples following its peak latency.

Spike template estimation

To estimate the spike template, we implemented a two-step procedure involving LDA–GMM and iForest. LDA–GMM is an iterative subspace learning method that combines LDA and GMM clustering [22]. In simple terms, it iteratively updates the discriminant prediction matrix using GMM until achieving maximum cluster separability. The objective function to be maximized can be expressed as follows:

$$\underset}^}^}}Trace\left(\frac}^}_\text}}^}_\text}\right)$$

(4)

where L represents the labels clustered by GMM, and Sb and SW are the covariance matrices representing the between-class and within-class variances, respectively. The superscript T denotes the transpose of the matrix. Thus, this algorithm can identify single unit spike candidates by detecting outliers. However, the efficacy of outlier detection may be affected by the initialization of GMM, owing to the intricate multimodal distribution of data within the subspace. Therefore, we additionally applied the iForest, which is one of the unsupervised anomaly detection methods, to potentially find spike templates close to single unit spikes, further enhancing the quality of the clustered spike waveforms. The iForest constructs an ensemble of isolation trees (iTrees) from the datasets, identifying anomalies as instances with shorter average path lengths in the iTrees [24]. With path lengths, h(u), anomaly score, sa, reflecting the degree of anomaly is given as follows:

$$_\left(\mathbf,n\right)=^\left(\mathbf\right)\right)}}$$

(5)

where E(h(u)) is the average of h(u) for collected isolation trees, and c(n) is the normalizing factor, which can be described as the average of h(u) given node n. The outlier score ranges between 0 and 1, with values closer to 1 indicating a higher likelihood of being an outlier.

Initially, we differentiated each detected spike waveform with respect to time (du/dt), followed by the application of LDA–GMM [22]. This approach was chosen because it is well-known that differentiated waveforms tend to outperform undifferentiated ones [22]. Each clustered waveform group was refined, by detecting and eliminating anormalies through the iForest. At this point, we optimized the contamination fraction of the iForest by calculating the unimodality of the feature distribution through four-fold cross-validation for each waveform group. The spike templates were then estimated by taking the median of the waveforms for each waveform group (see Fig. 2E).

The number of spike templates confirmed from synthetic spikes can be equal to or fewer than the predefined number (three units in this simulation) during the generation of synthetic spike data, depending on the specified changes in SNR. However, to estimate spike templates from real neural spikes, we needed to determine the number of potentially valid clusters (to be used as spike templates). We employed the LDA–GMM method, gradually increasing the number of clusters (from 2 to 5) according to the method proposed by Keshtkaran and Yang [22]. Subsequently, we examined whether the data distribution in the subspace was over-clustered using the Anderson–Darling test. This method aligns with optimizing the contamination fraction of the iForest. All these procedures were implemented based on the “iforest” function from the Matlab toolbox and the LDA–GMM toolbox available in [22].

Synthetic spike generation for classifier

When building a classifier, it is crucial to reconstruct synthetic spike data (or simulation data) with characteristics as similar as possible to the observed data. Initially, we estimated the parameters of the gamma distribution representing the ISI distribution of the observed data using maximum likelihood estimation, performed using the “gamfit” function in the Matlab toolbox. However, if the SNR is low, fewer spikes may be detected, potentially leading to fitting failure. So, in this scenario, we randomly sampled alpha between 1.01 and 2, the same as the simulation data generation process (see “Data generation for simulation” section). By utilizing the estimated parameters of the gamma distribution, we estimated the firing rate according to Eq. 2. Using these parameters of the gamma distribution and the spike templates estimated from the observed data, we generated synthetic spike data following the method mentioned in “Data generation for simulation” section. To minimize heterogeneity with the observed data, the proportion of overlapping spikes in the synthetic spike data was maintained to be similar to those of the observed data, and the SNR was set to the SNR estimated from the observed data.

Classification of overlapping spikes

To classify overlapping spikes from detected spikes, ground truth is necessary to train the classifier. In this section, we detail constructing the classifier using synthetic spikes to identify overlapping spikes and its testing with real detected neural spikes.

Feature extraction

Projecting detected spikes into a low-dimensional subspace is beneficial as it efficiently summarizes essential information. For classifier training, we first projected the detected spikes into a low-dimensional space using principal component analysis, which can be expressed as:

$$\mathbf=}^\overline }$$

(6)

where z is the score of the principal components, CT ∈ ℝD×M denotes the transposed loading matrix, and \(\overline }\) is the centralized waveforms. We determined dimensionality, D, by identifying the minimum number of components needed to account for at least 90% of the data variance. This determination was based on the eigenvalues’ contribution the total variance of the training data, which on average results in D being 15. To ensure dimensional consistency across both datasets, it is important to perform subspace learning when the SNR and single unit waveforms of synthetic and real neural spikes are similar.

Classification

Since the principal component scores for synthetic spikes with the label information of overlapping spikes exhibit inherently non-linear class distributions, we opted to construct a support vector machine (SVM) with a radial basis function (RBF) kernel [32]. We set the box constraint to 1, a value optimized to achieve the highest possible classification accuracy.

Particle swarm optimization algorithm

Based on the model illustrated in Fig. 1, the process of decomposing overlapping spikes requires the optimization of delay time, τ, and coefficients (also known as contribution index), w, for each single unit j. The objective function, aimed at minimizing the mean absolute error of this model, is formulated as:

$$\underset^^}}\frac\left|\sum_^\sum_^__\left(t-_\right)-_\left(t\right)\right|$$

(7)

where J denotes the total number of the spike templates, uo is the observed overlapping spike. If the number of single units is two, overlapping spikes could be modelled as simple combinations based on time delay changes and partitioned using cross-correlation, etc. However, as the number of combinations increases exponentially with the number of single units, we estimated the variables of the overlapping spike model using Particle Swarm Optimization (PSO), a metaheuristic search algorithm that iteratively explores the solution space [25]. PSO operates based on a set of candidate solutions, called particles, and optimizes the solution by adjusting the particles’ positions and velocities within the search space. In each iteration, PSO evaluates the movement of each particle to find a better position, which is considered the optimal solution. The estimated variables include τ and w, each of which is constrained within a certain range. Specifically, we limited the τ to be between − 44 and 44, and w between 0 and 1, depending on the number of estimated spike templates. To prevent excessive exploration of the search space, the inertia weight damping ratio was set as 1. Additionally, to limit the movement of particles in each iteration, the inertia weight, wχ, was set as 0.55. Additionally, the number of PSO particles was fixed at 100, and the algorithm was repeated 1000 times for each spike. This represents the experimental number of iterations in which the objective function can be minimized. If the positions (or variables) are no longer updated, the iteration stops.

Assessment

We evaluated the proposed method’ sensitivity to noise by changing SNR from 1 to 3 with 0.1 interval. To assess the impact of noise, we iteratively regenerated simulation data up to 1000 times. We then compared the proposed method with the following approaches: “outlier detection using LDA–GMM (LG)”, “outlier detection using iForest + LDA–GMM (IF+)”, and “testing the subspace scores of observed spikes using a synthetic spike subspace score-based classifier with PCA ”.

To determine iForest’s effectiveness on template spikes, we calculated the root-mean square error (RMSE) between the refined and the ground-truth spike waveforms and compared it to the distribution of all spike errors without iForest. This can be expressed in the following formula:

$$RMS_=\sqrt\sum_^}_-}}_\right)}^}$$

(8)

where \(}}_\) denotes the refined spike waveform of the m–th spike within the single unit j. Based on this, we quantified the relative iForest effect by defining it as the difference in RMSE between IF+ and LDA–GMM. The iForest effect is expressed as follows:

$$\text=\frac}_}^-}_}}}_}^+}_}}$$

(9)

where \(}_}^\) and \(}_}\) denote the RMSE of LG and IF+, respectively. This effect was represented as a function of SNR.

In addition, since overlapping spikes occur in about 20% of total spikes (see Figure S1 in Supplementary Material), we quantified precision, which represents the ratio of predicted overlapping spikes among actual overlapping spikes. Additionally, recall, which indicates the probability that the predicted overlapping spike contains an actual overlapping spike, and the F1-score, which quantifies performance by considering the trade-off between precision and recall, were also calculated. Statistical tests for classification performances were conducted using the Wilcoxon rank sum test and Kruskal–Wallis test, with a Tukey–Kramer correction applied for the post-hoc multiple comparison testing. In addition, we compared the effects of the proposed approach, and LDA–GMM across varying SNRs using the Friedman test.

To evaluate the decomposition of overlapping spikes, we calculated cluster-wise precision, recall, and F1-score between true spikes and those from reallocated spike trains of decomposed single units. We also confirmed the correlation between the time lags of actual spikes and their estimated counterparts. As matching the stochastic ISI distribution is crucial for assessment, we calculated the absolute difference between Gamma distribution parameters estimated from actual and estimated spike trains. In addition, we compared the proposed method to LDA–GMM by measuring event synchronization between single unit spike trains. Event synchronization is quantified by the number of nearly simultaneous occurrences of spike events, based on the relative timings of events within the time series, such as local maxima. We implemented this metric with Matlab toolbox for Event Synchronization, which is available at the [33, 34].

For real neural spikes, we have limited access to ground truth data. We thus perform a qualitative comparison by examining the change in ISI distribution of overlapping spikes. In the verification of the real neural data, we estimated and compared the Gamma distribution parameters of the ISI distribution to the results of LDA–GMM, by referencing the refractory periods of neural spike trains.