Deep Spike Decoder (DSD)
Authors: Emrah Adamey
Tarin Ziyaee
Nishanth Alapati
Jun Ye
Venue:  (2020)
Keywords: self-supervised
deep learning
spike sorting
EMG
sEMG
autoencoder
inductive bias

Abstract: Spike-sorting is of central importance for neuroscience research.  We introducea novel spike-sorting method comprising a deep autoencoder trained end-to-endwith a biophysical generative model, biophysically motivated priors, and a self-supervised loss function to training a deep autoencoder. The encoder infers the ac-tion potential event times for each source, while the decoder parameters representeach source’s spatiotemporal response waveform.  We evaluate this approach inthe context of real and synthetic multi-channel surface electromyography (sEMG)data, a noisy superposition of motor unit action potentials (MUAPs).  Relative toan established spike-sorting method, this autoencoder-based approach shows su-perior recovery of source waveforms and event times.  Moreover, the biophysicalnature of the loss functions facilitates interpretability and hyperparameter tuning.Overall, these results demonstrate the efficacy and motivate further developmentof self-supervised spike sorting techniques.

Under review as a conference paper at ICLR 2020
DEEP SPIKE DECODER (DSD)
Anonymous authors Paper under double-bling blind review
ABSTRACT
Spike-sorting is of central importance for neuroscience research. We introduce novel spike-sorting method comprising deep autoencoder trained end-to-end with biophysical generative model, biophysically motivated priors, and a self- supervised loss function to training deep autoencoder. The encoder infers the ac- tion potential event times for each source, while the decoder parameters represent each source's 's spatiotemporal response waveform. We evaluate this approach in the context of real and synthetic multi-channel surface electromyography (sEMG) data, a noisy superposition of motor unit action potentials (MUAPs). Relative an established spike-sorting method, this autoencoder- perior recovery of source waveforms and event times. Moreover, the biophysical nature of the loss functions facilitates interpretability and hyperparameter tuning. Overall, these results demonstrate the efficacy and motivate further development of self-supervised spike sorting techniques.
approach shows su-
INTRODUCTION
Spike sorting focuses on analyzing electrophyssilogical data, detecting occurrences of action poten- tials, and assigning them to distinct neurons. The resulting spikes-based representation is regarded as a complete, high-resolution, noise-free, and disentangled representation of the underlying neural activity. Spike sorting has been studied extensively for use with extracellular electrode (reviewed by Rey et al. 2015; Carlson & Carin 2019: Hennig et al. 2019) and surface clectromyography (sEMG; Farina et al. 2004| Holobar & Zazula 2007] Negro al. 2016). In this paper, we propose novel solution to the spike sorting problem that marries biophysical inductive biases and constraints with the power of self-supervision.
Our approach uses a deep autoencoder architecture where the encoder functions as spike sorter and the decoder implements generative model. The entire model is trained end-to-end in completely self-supervised fashion by leveraging a set of loss terms that enforce piophysicallymmotiaaie ors. In this paper. the proposed approach evaluated using both real-world and realistically simulated sEMG datasets. The results demonstrate the validity of our approach and motivate further develop- ment of the overall paradigm of unsupervised deep spike sorting.
pri-
PROBLEM FORMULATION
Here, we briefly discuss the problem formulation, and generative model giving rise to the measured sEMG signal (Merletti & Farina 2016).
A motor neuron innervates a set of muscle fibers, and together they are referred to as a motor unit (MU). As an MU emits spike, denoted by the variable s, it passes through layers of connective tissue, blood, and skin, before being picked up by the electrodes on the surface. This pattern is referred to as the
(MUAP). The received sEMG signal at the electrode, is the result of an MU spike being convolved with kernel, given by the tissue between the m 'th MU and the e'th electrode, he.m- The result of this convolution is the MUAP-waveform, which since the spike is binary is also the shape of the convolution kernel. The kernel shape is directly and physically related to the space between the MU and the electrode.
As other spikes from the same MU are transmitted, they also convolved with the same kernel h.... on the way to the electrode. Therefore, the signal that the o'th electrode sees due to the m 'th MU'S spike train sm (t), is the convolution of the two entities, given by:

---

Under review as a conference paper at ICLR 2020
Teld Σ hem [mmmlmmmll 1) no
where ne[t] is the thermal noise at the e'th electrode. Of course, there are more than one MUs in fact there are M of them, so the total result of what the e'th electrode sees is sum of all their
effects, and we have:
M-14-1 Σ m=0
xe
hemm [msm[ 11 + ne(t]
(2)
METHOD
The spike sorting problem can be cast an unsupervised representation learning problem (Good- fellow et af.] 2016) [Bnniio et al. 2013), where the multi-channel sEMG signal represented as individual MU spike trains and the corresponding MUAP spatiotemporal waveforms. Here, tackle this problem using a deep autoencoder where the encoder extracts spike pulse trains and the decoder reconstructs the sEMG signal according the biophysics of the problem. This network is then trained using loss terms that enforce the reconstruction of the sEMG signal, the sparsity of the spike trains, and the parsimony of the MUAPs, the uniqueness of motor units, and the refractory period of motor neurons.
3.1 MODEL ARCHITECTURE
The DSD model can be put into the following form where Equations 3 and represent the encoder and the decoder parts correspondingly. The encoder fe (-) takes in the SEMG data x E R RET where E is the number of electrodes and T is the number of timesteps, and the encoder weights wa. to infer the spike probabilities å E R NaT where N is the maximum number of MUs the network is allowed learn. In this representation, individual MU spike probabilities are inferred for all MUs for all timesteps.
$ fe (r; we)
c = fa (8: wa)
(4)
The decoder fa (-) then takes in the inferred MU spike probabilities and infer the reconstructed sEMG signal using the decoder weights wg. Since the decoder mimics the biophysical model closely, we can recover the MUAP waveforms corresponding to each MU from the decoder weights wa.
Here, the encoder model corresponds to the learned deep spike sorter and the decoder model corre- sponds the learned biophysical model. During inference time, we typically run only the learned encoder for spike sorting purposes.
DSD Encoder The encoder function fa (.) constitutes a deep convolutional neural network (CNN) which is tasked with solving the spike sorting problem. While traditional spike sorting approaches first identify peaks and then cluster them into motor units, our deep spike sorter solves both problems simultaneously in data-driven manner.
The encoder is implemented as a segmentation network followed by a classifier head. The segmen- tation network takes in the sEMG data R ExT and encodes it in a latent space of size (L, E, T)) where L is the dimensionality of the latent space. The classifier head then takes in this latent repre- sentation and calculates spike scores using linear classifier, with input size (I, E), for each motor unit that runs through each timestep, which are then passed through a Gumbel-soffnax function (Jang al. 2016) for spike binarization.
We have two variants of the DSD model based on the segmentation network used. In one variant, we densely connected neural network architecture (i.e. DenseNets) similar to Huang al. (007).
2

---

Under review as a conference paper at ICLR 2020
The other variant utilizes skip connections between dense blocks and is based on the hundred layer tiramisu architecture (Jégou et al. 2017).
DSD Decoder The decoder function (-) is strictly linear model that mimics the convolutive multi-channel sEMG mixing model. This biophysical model, described in detail by Merletti & Farina (2016), can be deseribed as single convolution operation where the convolutional kernels are the spatiotemporal waveforms corresponding to each motor unit action potential (MUAP). Using this biophysical model helps constrain the problem such that, together with the training losses, coincides with the spike sorting problem. And the interpretability of it allows us to easily recover the MUAP waveforms corresponding each MU.
In our implementation, we used the canonical tensor decomposition trick (Kolda & Bader 2009 Hitchcock 1927), alongside the spatial and the temporal dimensions to significantly reduce the number of learned parameters while maintaining the expressive power necessary to span the space of possible MUAP waveforms.
3.2 TRAINING LOSSES
Reconstruction Loss of the reconstruction error.
The reconstruction loss is given in Equation It penalizes for the I4 norm
Lrecomstraction
|| 111
(5)
Sparsity Loss MU spikes can be rare events relative to typical sampling rates of sEMG devices. We use the L1 norm the spike probabilities to enforce this prior knowledge, as shown Equa- tion 6
Espasity
||$||
(6)
Parsimony Loss One of the main challenges in using a deep autoencoder for spike sorting is the cardinality problem. Since we do not know the exact number of motor units, overparameterize the network based on an upper limit. Then, by penalizing for the proposal of new MUAP waveforms, are able to force the network to be parsimonious with respect to the set of MUAP proposals uses to explain given sEMG signal. We use the Lı-L2 norm to enforce this group sparsity as in Bach et al. 2012
Ow Øt
paasimony
gEg
Equation| 7 describes the parsimony loss. The partitioning for L1-L2 norm is defined such that each partition, G. corresponds to a tensor w that describes the local spatiotemporal waveform around a neighborhood the number of MUAPs and their spatial footprints. Also note that apply this loss not directly the spatiotemporal waveforms, but rather on their first derivatives with respect to time, in order to enforce temporal smoothness.
electrodes for a particular MU. This forces the network to minimize both
Uniqueness Loss The uniqueness loss forces the MUAP waveform proposals to focus on explain- ing different phenomena. It implements a joint prior on the set of MUAP waveforms and their corresponding MU spike trains. This is done by penalizing for the joint occurrence of the temporal similarity between MU spike trains and the spatial similarity between their MUAP waveforms.
Luniqueness ||vee ((352)
? (rmm") (J-1))ll
In Equation the terms and m correspond to dilated version of the spike trains and the spa- tial energy distribution of the MUAP waveforms. The matrices J and I correspond to all-ones matrix and an identity matrix respectively. The operators vec (.) and @ correspond to vectorization operation and the hadamard element-wise product.
3

---

Under review as a conference paper at ICLR 2020
Refractory period of motor umits describe the amount of time needed be- tween two spikes from a single motor unit. The refractory period loss encodes this prior knowledge by penalizing for multiple spike occurences from single MU within its time window. It can be implemented efficiently as Equation
Refractory Period Loss
Leefractory
|ReLU (sumpool (s) 1)ll1
(9)
Total Loss The total loss is combination of the above loss items shown in Equation 10 where a, B3 n, and are hyper-parameterd to be tuned.
Liotal
rcconstruction + sparsiiy 8 parrimony + Luniqueness mererretorr
(10)
RESULTS
4.1 EVALUATION WITH REAL-WORLL DATA
We first evaluate the DSD algorithm real sEMG data acquired from 16 differential sensors worn around the circumference of the arm. The sEMG data recorded is shown in Figure [ In this case, the human user was instructed to generate (imperceptible) tensing of muscles, first of the index finger and then of the small finger. Because our sEMG datasets lack simultaneously recorded intramuscular EMG serve as ground truth signals Negro al. 2016), this evaluation is limited to unsupervised analyses.
separate activity modes, corresponding
slight
Its subsequent blind factorization by the DSD shown Figure We can clearly see the following facets: MUAPs are spatially localized with a clear spatial peak and fall-off. Their smoothness and time-profiles look reasonable, as expected for differential recordings, and the MUAPs look like they roughly factorize the data. Furthermore, their corresponding spike-trains also seem to spike at the appropriate peak locations. Lastly, the cardinality of the triks looks reasonable as well, in this case parsimoniously explaining the dataset using two unique waveforms. The zoomed in versions of both results are shown for the orange and green waveforms in Figure the reconstruction errors more saliently, and the left hand side,
S On the right
hand side, some missed detections corresponding to what look like low amplitude spikes.
4.2 QUANTITATIVE RESULTS WITH SIMULATED DATA
To quantify the DSD performance, we generate a family of simulated datasets using real-data (see Appendix C similar the hybrid ground truth approach used by byPachitarii et (2016) and Rossant et al. (2016). The main goal here is to compare MUAP waveforms similarities and corresponding
EMG Sigral
Figure Recorded 16-channel wrist sEMG data with twitches index finger in the first active section, and twitches in pinkie in the second active section.
4

---

Under review as a conference paper at ICLR 2020
Figure 2: DSD-D factored sEMG data with 10 learned MUAP templates on the right, 10 corre- sponding spike trains on the upper left, and spike firing time overlaid the sEMG the lower left.
Figure 3: DSD-D factored sEMG data, zoomed-in sections of the spike trains and sEMG. On the left shows the pinkic active section, and on the right shows the index finger active section.
5

---

Under review as a conference paper at ICLR 2020
Table I: sEMG decomposition evaluation results on simulated datasets
Dataset
Metric
KiloSort
DSD-T
DSD-D
3 MUs
recall
0.6667
1.0000
1.0000
dissimilarity
0.7112
0.0802
0.0937
cardinality mismatch
-0.3333
0.0000
0.0000
overall_avg-pa
0.6240
0.7899
0.7818
overallavg-prr g-ppt
0.8089
0.7570
3 MUs w/ high noise
recall
0.6667
1.0000
1.0000
dissimilarity
0.6541
0.0802
1.36
cardinality mismmtch
-0.3333
0.0000
0.0000
overall_ave avg--pr
0.3532
0.8063
0.7428
overalLavg-prt
0.8126
0.7434
5 MUs
recall 0.6000
1.0000
1.0000
dissimilarity
0.6068
0.3247
0.2255
cardinality mismatch
0.0000
0.0000
0.4000
overall avg prt overall_avg_prg Lavg-prt*
0.5955
0.4539
0.7553
0,4835
0.6210
7 MUs
recall 0.1429
0.8571
1.0000
dissimilarity
0.7786
0.1484
0.2082
cardinality mismatch
0.0000
0.0000
0.1429
overall_avg-pp
0.6452
0.5849
0.6535
overallaav-ppr
0.5868
0.6537
9 MUs
recall
0.4444
0.8889
0.7778
dissimilarity
0,6012
0.1288
0.1188
cardinality mismmtch
-0.2222
-0.1111
-0.1111
overall_ave avg-prt 0.4625
0.6112
0,5492
overallLav--prt
0.6174
0.5412
indicates results on validation datasets for DSD
spike-train timing accuracies against this ground truth. Those are encoded in the following four metrics below:
cardinality mismatch: The normalized cardinality difference between inferred and ground-truth MUAPs.
2. recall: The percentage of MUAPs successfully associated to ground-truth MUAPs.
3. dissimilarity: The reconstruction accuracy of the associated MUAPs.
4. overall avg- prt: The timing accuracies the inferred spike trains.
See Appendix B for detailed description of these metries and how they' 're calculated.
We evaluate two DSD variants, DSD-T, and DSD-D, against each other and also against the estab- lished spike sorting algorithm KiloSort | Paachitariu et al. 2016). DSD-T uses the tiramisu encoder and is trained with reconstruction, sparsity, parsimony, and refractory losses. Similarly, DSD-D uses simply replaces replaces the refractory loss with the uniqueness loss, and has the DenseNet encoder front end.
The DSD-T models were trained with batch size of 32, using the Adam optimizer (Kingma & Ba 2014) with an initial learning rate of 0.001, and the Gumbel-Softmax exponential temperature decay for every 32000 batches from max temperature, 1.5 min temperature, 0.05. Dropouts were added in the tiramisu dense blocks, zeroing out channels with probability 0.3 during training. The DSD- models were trained with batch size of 16, using the Adam optimizer with an initial learning rate of 0.0006, and the Gumbel-Softmax inverse time constant of 0.0001 decreasing temperature every 1024 batches. Dropouts were added the output of dense layers, zeroing out channels with probability 0.01 during training. Both were run using the Pytorch 1.0 framework (Paszke et al. 2017).
K githuh repo: http::/giittb.cco//oneee//////ih

---

Under review as a conference paper at ICLR 2020
PLAP THe Duslay
Figure 4: 3 MUAPs with high noise dataset. DSD-T results of inferred MUAP templates and spike trains (red) and ground-truth (black) Note: The shifted MU-1 template doesn't 't affect the template score, but causes an offset in the corresponding spike train, which affects the overall_avg.pp _agpprt score
For those experiments, the cardinality range of a dataset's 's underlying MUAPs are known, and sets of hyperparameters for both KS-1 and DSD-X were tuned accordingly. This was done assess: i) How easy and meaningful hyperparameter tunings is for both algorithm families, and ii) to assess how well their performance is even if cardinality ranges are known ahead of time. Table summarizes those results.
to
A visualization of the results corresponding to the second of the Table is shown in Figures 5 Qualitatively, we can see that the DSD's ss recovered MUAPs and their corresponding spike- trains are much more aligned with those of the ground truth. (Note the temporal shift of one of the MUAP's signatures) This does not affect MUAP scoring, as this temporal shift is accounted for in the Thifted-Faarnn-diittae (see Appendix B). This does cause a corresponding misalignment of the spike-train however, and this can manifest as penalization in the overall avg prt score, as won't be until a larger window is considered before which the offset spike is accounted for. The kilosort result appears to have MUAPs that are much more oversegmented overall, and one is missing altogether. The corresponding spike-trains also have have one that is also missing.
In general, we can see that the DSD results seem to be much more performant across the board with respect to KiloSort, (at least on those short 30 second datasets). The DSD appears able to recover good amount of the underlying MUAP waveform shapes faithfully, and this appears to be fairly consistent across all datasets. In regards to the spike-train train timing accuracies, the DSD family appears to also perform very well relative KiloSort, although the number of active MUAPs increase, both the KiloSort and DSD families appear to show a beginning trend towards becoming more or less similar in performance. On the dataset with MUAPs for example, the DSD-D is just above KiloSort in performance on spike-train accuracies, and KiloSort beats out DSD-T.
As DSD-T and DSD-D were trained on different architectures and loss parameters, we can also gain insights into what good combinations are useful for which cases. Of important note, was the ease at which we were able to tune the hyperparameters of DSD this context, owing to the interpretability of the loss terms. For example, avoiding MUAP shape splittings was readily tunable via the Luniqueness loss, (which this loss was specifically designed for). Similarly, tuning for excessive spike-trains was facilitated by increasing the L1. We view this interpretability of the DSD's loss functions as an advantage that comes "for free"", since the DNN architecture incorporates biophysical constraints from the ground up.
All training sets were 30 seconds in length. The training time for KiloSort took seconds, and ran on Intel(R) Core(TM) i7-8750H CPU @ 2.20GHz. KiloSort does not have concept of inference time, as it runs directly on a dataset and can immediately recover spike-trains and the corresponding MUAP waveforms. For the DSD, the training time could very anywhere from 10 minutes to 1 hour, depending on the number of epochs utilized. Training the DSD took place on one Vidia GeForce GTX 1070 GPU. Inference time of the DSD, which also took place the GPU, was ms.
7

---

Under review as a conference paper at ICLR 2020
Figure 5: 3 MUAPs with high noise dataset. KiloSort results of inferred MUAP templates and spike trains (red) and ground-true ttuth (black)
CONCLUSION
We present a novel self-supervised learning algorithm for spike sorting of multichannel sEMG data. This model marries the power of deep end-to-end neural-netwworss with the inductive biases of bio- physical models. The architectural novelty is forcing the deep net to reconstruct sEMG through the bottleneck of the biophysical model that created it in the first place, while the loss novelties are around the explainability and interpretability of them vis-a-vis the biophysics. The architec- ture therefore forces the net learn the appropriate generative latent variables the spike-trains, and their corresponding MUAP waveforms, and the loss terms provide clear physical meaning for tuning the model. The algorithm is elegantly supervised, and requires no training labels whatsoever. The model once trained, can then be applied as is to similar data sets to provide spike- eatkeeexttattion. Relative to an established spike-sorting method, this autoencoder-baase recovery of source waveforms and event times at various degrees of superposition, (at least short duration data sets).
approach shows superior
The DSD framework opens up plethora of exciting avenues for the application of self-supervision to complex signals such as sEMG. More broadly however, showcases an exciting and concrete methodology for incorporating physics-based priors hand in hand with the power of end-to-end deep learning.
REFERENCES
Francis Bach, Rodolphe Jenatton, Julien Mairal, Guillaume Obozinski, et al. Optimization with sparsity-induchn penalties. Foundations and Trends in Machine Learning, 4(1):1-106, 2012.
Yoshua Bengio, Aaron Courville, and Pascal Vincent. Representation learning: A review and new perspectives. JEEE transactions on pattern analysis and machine intelligence, 35(8):1798-188 18288, 2013.
David Carlson and Lawrence Carin. Continuing progress of spike sorting in the era of big data. Current opinion in neurobiology, 55:90-96, 2019.
Dario Farina, Roberto Merletti, and Roger M Enoka. The extraction of neural strategies from the surface emg. Journal of applied physiology, 96(4):1486-1495 2004.
Goodfellow, Yoshua Bengio, and Aaron Courville. Deep leurning. MIT press, 2016.
Matthias H Hennig, Cole Hurwitz, and Martino Sorbaro. Scaling spike detection and sorting for next-generation electrophysiology In In Vitro Neuronal Networks, pp. 171-184. Springer, 2019.
Frank L Hitchcock. The expression of tensor polyadic as sum of products. Journal of Mathematics and Physics, 6(1 6(1-4):164-189, 1927.
Ales Holobar and Damjan Zazula. Multichannel blind source separation using convolution kernel compensation. IEFE Transactions Signal Processing, 55(9):448744496 2007.
8

---

Under review as a conference paper at ICLR 2020
Gao Huang, Zhuang Liu, Laurens Van Der Maaten, and Kilian Weinberger. Densely connected convolutional networks. Proceedings of the IEEE conference on computer vision and pattern recognition, pp.
2017.
Eric Jang, Shixiang Gu, and Ben Poole. Categorical reparameterizztioo with gumbel-softmax. arXiv preprint 1rXiv:1611. 1.01 144, 2016.
Simon Jégou, Michal Drozdzal, David Vazquez, Adriana Romero, and Yoshua Bengio. The one hundred layers tiramisu: Fully convolutional densenets for semantic segmentation. In Proceed- ings of the IEEE Conference Computer Vision and Pattern Recognition Workshops, pp. 11-19, 2017.
Diederik P Kingma and Jimmy Ba. Adam: A method for stochastic optimization. arXiv preprint 1rXiv::412.69666 2014.
Tamara G Kolda and Brett W Bader. Tensor decompositions and applications. SIAM review, 51(3): 455-500, 2009.
Harold W Kuhn. The hungarian method for the assignment problem. Naval research logistics quarterly, 2(1-2):83-97, 1955.
Roberto Merletti and Dario Farina. Suffcce electronyography: physioogg,, engineering, and appli- cations. John Wiley & Sons, 2016.
Francesco Negro, Silvia Muceli, Anna Margherita Castronovo, Ales Holobar, and Dario Farina. Multi-channel intramuscular and surface emg decomposition by convolutive blind source separa- tion. Journal of neural engineering. 13(2):22027,, 2016.
Marius Pachitariu, Nicholas Steinmetz, Shabnam Kadir, Matteo Carandini, and Kenneth Harris. Kilosort: realtime spike-sorting for extracellular electrophysiology with hundreds of channels. BioRxiv, pp. 061481, 2016.
Adam Paszke, Sam Gross, Soumith Chintala, Gregory Chanan, Edward Yang, Zachary DeVito, Zeming Lin, Alban Desmaison, Luca Antiga, and Adam Lerer. Automatic differentiation in PyTorch. In NIPS Autodiff Workshop, 2017.
Hernan Gonzalo Rey, Carlos Pedreira, and Rodrigo Quian Quiroga. Past, present and future of spike sorting techniques. Brain research bulletin, 119:106-118 17, 2015.
Cyrille Rossant, Shabnam N Kadir, Dan FM Goodman, John Schulman, Maximilian LD Hunter, Aman B Saleem, Andres Grosmark, Mariano Belluscio, George H Denfield, Alexander S Ecker, et al. Spike sorting for large, dense electrode arrays. Nature neuroscience, 19(4):634, 2016.

---

Under review as a conference paper at ICLR 2020
A GUMBEL-SOFTMAX BINARIZATION OF SPIKE PULSE TRAINS
The spike probabilities predicted by the DSD encoder are binarized using the Gumbel-softmax cat- egorization trick in order avoid partial spikes. This operation is given in Equation 11
(o
P(a))
(1I)
It is similar to the standard softmax operation with two terms added: the temperature parameter and a draw from a Gumbel distribution g G(O, 1). The Gumbel draw can be computed by g -log(-log(u)) where UO, 1) is draw from a uniform distribution.
During training time, as the temperature term т reduced, the distribution of spike probabilities begins to approximate one-hot vector, forcing the DNN to make spike non-spike decisions. During inference time, we drop the stochastic element altogether and round the resulting probabili- ties.
B EVALUATION FRAMEWORK
The evaluation framework evaluates the correctness of MUAP waveforms and Spike Trains recov- ered by spike decomposition algorithms. It does so by comparing ground truth and recovered MUAP waveforms, and also comparing ground truth and recovered spike trains.
B.1 EVALUATING MUAP WAVEFORM RECOVERY
MUAP Waveform recovery is evaluated by first associating each ground truth MUAP waveform to a recovered MUAP waveform and then computing metrics, described in the following sections.
B.1.1 MUAP WAVEFORM ASSOCIATION
Each ground truth MUAP waveform is associated to the best matching recovered MUAP waveform. Given the Normalized shifted Farina distance (sfd) metric, as in Merletti & Farina (2016). between pairs of ground truth and recovered MUAP waveforms, matches found using the Hungarian algorithm (Kuhn 1955).
sfd is computed as follows:
normalize(mu) mu /
(12)
mu
(mergy(mu) mm
(13)
2 energy(mul mu2 aligned)
sfd(mul, mu2)
(14)
energy(full) +
where is the per channel mean of the corresponding MUAP, mul and mu2 are normalized as above to yield mul and mu2, mu2 is temporally aligned to ruul to yield rhu2_aligned such that their correlation is maximum.
B.1.2 RECALL
Recall computed as the fraction of MUAP waveforms recovered correctly. A value of 1 is the best score and value of 0 is the worst.
Noorveetocooveies
recall
(15)
N muaps
Here, Nonnet recoveries is computed as number of instances where sfd between ground truth and it's associated MUAP waveform below a threshold. (default: 0.9). N muaps is the number motor units in the ground truth.
10

---

Under review as a conference paper at ICLR 2020
B.1.3 CARDINALITY MISMATCH
Cardinality mismatch is computed the difference in number of recovered and ground truth MUAP waveforms normalized by the number of ground truth waveforms.
N muuss
Mgummaps
cardinality_minn wwimmttcchh..ccccmmmmmmmmmmm.
(16)
Ngt mtaps
This measure reflects scenarios where the algorithm failed to recover the exact number of ground truth MUAP waveforms. value of 0 is the best score. A value < indicates that we failed recover some MUAP waveforms, and value indicates that recovered more MUAP waveforms than ground truth.
B.1.4 DISSIMILARITY
For all correct recoveries the dissimilarity computed the average sfd distance between ground truth MUAP waveforms and it's associated MUAP waveforms.
N Namettaommri sfd (gt.mu,, associated.mmu
dissimilarity
(17)
This measure indicates how accurately recovered MUAP waveforms were reconstructed. A value of O is the best score, anything higher is progressively worse
B.2 EVALUATING SPIKE TRAIN RECOVERY
Spike Train recovery evaluated by first associating each ground truth Spike Train to a recovered Spike Train and then computing metrics, as described in the following sections.
B.2.1 SPIKE TRAIN ASSOCIATION
Each ground truth spike train needs associated to the best matching recovered spike train. Ground truth spike trains are binary. The recovered spike-trains are probabilities.
Let us assume that there were no alignment issues between ground truth and recovered spike trains. In such case, one could simply measure the L1 distance between ground truth and recovered spike trains, and come up with the best association.
Realistically however, alignment issues between spike trains can exist, even though they may have the same inherent firing pattern.
Therefore we modify the association measure, to incorporate mis-alignments. We do this in the following way: The LI measure is first applied to all samples across ground truth and recovered spike trains, except for all regions around ground truth spikes.
Those regions, are defined by +/- 20 samples around any given ground truth spike. Within this region, we simply pick the max spike probability from the recovered spike train, and use this value, we ignore all other values within the region.
For every region, now have list of ground truth spikes and list of max probability values corresponding to the ground truth spikes. We simply compute the LI measure of the two lists and sum the Ll measure computed above. This is the associated spikes distance (asd).
Given the asd metric between pairs ground truth and recovered spike trains, matching spike trains found using the Hungarian algorithm (Kuhn 1955
B.2.2 OVERALL AVERAGE PRECISION- RLECISION-RECALL- TOLERANCE
Each ground truth spike train needs to be evaluated against the best matching recovered spike train. Ground truth spike trains binary. The recovered spike-trains are probabilities.

---

Under review as a conference paper at ICLR 2020
Let assume there alignment issues between ground truth and recovered spike trains. In such a case, one could measure algorithm performance by computing a precision-recal curve for ground truth and its associated spike train.
The average precision (ap) score computed from the precision-recall curve above could used to measure recovery performance for single spike train.
However, alignment issues between spike trains can exist. Therefore we modify the evaluation measure, to incorporate mis mis-alignments.
Let us consider regions, defined by +/- samples tolerance window around ground truth spike, we detect the max spike probability from the recovered spike train in these regions. We temporally align the max spike probability with the ground truth spike in that region. All the other values within this region left is.
We now have aligned spike trains, for which we can calculate a single precision-recall curve and its average precision score.
We can repeat this process for a range of tolerance window values, to obtain a bank of precision recall curves and list of average precision scores. This accommodates range mis-alignment errors.
We can average all the average precision scores to obtain the precision-recaul recall-tolerance (prt) score for one spike train. This score measures the accuracy of spike train recovery.
Overall average precision-recall tollacceeeeeeeee is simply the average of all prt scores for all spike trains. This is measure of how accurately we can recover all the spike trains from sEMG data. A value of is the best score and value of the worst.
E Naw
Di=l
prt(gt associated
average_prt
(18)
N Ngs
Σ Dwinso ap(stı, st2, win)
ptt(st),st2)
(19)
20
We pick maximum tolerance window (win) of 10ms around ground truth spike to represent half the minimum interval between consecutive spikes firing from the same motor unit in the forearm. We increment win sizes by 0.5ms.
12

---

Under review as a conference paper at ICLR 2020
C SIMULATED DATASET GENERATION
The simulated datasets are 30 seconds long, 60,000 samples in length, with 16 channels in width. They are generated using the extracted MUAP templates and spike trains from real-world sEMG datasets. It is a human involved extraction process to ensure the quality of the spike trains and MUAP templates. After we have bank of MUAP templates and spike trains, we select several ones we want to use in a dataset. The dataset generation process is convolution of the MUAP templates with the spike trains and then adding Gaussian noise. In this way, the simulated datasets closly resemble the real-world sEMG datasets, but with the ground-truth spike trains and MU templates for quantitative evaluation. We generate these datasets with different number MUs and superposition in order to vary the spike sorting difficulty level and thus test the algorithms thoroughly.
Figure 6: Simulated dataset with ground-truth MUAP templates and spike trains
More specifically, there are 5 simulated datasets we use in evaluation. The 3 MU dataset has 3 MU templates with channel offset which is the most simple and non-trivial case. The MU with high noise dataset has the same MUs as the 3 MU dataset except the noise power is increased by 4 times, which is test the performance
lower SNR datasets. The MU dataset has MU templates on the same channel and 2 MU templates on another channel, which is to test the capability of addressing the spatial superposition problem. The MU dataset has 5 MU templates with one channel offset, and 2 MU templates on the same channel. The 9 MU dataset has 5 MU templates with one channel offset and MU templates on the same channel. The last two datasets are the stress-test cases by increasing the amount of MUs and spatial superpostion. Figure shows the example of MU simulated dataset on the right, the MUAP templates used on the upper left, and the spike trains used on the lower left.
13


Reviews:
[{"review_id": "rJxSU_FbqB", "reviewer": "AnonReviewer1", "summary": "", "questions": "", "limitations": "", "rating": 1, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "This paper describes a machine learning model for learning spiking representations from EMG (electromyography) like data. The basic model is an autoencoder, where the bottleneck layer is designed to project the waveform inputs into \"spike\" representations. Specifically, there is a DSD (Deep Spike Decoder) encoder and decoder, where the encoder is implemented as either a DenseNet[1]-like or Tiramisu[2]-like architecture. The bottleneck seems to be implemented by a linear layer binarized with a Gumbel-Softmax activation. The decoder is a linear layer. Several losses are considered, including a L4-norm (!), a sparsity loss, a parsimony loss, a uniqueness loss, and a refractory period loss. The model is validated qualitatively on real EMGs, and quantitatively on synthetic data.\n\nThe strengths of the paper are the following:\n- Spiking models are very interesting class of models, which if attained would have a great impact on several other areas of machine learning. \n\n- I like the straightforward design that is fit to the purpose of generating spiking representations. A Gumbel-Softmax has proven its validity and is a logical fit to the problem setup.\n\n- I like the simpler setup the paper chooses.\n  -- First, the paper chooses not to take the route of trying to learn backpropagation through infinite time steps, as often happens in spiking methods. This is a beast on each own, which would be nice to set as an ultimate goal (in my opinion). For now, however, it's ok if we forego this requirement.\n  -- Second, the paper chooses not to assume that signals come asynchronously, which again makes things unnecessarily complex, given the state of the field.\n\nThe weaknesses of the paper:\n- While the problem setup is simple, perhaps it makes some assumptions that are too strong. In my opinion, the strongest one is that of batch-learning (for the lack of a better name). Often, spiking models are studied in an online setup where the data comes online (continuously), and then spikes are generated when sufficient changes in the inputs warrant a delta-difference (spike). When in batch mode, however, it should be quite much easier to obtain spikes that lead to good reconstructions, as there is no much need for a \"memory\" mechanism. While one could argue whether an online setup is necessary or not, in my opinion it is necessary and would make a spiking model challenging and interesting to learn. Otherwise, it looks like a conventional autoencoder, only with spikes instead of dense representations.\n\n- The model is unclear and writing vague.\n\n- First of all, after reading the abstract and the introduction I get the feeling that the model is probabilistic, as there is mention of priors and autoencoders. Also, later on a Gumbel-Softax is mentioned for binarization. Gumbel-Softmax is a continuous function and binarization makes sense when sampling, that is when assuming a generative process. However, the rest of the paper seems not to explain a probabilistic or generative model. There is no explicit prior discussed. There is no sampling procedure discussed. The losses are explained but not within a particular probabilistic or stochastic framework. If the model is not stochastic, one way or the other, how are the discrete spikes obtained and how is the system is trained?\n\n- All the loss functions, albeit logical when taken one by one, they do look ad hoc and appearing out of the blue. This perhaps relates to the previous point, where it is not clear if the model is stochastic or deterministic. If it is deterministic, it is ok to have all the loss functions appearing like that, but I would expect some more explanation on what purpose do they fill w.r.t. the spiking model. In the end, what does the deep spiker model try to achieve? Learn spikes as representations? Recover the original spikes? Be sparse? If yes, why and how sparse? Be energy-efficient?\n\n- Some design choices are quite unclear. Generally, it is fair to say that the experimental and design setups are rather simple: multiple 1-D waveforms and not much noise (from what I get). In that context, it is not clear why DenseNets, or even Tiramisu-Nets are used as an encoder; especially when the decoder is a simple linear model. \n\n- Also, phrases like \"what good combinations are useful ... the hyperparametres of DSD\" do not add to the clarity. There is no exploration of hyperparameters in the experiments and no individual examination of the contribution of each loss (unless I missed it somewhere).\n\n- Similarly, what does \"For the DSD ... between 10 minutes to an hour depending on ...\".  Such statements should be more precise, for instance plotting wall clock time vs training loss.\n\nAll in all, while I like the motivation and the original direction, I believe there exist a lot of questions  unanswered before acceptance.", "extra_scores": {"experience_assessment": "I have published in this field for several years.", "review_assessment:_thoroughness_in_paper_reading": "I read the paper at least twice and used my best judgement in assessing the paper.", "review_assessment:_checking_correctness_of_experiments": "I carefully checked the experiments.", "review_assessment:_checking_correctness_of_derivations_and_theory": "I carefully checked the derivations and theory."}, "extra_text": {}}, {"review_id": "BJeR3p3xir", "reviewer": "AnonReviewer4", "summary": "", "questions": "", "limitations": "", "rating": 1, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "This paper proposes a new algorithm for spike-sorting. It is implemented by a deep autoencoder with biophysically motivated loss functions. The encoder is the main module which conducts the spike-sorting task while the decoder is only used for training the model in an end-to-end fashion. \n\nThis paper should be rejected due to the following arguments:\t\n- The paper lacks a section on literature survey, to let the reader know how/where the proposed method fills the gap in the current state-of-the-art. They do compare their results with the KiloSort (Pachitariu et al., 2016) algorithms, however, no discussion is provided on how it works and why their method outperforms it. \n- It is unclear why the reconstruction loss is chosen to be an L4 norm as opposed to L2.\n- The authors claim that the parsimony loss as defined in Eq. (7) forces “the network to be parsimonious with respect to the set of MUAP proposals it uses to explain a given sEMG signal.” My understanding, however, is that the only functionality of the loss defined in Eq. (7) is to enforce temporal smoothness. More elaborate explanation is needed to support the authors claim.\n- I could not understand the functionality of the uniqueness loss. Specifically, why should “the joint occurrence of the temporal similarity between MU spike trains and the spatial similarity between their MUAP waveforms” be penalized? Isn’t that the case that same stimuli should result in similar response? It is unclear what this has to do with forcing to explain different phenomena.\n\nThings to improve the paper that did not impact the score:\n- The method (and the paper) is named deep spike “decoder” (DSD) while in fact the “encoder” part of the learned deep autoencoder actually conducts the spike-sorting task. This could be confusing!\n- Page 2, Sec. 3.1, line 2: Should use \\times in inline equations in Latex for the multiplication symbol, not character x. Fix everywhere in the text.\n- Page 6, Par. -2, line -2: The word “replicate” is repeated.\n- Non-legible plots axes.\n\n", "extra_scores": {"experience_assessment": "I do not know much about this area.", "review_assessment:_checking_correctness_of_experiments": "I assessed the sensibility of the experiments.", "review_assessment:_thoroughness_in_paper_reading": "I read the paper at least twice and used my best judgement in assessing the paper.", "review_assessment:_checking_correctness_of_derivations_and_theory": "I carefully checked the derivations and theory."}, "extra_text": {}}]

Author rebuttal:

Meta-review: The paper presents a model for learning spiking representations. The basic model is a a deep autoencoder trained end-to-end with a biophysical generative model and results are presented on EMG and sEMG data, with the aim to motivate further research in self-supervised learning.

The reviewers raised several points about the paper. Reviewer 1 raised concerns about lack of context on surrounding work, clarity of the model itself and motivating the loss. Reviewer 2 pointed out strengths of the paper in its simplicity and the importance of this problem, but also raised concerns about the papers clarity, again motivations on the loss function and sensibility of design choices. The authors responded to the feedback from reviewer 1, but overall the reviewer did not think their scores should be changed.

The paper in its current form is not yet ready for acceptance, and we hope there has been useful feedback from the reviewing process for their future research.
Decision: Reject