Introduction

Over the past 2 decades, neural interfacing devices have evolved from microfabricated microelectrode arrays to self-sufficient fully-implantable microsystems with bidirectional wireless connection to the external world1,2,3. Applications of brain-implantable microsystems are extended from basic and clinical neuroscience to neuro-prostheses, therapeutic applications, research and development in cognitive sciences, and brain-machine interfacing (BMI). What uniquely distinguishes intra-cortical neural recording from non-invasive recording/imaging approaches such as electro-encephalography (EEG), Magnetic- Resonance Imaging (MRI), and Near-Infrared Spectroscopy (NIRS) is the much higher temporal and spatial resolutions it offers. As the recording capacity of fully-implantable neural interfacing devices exceeded tens of concurrent channels in mid 2000s, the necessity of efficient on-implant data compression/reduction became inevitable4,5,6,7,8,9,10,11. This was mainly because of the restrictions in both the wireless transmission bandwidth and the power budget allocated for the telemetry of the recorded neuronal activities off the implant. In applications where action potentials are taken as the signal component conveying the key information, extensive data reduction is achieved through detecting and extracting action potentials (a.k.a. neural spikes) and discarding the background noise in between12. To further enhance the data compression and at the same time preserve individual spikes waveshapes, a wide spectrum of temporal and spatial compression techniques has been proposed. The temporal compression techniques, in general, reduce the volume of data associated with consecutive samples in a time window on a given recording channel. As examples of this class of spike compression techniques, one can name transform-based signal processing methods such as the discrete wavelet transform (DWT) with Daubechies4, Symlet4, and Haar basis functions13,14,15, the Walsh-Hadamard Transform (WHT)16, the discrete cosine transform (DCT)17, and the Discrete Tchebychef Transform (DTT)18. On the other hand, the spatial compression techniques take signal samples on multiple channels of concurrently-recorded neural signals, and look into the patterns existing in those parallel data streams. Compressive sensing methods19,20, inpainting-based compression21, the whitening transform22,23, and the MBED technique24 are among the effective special compression methods reported in the literature. Recently, salient progress has been reported in the microfabrication of neural recording microelectrode arrays with a channel count of around 1000 to 1600 recording channels5,6,7,8,9,10. With the introduction of such highly-dense microelectrode arrays, the bottleneck in the development of high-density neural recording brain implants shifts to the handling and subsequently wireless transmission of the huge amount of recorded data. Live streaming of such a high volume of concurrent neuronal data via a wireless link with limited bandwidth requires innovative spike compression ideas, which are saliently more efficient than the existing methods.

Methods

In this work, a spike compression method dedicated to high-density brain-implantable microsystems is proposed. As a preliminary step, action potentials are extracted and the rest of the signal (which is the inter-spike background noise) is discarded. This results in considerable data reduction, especially for low to moderate spike firing rates. As shown in Fig. 1, it is assumed that the neural recording system under study comprises an implantable module and an external module, which are linked through wireless connection. It is assumed that the recorded neural signal is first preconditioned (i.e., amplified and bandpass filtered to remove local field potentials) in the analog domain and properly converted to digital. The signal then goes through a spike detection and extraction process, and the extracted spikes are delivered to a neural signal processor for spike compression based on the method proposed in this work In the proposed method, a neural spike is represented using the start and end points as well as global and local extremum points of its waveshape, which are all referred to as salient samples, hereafter. For a given spike, first, the salient samples are extracted on the implant module, which can be taken as a ‘selective downsampling’ step. Then, the key attributes of the salient samples are transferred to the external module. On the external side, some predefined functions are fit to the received salient samples in order to reconstruct the spike waveshape. Through sending the timing (sample index) and amplitudes of only the salient samples (rather than telemetering the amplitudes of all the spike samples) from the implant to the external module, notable data reduction is achieved.

Fig. 1
figure 1

Functional diagram for the implementation of the proposed data reduction approach on a high-density brain implant.

Full size image

The proposed method offers three main advantages: (a) the volume of the data transmitted off the implant is substantially reduced, (b) the reconstruction of the spike waveshapes using predefined smooth curves allows for the removal of amplitude noise contamination of the recorded spikes, and (c) the low complexity of the computations on the implant side, and shifting the major part of the computations to the external side are the key aspects of the proposed idea, which make it a proper candidate for hardware-efficient, on-implant, neural signal compression.

In general, when downsampling a signal that is already sampled according to the Nyquist criterion25, the preservation of the key signal information becomes a critical concern. In the idea proposed in this work, part of the signal information is conveyed by the fitting function(s) used for spike waveshape reconstruction. If properly chosen, the curvature of the fitted reconstructive functions will retrieve part of the inter-sample variations that are lost as a result of downsampling.

There are two key issues that contribute to defining the details of the proposed spike compression approach: (a) The typical spike waveshape is too complicated to be modeled using ordinary algebraic or trigonometric functions, and (b) spike waveshapes vary in so many morphological aspects from one neuron to another. Therefore, as illustrated in Fig. 2a, spike waveshapes are broken into multiple primitive monotonic segments, PiPj, with the junction points, Pi and Pj, being salient samples. Attributes of the salient samples are then framed and transmitted to the external side of the system. On the external side, each segment is fitted by a smooth curve, which is formulated using a primitive fitting function, yij. Finally, the waveshape of the received spike is reconstructed by concatenating the retrieved spike segments.

Fig. 2: Illustration of the proposed idea.
figure 2

a Waveshape of a typical neural spike (dotted line), the concept of salient samples (circles), and the fitting function used to model the segment waveshape (solid line) (b) sample triplets, triplet representative samples (TRS), and TRS slopes on a typical spike waveshape. c Some of the waveshape types (with and without zero slopes at either end) that can be generated using a third-degree fitting function (d) an extremum point (Pi) at which the slope cannot (left) or can (right) be approximated to zero (e) Adding a ‘rivet sample’ to a segment to substantially enhance fitting accuracy.

Full size image

Finding salient samples

Salient samples in a neural spike were defined as the start and end of the spike, and the samples corresponding to the extremum points of the spike waveform. Finding an extremum point is done by detecting a change in the sign of the slope of the spike. To identify the extremum points with rather low computational complexity, every three consecutive samples of the spike under study is taken as a sample triplet as illustrated in Fig. 2b. Each sample triplet is represented by a triplet representative sample (TRS). The amplitude of the TRS is basically the average of the triplet samples:

$${trs}(m)=[s(3m)+s(3m+1)+s(3m+2)]/3$$
(1)

in which s(i) is the i-th sample of the neural spike, and m is the index of both the sample triplet under study and the associated TRS amplitude. A change in the sign of the slope of the TRS-spike at index m is an indication of the existence of an extremum point there. Comparison of the amplitudes of the three samples in that triplet determines which one is the extremum point. When trs(m)-trs(m-1) > 0 and trs(m + 1)-trs(m)<0, the m-th sample triplet contains a maximum point, and that is the sample with the largest amplitude in that triplet. Otherwise, a change of slope sign from negative to positive indicates that there is a minimum point in the m-th sample triplet, which is the sample with the smallest amplitude in that triplet. It is worth noting that, to further simplify the computational complexity of finding spike extremum points, the divide-by-3 is omitted from (1) as it has no effect on locating extremum points.

Primitive fitting functions

The fitting functions used in the proposed approach are chosen from a library of polynomial functions. In this work, it is suggested that each segment is interpolated primarily using a third-degree polynomial function expressed as:

$$y\left({kT}\right)={C}_{3}{({kT})}^{3}+{C}_{2}{({kT})}^{2}+{C}_{1}({kT})+{C}_{0},$$
(2)

where k is the sample index and T is the sampling time of the digitized signal. To simplify the discrete-time equations, as it is common, the sampling time is taken equal to 1 (T = 1) hereafter, and the signals are represented in terms of the sample index k. The coefficients C0 ~ C3 in (2) can, therefore, be determined according to the attributes of the start and end samples of the segment, Pi (i, Ai) and Pj (j, Aj), as well as the slopes of the signal at those samples. It is worth mentioning that, to simplify calculations on the implant side, instead of sending the aforementioned slopes, what is reported to the external setup is the ‘forward difference’ between the amplitudes of the start sample and the sample next to it (i.e., ∆F,i = Ai+1 – Ai), and the ‘backward difference’ between the end sample and the sample prior to it (i.e., ∆B,j = Aj – Aj-1).

Knowing that a typical spike has a continuous and derivable waveform, it can be said that the slope of spikes at extremum points is always zero. Assuming that a spike has a gentle rise from the baseline level at the beginning and also gradually approaches the baseline in the end, the slopes at both ends of spikes can be taken equal to zero as well. This means that the slopes at the start and end of each and every segment of a typical spike are always equal to zero. In this situation (shown as the left-most segment in Fig. 2c) the third-degree fitting function is simplified to:

$$y\left({kT}\right)={C}_{3}{({kT})}^{3}+{C}_{2}{({kT})}^{2},$$
(3)

in which

$${C}_{3}=\frac{2{A}_{{ji}}}{{(\Delta T)}^{3}}\& {C}_{2}=-\frac{3{A}_{{ji}}}{{(\Delta T)}^{2}}$$
(4)

where \({A}_{{ji}}={A}_{j}-{A}_{i}\) and \(\Delta =j-i\). However, even though the slopes at extremum points are always ‘mathematically’ equal to zero in continuous time, this might not hold true in the discrete time (i.e., for a sampled spike). The situations where the slope at an extremum point in discrete time can or cannot be taken equal to zero are shown in Fig. 2d. Specifically speaking, when a given segment, PiPj, spans a rather large amplitude during a short period of time, (the absolute value of) the slope of the spike when going from the segment start sample to the sample immediately next to it (i.e., from Pi to Pi+1) or from the segment end sample back to the sample immediately prior to it (i.e., from Pj back to Pj-1) is likely to be considerably greater than zero. The criterion for a ‘large segment amplitude span’ is when it is greater than half the full-scale amplitude range, VFS (i.e., |Aji | >VFS/2). In this case, in order for the suggested fitting function to fit such a spike segment more accurately especially at both ends, the original 4-term polynomial of (2) is used in full.

The curve fitted to a spike segment lies exactly on the original segment at both ends (i.e., at Pi and Pj), and might be to some extent off elsewhere. Traditionally, the degradation of the spike waveshape caused by a compression-reconstruction process is assessed by the reconstruction error. This error is usually quantified as the root-mean-square of the sample-to-sample difference between the original spike (\({{{\boldsymbol{s}}}}\)) and the associated reconstructed spike \(({{{\boldsymbol{r}}}})\)15,26:

$$e=\sqrt{\frac{1}{K}{\sum }_{k=1}^{K}{[{{{\boldsymbol{r}}}}\left(k\right)-{{{\boldsymbol{s}}}}\left(k\right)]}^{2}}$$
(5)

This error is also normalized to the effective value of the noiseless version of the same spike (\({{{{\boldsymbol{s}}}}}^{{{{\boldsymbol{{{\dagger}}} }}}}\)), and referred to as the Normalized RMS of reconstruction Error, NRMSE, formulated as:

$${NRMSE} \% =\frac{e}{\sqrt{\frac{1}{K}{\sum }_{k=1}^{K}{[{{{{\boldsymbol{s}}}}}^{{{{\boldsymbol{{{\dagger}}} }}}}\left(k\right)]}^{2}}}$$
(6)

Rivet samples

As a result of fitting a curve to the spike segment of interest, the reconstruction error is equal to zero at both ends of spike segments and non-zero otherwise. Even though for most intra-cortical neural spikes with typical waveshapes the 3rd-order polynomials of (2) or (3) fit all spike segments sufficiently accurately, there are still exceptions where the reconstruction error is unacceptably high. As a remedy, as illustrated in Fig. 2e, the fitted curve is rivetted to the original spike segment at a sample in the middle of the segment, almost where the fitting error peaks. The fitting error is, consequently, zeroed at the rivet sample, PR, and the maximum error is therefore meaningfully lowered. In this situation, to allow for a more accurate curve fitting, a five-term 4th-order polynomial function is used:

$$y\left({kT}\right)={C}_{4}{\left({kT}\right)}^{4}+{C}_{3}{\left({kT}\right)}^{3}+{C}_{2}{\left({kT}\right)}^{2}+{C}_{1}\left({kT}\right)+{C}_{0}$$
(7)

The coefficients C4-C0 in (5) are calculated on the external side of the system using the attributes of the segment start and end samples and the rivet sample. The location of the rivet sample cannot be determined on the implant side according to where the fitting error is maximized as it would require doing the segment reconstruction computations on the implant side. Therefore, always the 6th sample after the segment start sample is experientially taken as the rivet sample. This way, there will be no need for reporting the timing of rivet samples. It is interesting to note that the approach proposed in this work does not rely on any specific spike pattern nor on identified spike classes. Therefore, the proposed compression method does not care whether the detected amplitude variation is for a single neural spike or it is the overlap of two different spikes.

The neural spikes studied in this work are of specific waveshapes and contaminated with background noise. As a result of the spike compression technique introduced in this work, not only the spike waveshape is altered by curve fitting to some extent, most of the spikes are also subject to considerable noise reduction. Therefore, before continuing with the details of the proposed spike compression approach, it is of crucial importance to introduce more appropriate measures for the assessment of the quality of the signals (neural spikes) before and after they are processed.

Signal quality assessment

Introduction of a proper measure for reconstruction performance assessment requires a brief neuroscientific background, which is provided as follows:

In neuroscience, it is believed that the action potentials generated by a given neuron are almost the same in waveshape at the origin27. When they are recorded extracellularly, however, the background noise causes slight amplitude differences, which turn those identical spikes into a ‘class’ of spikes. The ‘within-class variability’ caused by the background noise is usually ignored and all the spikes in the class under study are considered of the same origin when it comes to interpreting the information conveyed by those spikes. Suppose that the neural spike extra-cellularly acquired by a brain-implantable neural recording device, \({s}_{i}\left(k\right)\), is the i-th firing of the neuron under study, where k is the sample index, ranging from 1 to K. This spike, as explained above, is assumed to be a noiseless spike at the origin, \({s}^{{{\dagger}} }\left(k\right)\), contaminated by the i-th observation of a random noise process, \({n}_{i}\left(k\right)\):

$${s}_{i}\left(k\right)={s}^{{{\dagger}} }\left(k\right)+{n}_{i}\left(k\right)$$
(8)

The noise is assumed to have a mean value of zero and a standard deviation of \({\sigma }_{N}\). The dissimilarity between the recorded noisy spike and the noiseless spike, which is referred to as their dissimilitude hereafter, is defined as:

$${d}_{i}=\sqrt{\frac{1}{K}{\sum }_{k=1}^{K}{[{s}_{i}\left(k\right)-{s}^{{{\dagger}} }\left(k\right)]}^{2}}={\sigma }_{N}$$
(9)

As illustrated in Fig. 3, as long as the noise standard deviation remains the same, the locus of all the spikes of the same origin recorded under the same conditions is a circle of radius \({\sigma }_{N}\) with the associated noiseless spike, \({s}^{{{\dagger}} }\), at the center. In other words, despite their sample-to-sample, noise-induced amplitude differences, all such spikes have the same dissimilitude with the noiseless spike at the center. Let us assume that a spike processing task (smart downsampling followed by proper curve fitting in our case) converts recorded noisy spikes, s, to reconstructed spikes, r. As a result, the dissimilitude of the spike under study with the associated noiseless spike at the center changes from \(d(={\sigma }_{N})\) to \(\hat{d}\). The degradation of the spike waveshape caused by the compression-reconstruction process, which was already quantified in (5) by the reconstruction error and can also be expressed in terms of the dissimilitudes of the original and reconstructed spikes with the noiseless spike (i.e., \(d\) and \(\hat{d}\), respectively) as:

$$e=\sqrt{|{(\hat{d})}^{2}-{(d)}^{2}|}$$
(10)
Fig. 3: Assessment of signal quality.
figure 3

a The traditional reconstruction error, e, versus the proposed ‘added dissimilitude’, ε. The former measures the dissimilarity of the reconstructed spike to the original noisy spike, and the latter quantifies the impact of spike processing on the dissimilarity of the spike under study with respect to the associated noiseless spike (as a reasonable reference). b A spike processing task, in general, might be able to reduce the noise content of a recorded spike. In this case, while the traditional reconstruction error is always positive and misleadingly reports signal quality degradation, the added dissimilitude will be negative (ε<0) indicating signal quality enhancement (i.e., denoising).

Full size image

The important point here is that while the reconstructed spike is away from the original noisy spike by \(e\), what matters is still the resemblance of the reconstructed spike to the associated noiseless spike, not to the noisy spike! Therefore, to assess the real impact of the compression-reconstruction process on the waveshape of the spike being processed, it is proposed to measure the change it causes in the dissimilitude between the that spike and the noiseless spike. Referred to as the relative added dissimilitude, this is simply the relative difference between the dissimilitudes of the original and reconstructed spikes with the noiseless spike, denoted in Fig. 3a as \(\varepsilon\):

$$\varepsilon \% =[(\hat{d}-d)/d]\times 100 \%$$
(11)

There is yet another reason for the superiority of the proposed ‘added dissimilitude’, \(\varepsilon\), to the traditional ‘reconstruction error’, \(e\) (given in (10)), which is highlighted when the signal processing task to some extent denoises the spike (as it is the case in this work). According to the definition provided precedingly, denoising reduces the dissimilitude between the spike being processed and the noiseless spike, and the added dissimilitude, therefore, takes on negative values. This is while the traditional measure (\(e\)) is blind to denoising and misleadingly reports an always positive reconstruction error in such scenarios.

It is true that the proposed spike compression technique fits smooth curves to the salient samples of the spikes being compressed, and that provides the opportunity of reducing the unwanted amplitude fluctuations caused by the noise that is mixed with the signal upon recording. The curve fitting error together with the random nature of the added noise, however, make it difficult to claim that our spike compression always brings the recorded spikes closer to the associated noiseless spike than they were prior to compression. Therefore, the spike denoising property of the proposed idea is statistically studied for a class of M recorded noisy spikes of the same origin. As Fig. 3b illustrates, let us assume that some of the reconstructed spikes fall inside the dotted circle (meaning that they have come closer to the associated noiseless spike), and some others are not. To gain an overall sense about the denoising property for the entire spike class under study, the Net Added Dissimilitude (NAD) is defined as:

$${NAD} \% =\frac{1}{M}{\sum }_{m=1}^{M}{\varepsilon }_{m}\times 100 \%$$
(12)

where M is the number of spikes in the class under study, and εm is the added dissimilitude for the m-th spike. A negative NAD% is an indication of the reconstructed spike ensemble being of smaller average distance to the noiseless spike than the original noisy spikes are. In this case, it can be claimed that, overall, spike denoising has been accomplished. It is worth noting that the neither the concept of the NAD metric proposed in this paper, nor its definition and applicability is specific to neural spikes and their properties. Therefore, this measure can be used for any similar situation in other applications as well.

Impact of noise

According to the above discussions, so long as the spike under study is noiseless, selective downsampling and subsequently fitting smooth curves to the resulting salient samples works quite well. In reality, however, the non-negligible amplitude noise that is unavoidably mixed with the spike inevitably affects the idea proposed in this work. This section studies two major effects of the existence of noise: (a) extremum point displacement, and (b) the impact of noise-related salient sample amplitude fluctuations on the fitted functions.

Assuming that the locations of the start and end samples are not affected by noise, our discussion here is focused on how the introduction of amplitude noise affects the locations of extremum points. Perhaps the most critical impact of the background noise is the displacement of extremum samples for the spike under study. The basis for the analysis of this displacement is presented in Fig. 4. A part of a noiseless spike containing an extremum point (point G) is shown as a solid line on the right-hand side of Fig. 4a. Let us assume that each and every sample of this spike is subject to Gaussian amplitude noise with zero mean value (\({\mu }_{N}=0\)) and standard deviation of \({\sigma }_{N}\). As shown on the left-hand side of Fig. 4a, the amplitudes of the extremum sample G and all other samples including the neighboring samples F and H, will therefore have Gaussian distributions with the same standard deviation (i.e., \({\sigma }_{F}={\sigma }_{G}={\sigma }_{H}={\sigma }_{N}\)) around the associated noiseless amplitudes AF, AG, and AH, respectively (i.e., \({\mu }_{F}={A}_{F}\), \({\mu }_{G}={A}_{G}\), and \({\mu }_{H}={A}_{H}\)). As illustrated in Fig. 4a, such independent noise-induced amplitude variations can slightly change the spike waveshape (from the solid curve to the dashed curve). As a result, there is a chance one of the neighboring samples F or H becomes the extremum point rather than the sample G. The probability of the displacement of the extremum point from G(\({T}_{G},{A}_{G}\)) to a neighboring point, say H(\({T}_{H},{A}_{H}\)), can be calculated by finding the probability of a sign change for the amplitude difference: \({A}_{D}={A}_{G}-{A}_{H}\). The difference between two random variables each being of Gaussian distribution (e.g., \({P}_{G}({\mu }_{G},\,{\sigma }_{G})\) and \({P}_{H}({\mu }_{H},\,{\sigma }_{H})\) in our case) is of Gaussian distribution itself28, \({P}_{D}({\mu }_{D},{\sigma }_{D})\), where \({\mu }_{D}={\mu }_{G}-{\mu }_{H}={A}_{G}-{A}_{H}\) and \({\sigma }_{D}=\sqrt{{\sigma }_{G}^{2}+{\sigma }_{H}^{2}}=\sqrt{2}{\sigma }_{N}\). Assuming that the sample G is a maximum point, as shown in Fig. 4a, the difference \({A}_{D}\) is positive for the noise-free scenario, and the probability of a noise-induced extremum point displacement is equal to the area under the PDF curve for \({A}_{D} < 0\). Shown in Fig. 4b as the hashed area, this is the probability of the sample ‘H’ becoming the extremum point (rather than G), and called the probability of extremum point displacement (PXD), hereafter:

$${P}_{{XD}}= \int _{-\infty }^{0}D\left({\mu }_{D},{\sigma }_{D}\right){{{\rm{d}}}}x\\ = \frac{1}{2}[1-{erf}(\frac{{A}_{D}}{2{\sigma }_{N}})]$$
(13)
Fig. 4: Displacement of an extremum sample as the consequence of adding noise to the spike amplitude.
figure 4

a The extremum sample and the neighboring samples are all subject to amplitude noise. As a result, the location of the extremum sample might accordingly change. b The probability of extremum points displacement (PXD) can be calculated using the PDF of the amplitude difference between the extremum sample and the neighboring samples. c The impact of the spike slopes around extremum samples on the PDFs of extremum sample displacement in a typical spike. d The impact of amplitude noise at both ends of a segment on the fitting functions used to reconstruct the segment waveshape.

Full size image

In (13), erf(.) is the Gaussian distribution error function, and the probability of extremum point displacement is noticeable (i.e., PXD > 0.017) if AD is smaller than \({3\sigma }_{N}\). The chance of extremum point relocation is, therefore, contingent on both the extent of noise contamination superposed on the spike and the slope of the spike around the extremum point under study. Figure 4c illustrates the PDFs of displacements for two extremum points in a typical spike. So long as the noise amplitude around an extremum point is smaller than \({3\sigma }_{N}\), the displacement of that point is probable with a Gaussian distribution. The higher (the absolute value of) the spike slope is around a given extremum point, the narrower the bell-shape displacement PDF will be, meaning that notable displacement of that point will be less probable.

Aside from the displacement of extremum points in time, noise-induced fluctuations in the amplitude of salient samples also cause changes in the fitted functions. To study this effect, as illustrated in Fig. 4d, let us assume that the amplitudes of the salient samples at both ends of a given spike segment, Ai and Aj, are contaminated with uncorrelated random Gaussian noise with a mean value of zero and standard deviation of \({\sigma }_{N}\) (\(N\left(0,{\sigma }_{N}\right)\)). Therefore, the amplitudes of the noisy salient samples are expressed as \({\widetilde{A}}_{i} \sim N\left({A}_{i},{\sigma }_{N}\right)\) and \({\widetilde{A}}_{j} \sim N\left({A}_{j},{\sigma }_{N}\right)\). Since the polynomial coefficients are functions of a linear combination of salient sample amplitudes (e.g., ref. (4)), this noise causes Gaussian noise in the coefficients of the associated fitting function. It can be shown that the coefficients of the two-term fitting function in the presence of noise are expressed as:

$${\widetilde{C}}_{3} \sim N(\frac{2}{{\Delta }^{3}}{A}_{{ji}},\frac{2\sqrt{2}{\sigma }_{N}}{{\Delta }^{3}})$$
(14)

and

$${\widetilde{C}}_{2} \sim N(\frac{-3}{{\Delta }^{2}}{A}_{{ji}},\frac{3\sqrt{2}{\sigma }_{N}}{{\Delta }^{2}})$$
(15)

where \({A}_{{ji}}={A}_{j}-{A}_{i}\) and \(\Delta =j-i\). In this case, it can be shown that the amplitude of a given reconstructed sample (say sample #k) calculated using the resulting fitting function is subject to Gaussian noise, formulated as:

$${\widetilde{R}}_{k} \sim N({R}_{k},{\sigma }_{R,k})$$
(16)

In (16), k is the sample index ranging from i to j (which are the indices for the start and end samples of the segment under study, respectively), \({R}_{k}\) is the reconstructed amplitude for sample #k prior to adding noise, expressed as:

$${R}_{k}=\left(\frac{2{k}^{3}}{{\Delta }^{3}}-\frac{3{k}^{2}}{{\Delta }^{2}}\right){A}_{{ji}}$$
(17)

and \({\sigma }_{R,k}\) is the standard deviation of amplitude variations for the reconstructed sample #k, formulated as:

$${\sigma }_{R,k}={\sigma }_{N}\sqrt{{\left(\frac{2{k}^{3}}{{\Delta }^{3}}-\frac{3{k}^{2}}{{\Delta }^{2}}+1\right)}^{2}+{\left(-\frac{2{k}^{3}}{{\Delta }^{3}}+\frac{3{k}^{2}}{{\Delta }^{2}}\right)}^{2}}$$
(18)

To analyze the effect of amplitude noise on the reconstruction error, the Normalized Reconstruction Error Power (NREP) is defined as:

$${{NREP}={NRMSE}}^{2}=\frac{\frac{1}{\Delta }{\sum }_{k=i}^{j}{({\widetilde{R}}_{k}-{\widetilde{s}}_{k})}^{2}}{{A}_{P-P}^{2}}$$
(19)

Here, \({\widetilde{s}}_{k}\) is the noisy amplitude of sample #k:

$${\widetilde{s}}_{k} \sim N\left({s}_{k},{\sigma }_{N}\right)$$
(20)

in which \({s}_{k}\) denotes the amplitude for sample #k prior to adding noise, \({\widetilde{R}}_{k}\) is the noise-induced reconstructed function at sample #k, and \({\sigma }_{N}\) is the standard deviation of added noise. Hence, the reconstruction error, \({X}_{k}\), is formulated as the difference between the noise-induced reconstructed function (\({\widetilde{R}}_{k}\)) and the noisy spike (\({\widetilde{s}}_{k}\)), and is of Gaussian distribution:

$${X}_{k}=\left({\widetilde{R}}_{k}-{\widetilde{s}}_{k}\right) \sim N\left({E}_{X,k},{\sigma }_{X,k}\right)$$
(21)

in which

$$\left\{\begin{array}{c}{E}_{X,k}={R}_{k}-{s}_{k} \hfill \\ {\sigma }_{X,k}=\sqrt{{{\sigma }_{R,k}}^{2}+{\sigma }_{N}^{2}}\end{array}\right.$$
(22)

Equation (21) can also be written as a Gaussian distribution in the standard form (i.e., with unity standard deviation):

$${X}_{k}=\sqrt{{{\sigma }_{R,k}}^{2}+{\sigma }_{N}^{2}}(N(\frac{{E}_{k}}{\sqrt{{{\sigma }_{R,k}}^{2}+{\sigma }_{N}^{2}}},1))$$
(23)

Replacing the reconstruction error (\({X}_{k}=\left({\widetilde{R}}_{k}-{\widetilde{s}}_{k}\right)\)) by the expression provided in (23), the NREP introduced by (19) can be expressed as:

$${NREP}= \, \frac{1\,}{{\Delta A}_{P-P}^{2}}{\sum }_{k=i}^{j}{X}_{k}^{2}\\ = \, \frac{{{\sigma }_{R,k}}^{2}+{\sigma }_{N}^{2}\,}{{\Delta A}_{P-P}^{2}}{\sum }_{k=i}^{j}{(N(\frac{{E}_{k}}{\sqrt{{{\sigma }_{R,k}}^{2}+{\sigma }_{N}^{2}}},1))}^{2}$$
(24)

which is of a noncentral Chi-squared distribution28,29. Noncentral chi-squared distribution has two parameters: the number of \({X}_{k}\) terms in the summation in (24) called the degrees of freedom for the distribution (i.e., \(n=j-i+1\)), and \(\lambda\) referred to as the non-centrality parameter of the distribution written as:

$$\lambda ={\sum }_{k=i}^{j}{(\frac{{E}_{k}}{\sqrt{{{\sigma }_{R,k}}^{2}+{\sigma }_{N}^{2}}})}^{2}$$
(25)

It can be shown that the total mean value and variance of the \({NREP}\) can be expressed as \(n+\lambda\) and \(2(n+2\lambda )\), respectively.

Number of single-channel engines

It is believed in neuroscience that spontaneous neuronal activities occur with the Poisson distribution28. This is taken as the computational basis for the calculation of the number of spike compression engines (N) in the 128-channel spike compressor. The probability of having a given number of events occurred in a fixed interval of time with the Poisson distribution is formulated as:

$$P\left(X=x\right)=\frac{{\lambda }^{x}{e}^{-\lambda }}{x!}$$
(26)

where λ is the average number of events, and x is the number of events. The average number of events (occurrence of neural spikes in our case), λ, is calculated in our case as:

$$\lambda ={FR}\times {T}_{{SPK}}\times {N}_{{ch}}$$
(27)

where FR is the spike firing rate, TSPK is the average time course of a typical spike, and Nch is the total number of channels. For a neural signal with average spike firing rate of 100 spike/s and typical spike duration of 2 ms, the average number of events on Nch = 128 channels equals 25.6 spikes at any given time.

Hardware design and implementation

To realize a multi-channel neural spike compressor based on the proposed idea, the associated hardware needs to comply with both signal processing requirements and hardware implementation restrictions. From the perspective of signal processing, the compressor circuit is expected to do the processing in the real time with acceptable precision. On the other hand, from a physical implementation standpoint, the hardware is required to be small in physical size and efficient in power consumption.

In order to be embedded in a high-density neural recording microsystem, a single-channel spike compression engine is designed based on the proposed approach. Figure 5 presents a functional block diagram and timing diagram of this engine. Incoming samples of the input neural signal are grouped into consecutive non-overlapping triplets. The sum (average) of the amplitudes of the three samples in each triplet form a low-pass filtered version of the spike (shown as TRS), which is of less noise contamination compared to the original spike. The slope of the TRS signal is measured in the discrete time and is taken as the basis to identify the extremum samples of the spike. Amplitudes and timings of the extremum samples are then stored in a local memory referred to as the ‘Salient Sample Attribute Storage’. As mentioned before, the part of the spike that lies between every two consecutive salient sample are taken as a spike segment. If the segment spans a rather wide amplitude range (half the full-scale range in this work), slopes at both ends of the segment are calculated and stored in the Salient Sample Attribute Storage. After the completion of the spike, attributes of all the salient samples are framed in the form of a serial data packet. Each data packet conveys the salient sample information for all the segments of the spike being compressed.

Fig. 5: The single-channel spike compression engine implementing the proposed spike compression approach.
figure 5

a Functional block diagram, b Simplified timing diagram, c Simplified block diagram of the 128-channel spike compressor.

Full size image

Based on the proposed spike compression technique, a 128-channel spike compressor is designed. Figure 5c shows a simplified block diagram for this spike compressor, which comprises N single-channel spike compressors of the type shown in Fig. 5. Prior to this circuit, neural signals on all the 128 channels are preconditioned in the analog domain, converted to digital, and then undergo a spike detection process. What is delivered to the circuit is a time-multiplexed stream of extracted neural spikes in the real time. The channels having spikes on them are referred to as ‘active channels’. An active channel router receives the isolated spikes, and directs them to one of the single-channel spike compressors. In the end, one shared data framing block collects and frames the information prepared by the spike compressor engines. In the case of having concurrent spikes on different channels, they are treated on a first come, first packed basis.

According to the analysis provided above, for λ = 25.6, the probability of having more than 32 concurrent active channels (X = 32) is as low as 12.36%. Therefore, 32 salient sample detectors are envisioned in the 128-channel compressor being designed. It should be added that each data packet generated by the 128-channel compressor carries the attributes of all the segments of only one spike on a given channel.

Real-time operation of the proposed spike compressor is of crucial importance in applications such as prosthetic devices and brain-machine interfaces30,31. The fact that the single-channel salient sample detector functions in real time when operated at a clock rate equal to the basic front-end sampling rate (i.e., 25 k Sample per Sec.) is translated into extremely low dynamic power consumption, especially when the system is designed for high channel counts. The only parts of the 128-channel spike compressor that is operated at a high clock rate are the active channel router and the data framing block. To guarantee live streaming of the compressed neuronal activities to the outside, these blocks receive a clock rate 128 times higher than the basic sampling clock rate.

The attributes of the salient and rivet samples acquired by the 128-channel spike compressor are framed in the form of serial data packets. Each packet starts with a start pulse and ends with a stop pulse for synchronization purposes. To allow for the live streaming of spike information, the data framing block is operated at a 128 times higher rate than the basic sampling clock.

Aside from signal- and system-level considerations, the efficient design of the hardware realizing the proposed compression technique is of crucial importance in achieving implant-appropriate performance. The efforts put in this direction can be summarized as follows:

  • At the architecture level, multiple ‘Salient Sample Detector’ blocks are available in parallel in order to process incoming spikes on active channels concurrently. This allows for the real-time operation of the multi-channel spike compressor with no enhancement in the clock rate, hence avoiding unnecessary penalty in dynamic power consumption. Moreover, the hardware used for ‘Averaging’ and ‘Slope Measurement’ in the block diagram of Fig. 5a are shared among all the 32 Salient Sample Detectors. Given that the Salient Sample Detectors are operated sequentially with a duty cycle of 1:32, this makes the designed processor efficient in silicon area (13%) with no speed or power consumption penalty.

  • At the building block level, area-/power-efficient computational building blocks are used to realize the processing tasks. As examples, multiple additions in the Low-Pass Filtering blocks (for the calculation of the amplitude of triplet average samples, TRSs) are simultaneously performed using carry-save adders, gray-code counters are used (for purposes such as sample time stamping and Salient Sample Detector address generation) to lower the consumed power.

  • At the signal level, the entire algorithm proposed in this work is realized with no need for multipliers, resulting in savings in power and area. To avoid extra power and area consumption, the length of the data conveying amplitudes is always truncated to 8 bits. Instead of calculating spike slopes, the differences between the amplitudes of consecutive samples are calculated and transmitted off the implant module. The temporal spacing between salient samples is transmitted to the outside rather than their absolute time stamps. This results in 9% enhancement in the compression rate. The Compression Rate (CR) in this work is defined as the ratio of the original number of bits (prior to compression) to the total number of bits conveying the information of the same spike after applying the proposed compression technique.

  • At the protocol level, outgoing data packets convey the minimum possible information, using which the data required for curve fitting and spike reconstruction are retrieved. For instance, salient sample type (start sample, end sample, or spike extrema), salient sample absolute time stamps, the polynomial orders for spike segments are all retrieved on the external side.

Once the attributes of salient samples of a given spike arrive at the external side of the system, the associated waveshape of that spike is reconstructed using a rule-based procedure. Figure 6 shows the flowchart of the spike reconstruction procedure. The procedure first locates the salient samples according to the amplitudes and timings received from the implant module. It then identifies the number of spike segments, and calculates the associated heights. Based on the height of each spike segment, the degree of the primitive fitting function is decided upon. If the segment height is greater than half the full-scale range, VFS, the 4th-degree polynomial of (7) is used as the fitting function; otherwise, the two-term 3rd-degree polynomial given in (3) is employed. In the former case, the polynomial coefficients are calculated using the timings, amplitudes, and slopes of the salient samples and the amplitude of the rivet sample accompanying the segment information. In the latter case, coefficients of the polynomial are simply computed using the timings and amplitudes of the segment start and end. Finally, the segments are concatenated and a unified spike is derived.

Fig. 6
figure 6

Flowchart of the rule-based spike waveshape reconstruction procedure.

Full size image

Results

The microfabricated neural spike compressor

The 128-channel neural spike compressor designed in this work was microfabricated in TSMC 130-nm CMOS process. Figure 7 shows the physical layout of the core circuit measuring 350 μm by 1050 μm (part (a)), the photograph of the fabricated and wire-bonded chip (part (b)), and the photograph and block representation of the experimental setup (parts (c)&(d)). As presented in Fig. 7d, the prototyped ASIC board is configured and receives the input data using an FPGA-based interface board, and delivers the outgoing serial data packets to an external host computer through a Data Acquisition (DAQ) box.

Fig. 7: The microfabricated 128-channel spike compressor.
figure 7

a physical layout of the core circuit (b) photograph of the fabricated and wire-bonded chip, c the experimental setup, and (d) the block representation of the experimental setup.

Full size image

Data set to verify the proposed spike compression idea

To verify the efficacy and evaluate the performance of the proposed spike compression technique, a dataset is composed using 8 different real spike waveshapes acquired from in-vivo, extra-cellular recordings32,33,34. Figure 8a presents the 8 spikes used in this work with a typical duration of 2 ms, and sampled with a resolution of 8 bits at 25 kSample per Sec. These spikes are first completely denoised by averaging 100 recorded real spikes from each spike class. To exhibit the success of the spike compression technique proposed in this work, in the spikes used for the tests, the number of extremum points range from 1 to 6, and the order and amplitudes and also timings of the extremum points are all different. To study the impact of noise on the performance of the proposed technique, the 8 noiseless spikes are then contaminated by white Gaussian noise, resulting in a dataset of 4800 spikes with a rather wide signal-to-noise (SNR) range, extended from 5 dB to 25 dB.

Fig. 8: Performance of the proposed technique for prerecorded neural spikes.
figure 8

a Eight different spike waveshapes acquired through in-vivo recording31,32,33. These spikes are used to verify the functionality and assess the performance of the proposed spike compressor. b Spike segmentation and curve fitting for three sample spikes without noise contamination. c Effect of adding rivet samples on curve fitting. d Reconstruction of Spike #1 for three SNR values: 20 dB, 25 dB, and 30 dB. e The dissimilitude diagram for a class of spikes with waveshape #5 (according to Fig. 8a) at SNR = 15 dB. f The NADs for all the 8 spike waveshapes as a function of the SNR. g The normalized root-mean-square of reconstruction error (NRMSE%) for all the 8 spike waveshapes as a function of the SNR.

Full size image

Figure 8b shows the performance of the segmentation and curve fitting for three sample spikes. In this test, only noiseless spikes are used to report merely the contribution of the curve fitting error in the reconstruction of spikes.

It is worth mentioning that the three spike waveshapes reported in Fig. 8b are so chosen to represent the typical (spikes #1&2) and maximum (spike #8) curve fitting errors. To demonstrate the impact of the rivet point in reducing the fitting error, the graph in Fig. 8c compares the curve fitting errors for all the 8 spike waveshapes with and without rivet points. Averaged on all the noiseless spikes in the dataset, the overall curve fitting error calculated based on normalized sample-to-sample error26 reduces from 5.98% to 1.8% as a result of using rivet samples. The price of this substantial improvement in the accuracy of the curve fitting is paid by sacrificing a fraction of the compression rate.

To be able to fairly compare the compression rate in this work with those reported in the literature, the compression rates need to be normalized to the firing rate. This gives a measure introduced in the literature as the True Compression Rate (TCR)15. In this work, the improved fitting accuracy (as a result of adding rivet samples) is achieved at the price of degradation in the compression rate (CR) from 390 down to 272 at a firing rate of 8 spikes/sec. Normalized to the firing rate, the true compression rate is 2176 (\(=272\times 8\)), which outperforms all other works so far reported in the literature (ref. Table 1).

Table 1 Performance summary of the proposed neural signal processer in data reduction in comparison with other works
Full size table

To study the denoising property of the proposed spike compression method, for each one of the waveshapes shown in Fig. 8a, classes of 100 prerecorded extracellular neural spikes with different signal-to-noise ratios (SNRs) are formed. For each spike waveshape, the SNR ranges from 5 dB to 25 dB. The results of waveshape reconstruction for Spike#1 for SNR = 20 dB, 25 dB, and 30 dB are presented in Fig. 8d. Figure 8e shows the ‘dissimilitude diagram’ reflecting the results of this study for Spike#5 at \({SNR}=15\). Both horizontal and vertical axes are normalized to the noise standard deviation, with all the 100 noisy spikes aligned on the same spot shown as a small blue filled circle, and the reconstructed spikes shown as small white filled circles. Interestingly, 85% of the reconstructed spikes are located inside the red circle, indicating negative added dissimilitudes. In other words, in 85% of cases, the compression-reconstruction approach exhibits the side benefit of spike denoising as well. The NAD for this ensemble of spikes is -27%, meaning that for this spike class overall, the proposed spike compression technique not only introduces no additional error, but it also even notably reduces the overall noise contamination of the whole spike class (by 27%). Figure 8f shows the NADs calculated for all the spike waveshapes in the dataset at a wide SNR range, from 5 dB to 25 dB. The NADs that fall below the zero level (shown using the dashed line) indicate denoising. According to this plot, the noise reduction property of the proposed technique is meaningful at low SNR scenarios (i.e., below 15 dB), which is highly desirable. As the SNR increases, the radius of the circle in the dissimilitude diagram shrinks down, and that causes an increasingly higher number of reconstructed spikes lie outside the circle. It is true that this is an indication of degradation in noise reduction, but it evidently happens at high SNRs (i.e., above 15 dB) where the noise content is negligible and failing to reduce it is not considered as a weakness. The positive NAD in such SNR scenarios is predominantly due to the fitting error rather than the noise content of the reconstructed spike. This error is, however, still so small that keeps the NRMSE% well below the commonly acceptable upper limit of 10 as presented in Fig. 8g. To conclude from the plots presented in Fig. 8f, g, at low SNRs, the proposed compression technique efficiently denoises the spikes. The high NRMSE for such SNR scenarios is therefore misleading and should not be taken as a drawback. This is because the NRMSE merely measures the distance between the reconstructed spike and its original noisy form and does not reflect the fact that the reconstructed spike is closer in waveshape (compared to the original noisy spike) to the associated noiseless form.

Figure 9 shows the displacement histograms (shown in blue, green, and red separately) for the three extremum points of 100 occurrences of an actual action potential with a signal-to-noise ratio (SNR) of 15 dB. It can be shown that such histograms would be of Gaussian form should the spike slopes on both sides of the extremum point were constant and equal in absolute value. In reality, as the histograms show, the profile of each one of the displacement distributions is skewed as a result of the non-constant and unequal spike slopes around the associated extremum points. Moreover, there is an inverse correspondence between the widths of the histograms and the spike slope around the associated extremum point.

Fig. 9
figure 9

Noise effect on time stamps of P2, P3, and P4 for the spike waveshape shown in gray for SNR = 15 dB.

Full size image

According to our measurements, operated at a supply voltage of 1 V and a master clock rate of 3.2 MHz, the 128-channel processor consumes a total power of 21 μW (0.164 μW per channel). The breakdowns of the overall silicon area and measured power consumed by the 128-channel spike compressor are shown in Fig. 10. As the pie charts show, as low as only 10% of the area and 26% of the power is spent on computations in the proposed algorithm, and the rest is consumed by storage elements.

Fig. 10
figure 10

Breakdowns of silicon area and measured power for the 128-channel spike compressor.

Full size image

Table 1 summarizes the design, operational, and performance specifications of the spike compressor in this work and compares them with those of other similar works appeared in the literature. Thanks to the computationally-light technique proposed in this work, our measured per-channel power consumption is notably lower (by 62%) than the best spike compressor reported so far, and the per-channel silicon area occupied by our processor is less than two thirds of the most compact work appeared in the literature. In addition to the superiority of the proposed spike compressor to other works in hardware implementation, its signal compression performance shows considerable improvement compared to the best spike compressors ever reported. As a normalized measure for signal compression, the TCR in this work is more than twice the highest TCR appeared in the literature (2176 vs. 915).

Conclusions

This paper proposes an innovative technique in spike compression, which is based on the segmentation of neural spikes, transmission of the key spike segment attributes off the neural recording implant module, and reconstruction of the spike waveshape by curve fitting. This and several other design techniques at signal and circuit levels help achieve the highest compression performance ever reported.

Moreover, it is explained in this paper that the reconstruction error traditionally used for measuring the extent of spike waveshape deformation caused by signal processing is a misleading measure. We have shown that this measure is blind to the denoising property of spike processing techniques, if any.

Instead, we have proposed ‘relative added dissimilitude’ and ‘net added dissimilitude’ as innovative measures for the assessment of the impact of signal processing on the quality of the neural spikes under study and have shown that the latter can serve as a quantitative figure for measuring spike waveshape denoising.