Bypassing spike sorting: Density-based decoding using spike localization from dense multielectrode probes
Authors: Yizi Zhang
Tianxiao He
Julien Boussard
Charlie Windolf
Olivier Winter
Eric M. Trautmann
Noam Roth
Hailey Barrel
Mark M Churchland
Nick Steinmetz
Erdem Varol
Cole Lincoln Hurwitz
Liam Paninski
Venue: NeurIPS 2023 spotlight (2023)
Keywords: neural decoding
brain-computer interfaces
spike sorting
variational inference
generative models

Abstract: Neural decoding and its applications to brain computer interfaces (BCI) are essential for understanding the association between neural activity and behavior. A prerequisite for many decoding approaches is spike sorting, the assignment of action potentials (spikes) to individual neurons. Current spike sorting algorithms, however, can be inaccurate and do not properly model uncertainty of spike assignments, therefore discarding information that could potentially improve decoding performance. Recent advances in high-density probes (e.g., Neuropixels) and computational methods now allow for extracting a rich set of spike features from unsorted data; these features can in turn be used to directly decode behavioral correlates. To this end, we propose a spike sorting-free decoding method that directly models the distribution of extracted spike features using a mixture of Gaussians (MoG) encoding the uncertainty of spike assignments, without aiming to solve the spike clustering problem explicitly. We allow the mixing proportion of the MoG to change over time in response to the behavior and develop variational inference methods to fit the resulting model and to perform decoding. We benchmark our method with an extensive suite of recordings from different animals and probe geometries, demonstrating that our proposed decoder can consistently outperform current methods based on thresholding (i.e. multi-unit activity) and spike sorting. Open source code is available at https://github.com/yzhang511/density_decoding.

Bypassing spike sorting: Density-based decoding using spike localization from dense multielectrode probes
Yizi Zhang
Tianxiao Hel.2+ Julien Boussard' Charlie Windolf'
Olivier Winter
Eric Trautmann
Noam Roth' Hailey Barrell
Mark Churchland
Nicholas A. Steinmetz
The International Brain Laboratory's
Erdem Varol
Cole Hurwitz
Liam Paninski
Columbia University
New York University
Thee International Brain Laboratory
University of Washington in Seattle
Abstract
Neural decoding and its applications to brain computer interfaces (BCI) are es- sential for understanding the association between neural activity and behavior. A prerequisite for many decoding approaches is spike sorting, the assignment action potentials (spikes) individual neurons. Current spike sorting algorithms, however, can be inaccurate and do not properly model uncertainty of spike assign- ments, therefore discarding information that could potentially improve decoding performance. Recent advances in high-density probes (e.g., Neuropixels) and computational methods now allow for extracting rich of spike features from unsorted data; these features can in turn be used to directly decode behavioral cor- relates. To this end, we propose a spike sorting-free frren decoding method that directly models the distribution of extracted spike features using mixture of Gaussians (MoG) encoding the uncertainty of spike assignments, without aiming to solve the spike clustering problem explicitly. We allow the mixing proportion of the MoG change over time response to the behavior and develop variational inference methods to fit the resulting model and to perform decoding. We benchmark our method with an extensive suite recordings from different animals and probe geometries, demonstrating that proposed decoder can consistently outperform current methods based on thresholding (i.e. multi-unit activity) and spike sorting. Open source code is available at https://gihhbb.co///hanng5// l/densstt decoding
1 Introduction
Decoding methods for large-scale neural recordings are opening up new ways to understand the neural mechanisms underlying cognition and behavior in diverse species (Urai al. 2022). The emergence of high-density multi-electrode array (HD-MEA) devices introduced a tremendous increase in the number of extracellular channels that can be recorded simultaneously (Jun al. 2017 Steinmetz et al. (021),, leading to scalable and high-bandwith brain computer interfaces (BCI) systems (Musk al. 2019 Paulk et 2022).
Traditional neural decoding methods assume that spiking activity has already been correctly spike-sorted. As result, these methods are not appropriate for situations where sorting cannot be performed with high precision. Despite intensive efforts towards automation, current spike sorting algorithms still require manual supervision to ensure sorting quality (Steinmetz et al. 2018). Even after careful curation, current spike sorters suffer from many sources of errors including erroneous spike assignment (Deng et al. 2015). The dense spatial resolution of HD probes makes some known issues of spike sorting more evident. With the increased density of the recording channels, the
37th Conference on Neural Information Processing Systems (NcurIPS 2023).

---

probability of visibly overlapping spikes (spike collisions) is higher (Buccino et al. 2022 Even for the same HD dataset, different spike sorters have low agreement on the isolated units and can find a significant number of poorly sorted and noisy units (Buccino et al. 2020). Consequently, only sin- gle units that pass quality control metrics are included in many neural coding studies (IBL. et 2022).
Because the soike-sorting problem remains unresolved, alternative approaches that do not rely on sorted single-units for decoding have been proposed. A popular choice is multi-unit threshold crossing that uses spiking activity on each electrode for decoding (Fraser et al., 2009 Trautmann et al. 2019). However, this approach ignores the fact that the signal on each electrode is combination of signals from different neurons, thus making inefficient use of the data (Todorova et al. 2014).
Ventura (2088) proposed a spike-sorting free decoding paradigm that estimates neuronal tuning curves from electrode tuning curves and then infers the behavior of interest using the estimated tuning curves and newly observed electrode spike trains. More recently, [Chen et al. (2012); K Kloosterman
jtt al. (2014): : Deng et al. (2015); Rezaei et al. (202T) developed spike feature decoding methods that use marked point processes to characterize the relationship between the behavior variable and features of unsorted spike waveforms. However, these decoders based on state-space models make explicit assumptions about the underlying system dynamics which reduce their flexibility in capturing complex relationships in the data. Moreover, these methods mainly utilize simple waveform features for decoding such as the maximum amplitude on each electrode and do not take advantage of HD spike features such as the estimated spike location.
To leverage the spatial spread and density of HD probes, Hurwitz et al. (2019); Boussard et al. (2021) developed spike localization methods. These methods estimate the source location of detected spike; this is a low-dimensional feature that is informative about the firing neuron's identity. We propose probabilistic model-based decoding method that scales to HD-MEA devices and utilizes these novel localization features in conjunction with additional waveform features. We use a mixture of Gaussians (MoG) model to encode the uncertainty associated with spike assignments in the form of parametric distributions of the spike features. Unlike traditional MoG models with a fixed mixing proportion, method allows the mixing proportion to depend on the behavior of interest and change over time. This is motivated by the theory that behavioral covariates that modulate neurons' firing rates also contain information about spike identities and that such tuning information should be incorporated into spike sorting and neural decoding in order to obtain unbiased and consistent tuning function estimates (Ventura 2009). To infer the functional relationship between spike features and behavioral correlates, we employ automatic differentiation variational inference (ADVI) (Kucukelbir al. 2017) and coordinate ascent variational inference (CAVI) (Blei al. 2017), which enable us to perform efficient and accurate inference while considering the behavior-modulated MoG model.
We apply our method to a large number of HD recordings and decode various types of behavioral correlates. Experimental results show that our decoder consistently outperforms decoders based on multi-unit threshold crossings and single-units sorted by Kilosort 2.5 (Pachitariu et al. 2023). We further validate the robustness of our method by applying it to recordings with different levels of sorting quality, HD probes with varying geometry, and recordings from multiple animal species. Con- sistent with previous results, our findings indicate that relying solely on "good" units, as determined by sorting quality metrics, leads to information loss and suboptimal decoding performance. This observation motivates our transition spike sorting-free decoding framework which enables to extract more information from the spiking activity and improve decoding performance.
2 Method
Consider an clectrophysiological recording comprised of K trials, where each trial is divided into T equally spaced time bins. Let {suR li Sat E R denote set of spike features, where i indexes the i-th spike, ny represents the number of spikes collected in the t-th time bin of the k-th trial, and D is the dimension of spike features. For example, the spike feature Sak = (Zuk, Ritk, axek) E R includes the spike location along the x- and z-axis of the probe, and its maximum peak-to-peak (max ptp) amplitude. Let yı: E RT be the observed time-varying behavior in the trial k; e.g., the speed of rotating wheel controlled by a mouse. When the behavior in the trial i does not vary over time, can take on either binary (y 13, or scalar (y E R) value, e.g., the mouse responds to
2

---

scalar- scalar-valued valued stimulus by making a binary decision.
The proposed decoding method comprises an encoder and a decoder model. During the training of the model, the encoder captures the relationship between the observed spike feature distribution and the observed behavior. During the testing phase, the decoder utilizes the learned relationship from the encoder predict the unobserved behavior based on newly observed spike features. In the following section, we present general formulation of the encoding-decoding paradigm and provide more detailed implementations the supplementary materials.
a), NP1
50c
1000
1500
(um)
2500
5502
3500
s0n
600
(ur))
12:00
14h0
1600
b)
training lmp
lime
splike feantres
decoded behaaiir
Isal trials
decader
Oc
0g
noural recordings
estimated behavior from threshokding
Art
Tate
encoddr
Bith
Pi
Pr
Ms
Falk
En
Ex
fring raleo
X
weightt matri
Figure 1: Decoding paradigm and graphical model. (a) Spike localization features, (ar, 2), the locations of spikes along the width and depth of the NP1 probe, and waveform features, a, the maximum peak-to-peak (max ptp) amplitudes of spikes. Amplitude is measured in standard units (s.u.). Spike features from the entire probe are shown, and we focus on a specific segment of the probe. (b) During the training phase, the encoder takes the observed spike features and behavior y from the train trials as inputs and then outputs the variational parameters 0 which control the dependence of the firing rate À on the behavior y. At test time, the decoder utilizes the learned model parameters obtained from the encoder and the observed spike features from the test trials to predict the corresponding behavior in the test trials. To ensure reliable decoding of behaviors, initially calculate the during training using the learned @ and observed behaviors y from the train trials. Then, compute the À during test time using the learned and the estimated behavior obtained through multi unit thresholding from the test trials. Finally, we generate the weight matrix W for both the train and test trials input to the final behavior decoder, e.g., linear regression or neural networks (Glaser et al. 2020; Livezey and Glaser (021). (c) In the encoder, the firing rates of each MoG component Aath modulated by the observed behavior yeu in the train trials. This modulation affects the MoG mixing proportion Teck. which in turn determines the spike assignment Zitk that generates the observed spike features sitk the train trials. the decoder, the behavior Yek in the test trials is unknown and considered as latent variable to be inferred. The decoder uses the observed spike features Sitk from the test trials along with the fixed model parameters O, learned by the encoder to infer the latent behavior Ytk.
3

---

Encoder. The multivariate spike feature distribution is modeled using a mixture of Gaussian (MoG).
The encoder generative model is as follows:
Aeth = Mt,yk. 0c), Pg p((0),
(1)
Natk
Pitk ~ Categorical Cask; Tek), Tek = {Tetk 1 Tetk = Ee Aerek"
(2)
Sitk ~ (N8ti Tess)) T) Ee),
(3)
where X(-) is function that describes the firing rate's dependence on behaviors yk. while p(0x) represents a general prior on 9. encompassing the model parameters for the mixture component c. Intuitively, the behavior-dependenn À governs the mixing proportion π of the MoG, which determines the specific component e from which a spike feature 8 is generated. As varies over time in response y. spikes that spatially close and share similar waveform features may originate from different MoG components at different time points within a trial. In our implementation, we parameterize using a generalized linear model (GLM), but alternative models such as neural networks (NN) can also be used; see the supplementary material for the GLM configuration.
We employ variational inference (VD) to learn the unknown quantities. In the standard MoG setting (Blei et al. (2017), our goal is to infer the spike assignment which indicates the latent component from which the observation s originates. However, unlike the standard MoG, our spike assignment influenced by the firing rates Ả of the neurons which are modulated by the behavior y- Consequently, learning the association between À and necessitates the estimation of the unknown model parameters 0. Our objective is to simultaneously learn both the latent variables 2 and model parameters based on the observed spike features and behavior To accomplish this, we posit a mean-field Gaussian variational approximation
((2, 0) = II (E<)))((O.) c,t
for the posterior p(2. 0 s. Subsequently, we employ the CAVI or ADVI methods to maximize the evidence lower bound (ELBO) and compute updates for z and 9. Analogous to the standard Detectation-Maximizuee
(EM) algorithm for MoG, the proposed CAVI and ADVI procedures consist of an E step and a M step. The step, which updates 2, closely resembles that of the ordinary MoG, while the M step, responsible for updating 0, differs. CAVI utilizes coordinate ascent to find 0 that maximizes the ELBO, while ADVI employs stochastic gradient ascent for updates. For detailed information on the CAVI and ADVI model specifications, refer to Supplement I and 2.
Decoder. The decoder adopts the same generative model in Equation | as the encoder with two distinctions: 1) ys is unobserved and considered latent variable we aim to estimate, i.e., yk p(ux.) where p(yk) general prior, and 2) the model parameters 0 are obtained from the encoder and kept constant. In practice, the choice of prior relies on the nature of yk- For instance, if yı is binary, we can sample from a Bernoulli distribution while a Gaussian process prior can capture the temporal correlation between time steps if yı: E RR.." The posterior p(z, y 8) is approximated using a mean-field field Gaussian variational approach
q(2, v) II ((<)))(v). c,t
(5)
We employ standard CAVI or ADVI methods to infer and decode y.
Robust behavior decoding. In practice, we found that direct decoding of y using the approximated posterior q(y) in Equation was not robust, leading to decoding results of inconsistent quality across different datasets. Although the factorization described in Equation 5 may not fully exploit the available information for predicting y. it is useful for learning about the spike assignment 2. To enhance decoding robustness, we compute a weight matrix from the MoG outputs as input to the final behavior decoder. The weight matrix, denoted as W, has dimensions K x C' x T and entries Wknt : Σee (( which capture the posterior probability of assigning spike collected at time I of the trial k into the component c. In scenarios such as spike sorting multi-unit thresholding, spikes are assigned deterministically to one of C' sorted single units or thresholded channels. this case, W has one-hot rows, and each entry Wkot represents the number of spikes belonging to trial

---

k, time t and unit (channel) To obtain ((asket), we rely on estimating the posterior Tetk, which requires either the observed Ytk or the estimated fu.- During model training, we can substitute the observed Yek into Equation to calculate the posterior watk for the train trials. At test time, we use the estimated ftk obtained from multi-unit thresholding along with the learned encoder parameters to calculate the posterior Teck for the test trials. With the posterior Tetk in hand, we then estimate q(2) to compute the weight matrix W for both the train and test trials, which serves as input to the final behavior decoder. The choice of the behavior decoder depends on the user's preference. For instance, we can use logistic regression as the behavior decoder for binary ye. and ridge regression for y E R that exhibit temporal variations. Additional information regarding the selection of the behavior decoder can be found in Supplement 5.
3 Experiments
We conducted experiments using both electrophysiological and behavior data obtained from the International Brain Laboratory (IBL) eI al. (02T). The electrophysiological recordings were acquired using Neuropixels (NP) probes that were implanted in mice performing decision-making task. Each recording comprises multiple trials with several behavioral variables recorded during each trial including the choice, face motion energy, and wheel speed; see Figure 2 (a) for details. Each trial has a duration of 1.5 seconds and is divided into 30 time bins of 50 milliseconds length. The NP probe spans multiple brain regions; an example of the brain parcellation can be seen in Figure 3 To prepare the recordings for decoding, we first applied IBL's standard destriping procedure (Chapuis) spike detection and denoising method described in Boussard et al. (2023). After preprocessing, we computed the spike locations (Boussard et al. 2021) to acquire spike features for decoding and then utilized registration techniques (Windolf et al. 2022) to correct for motion drift in the recorded data. Further details about data preprocessing can be found in Supplement In all experiments, we used spike locations along the width and depth of the NP probe, and maximum to-k-t amplitudes of spikes for decoding. We selected this set of spike features based on empirical evidence from our experiments, which showed their good decoding performance. Furthermore, previous studies (Boussard et al. 2023) Hilgen et al., 2017) have also recognized that these features were highly informative about unit identity. Figure a illustrates the spike localization and waveform features that were utilized for decoding.
et al. 2022) reduce artifacts. Then, we used subtraction-basee
We evaluate the performance of our decoding method by comparing it to the following baselines: (1) thiee-thresholded decoders which utilize the spiking activity on each electrode after voltage-based
detection step. (2) Spike-sorted decoders that utilize all single-units found using Kilosort (KS) 2.5 (Pachitariu et al. 2016). (3) Spike-sorted decoders based on "good" units which consist of KS units that have passed IBL's s quality control procedure (IBL et al. 2022). The parameters used for KS were tuned across multiple IBL. datasets as described by IBL's spike sorting white paper (Chapuis et al.. 2022
To assess the quality of decoding, we perform 5-fold cross validation (CV) and compute relevant decoding metrics. The coefficient of determination (2²²) is used to evaluate continuous behavior decoding (e.g., motion energy and wheel speed) while accuracy is used for discrete behaviors (e.g., choice). To demonstrate the efficacy of our approach in a wide range of settings, we conduct the following experiments.
Varying levels of spike sorting quality. Our objective is to compare the proposed decoding method to spike-sorted decoders using datasets which have varying levels of spike sorting quality. We apply our method to two datasets with high sorting quality ("good" sorting) and two with low sorting quality ("bad" sorting). The quality of sorting is assessed using IBL's quality metrics (Chapuis et al., 2022). Although motion registration has been performed, we find that the recordings that have "bad" sortings are more affected by motion drift then the recordings which have "good" sortings; see supplementary materials.
Different brain regions from 20 datasets. To demonstrate the efficacy of our method across many different datasets, we decode 20 IBI. datasets (IBI. et al. 2022). In these datasets, mice perform a behavioral task while NP1 probes record activity from multiple brain regions. These brain regions
5

---

a)
stimulus
molion eraraf
onergy Amm
Mm
mmm
choice (left / right)
wheel speed
good vs. bad sorting example sessions "Dad"
15
1mo (600)
d)
15
sessions
15
R2
good"
"good"
"taad
accuracy "good"
"bad"
R2
eeraee x M8
0.6
0.6
0.8-
0.4
0.5
0.4
motion
wheel speed
0.7-
Fren boed
0.3
motoo energy
abserved
KS
018
2.5 lopicon tammittioo (%6)
0.0
0.7
1.5
1.5
1.5
1.5
motion energy
iime (sec)
0.5
10
200
9401 R2 2010
0.3
15
155
1.5 0 1ime (sec)
1.5
1.5
3.5
5.0
wheel speed
missed sp Kos (16
2.5
0.7
0.0
1.5
0
1.5
1.5
0.5
wheel speed
2-
12 No 0.00
X
1.5
1iie (eec
1.5
155
cumuuativee wboollee (umbbeec)
25
Figure 2: Density-based decoding is robust to varying levels of spike sorting quality. (a) We decode various behaviors including choice, motion energy and wheel speed. In experimental setup, the mouse detects the presence of visual stimulus to their left or right and indicates the perceived location (choice) by turning steering wheel in the corresponding direction. Motion energy is calculated within square region centered around the mouse's whiskers. The example behavior traces are distinguished by different colors for each trial. (b) We compare decoders using two experimental sessions with "good" sorting quality (represented by the color green) and two sessions with "bad" sorting quality (represented by the color purple) based on IBL's quality metrics. The box plots display various statistical measures including the minimum, maximum, first and third quartiles, median (indicated by gray dashed line), mean (indicated by red solid line), and outliers (represented by dots). These decoding metrics obtained from 5-fold CV and averaged across both "good" and "bad" sorting example sessions. (c) We compare the traces decoded by spike-sorted decoders and method example sessions with "good" sorting quality (indicated by green) and "bad" sorting quality (indicated by purple). (d) The scatter plots depict the decoding quality motion energy, measured by Rt, with respect various spike-sorting quality metrics. Each point represents one of the 20 IBL sessions, and different colors and shapes are used to distinguish between the type of decoder and sorting quality. The sorting quality metrics include "contamination." which estimates the fraction of unit contamination (Hill et al. 2011), "drift." which measures the absolute value of the cumulative position change in micrometers per second (um/sec) of a given KS unit, and "missed spikes," which approximates the fraction of missing spikes from a given KS unit (Hill et al. 2011). These metrics are averaged across KS units in session. The scatter plots demonstrate that decoding quality tends to decrease when sorting quality is compromised. However, our method outperforms spike-soree ported decoders even in the presence of these sorting issues.

---

are repeatedly targeted across all the datasets. To explore how behaviors are linked to specific brain regions, we use spikes that are confined to a particular area of the mouse brain for decoding. We collect decoding results from the posterior thalamic nucleus (PO), the lateral posterior nucleus (LP), the dentate gyrus (DG), the cornu ammonis (CA1) and the anterior visual area of the visual cortex (VISa).
Different probe geometry. Our method is capable of decoding electrophysiological data from a variety of HD probes. To demonstrate this, we apply our method on Neuropixels 2.4 (NP2.4) and Neuropixels 1.0 in nonhuman primates (NPI-NHP) datasets. The NP2.4 and NP1-NHP recordings are preprocessed using an identical pipeline as employed for NP1; see Supplement 4 for details. For spike sorting, the NP2.4 and NP1 adopt the same KS parameters as outlined in IBL's spike sorting procedure. Different KS parameters are utilized for NP1-NHP probes which are detailed in Trautmann et al. (002))
NP2.4: Each NP2.4 probe consists of four shanks and total of 384 channels (Steinmetz) et at. (02T). NP2.4 probes are more dense (i.e., more channels in a given area) than NPT. The mice were trained in accordance with the IBL experiment protocols to perform visual decision-making task. The behavioral correlates we decode are choice, motion energy, and wheel speed.
NPI-NHP: NP1-NHP is designed for nonhuman primate species, such as macaques. The NPI-NHP probe maintains the same number of channels NP1 (384), but its overall length is extended, resulting in a sparser configuration compared to NP1 (Trautmann et al. (023). During the experiment, the macaque underwent training in a sequential multi-tareet reaching task (Marshall et al.) (022). The behavioral correlate we decode is the monkey's 's arm force. The probe was implanted the macaque's 's motor cortex.
Comparison to a state-of-the-art clusterless decoder. Although the lack of available code for prior methods make comprehensive comparisons difficult, we benchmark our density- decoder against a state-of-the-a clusterless decoder on datasets from both HD probes and multiple tetrodes. We compare our method the clusterless point process decoder Denovellss et al. (202T), which utilizes marked point process to connect spike features with behaviors. For more details of this comparison, see Section 9 of the supplementary materials.
To decode the binary choice variable, we utilize the CAVI algorithm described in Supplement 2. For continuous behaviors like motion energy, wheel speed, and arm force, we employ the ADVI algorithm outlined in Supplement 1. We specify the maximum number of iterations as hyperparameter in the CAVI model as it requires analytical updates for the model parameters. Running the CAVI encoder and decoder for fewer than 50 iterations yields satisfactory decoding outcomes. As for the ADVI algorithm, we implement it in PyTorc and update the model parameters using the Adam optimizer with learning rate of 0.001 and batch size of 6. The ADVI model is rum for 1000 iterations. Open I/density _dccoding
source code is available at https:/github..om//yaaaa
4 Results
Varying levels of spike sorting quality. The performance of our method in comparison to the spike-thresholded and spike-sorted decoders for both the "good" and "bad" sorting examples is summarized in Figure 2] For the "good" sorting examples, our method has the highest decoding performance for motion energy and wheel speed. For choice decoding, our approach is comparable to decoding based on all KS single-units and better than decoders based on multi-unit thresholding and "good" KS units. For the "bad" sorting sessions, the gap in decoding performance between our method and other decoders is more pronounced. Example traces are illustrated in Figure (c), which demonstrate that the behavior traces decoded by our method closely match the observed traces compared to decoded traces from the sorted decoders. In Figure (d), we quantify the relationship between sorting quality and decoding quality using data from 20 TBL sessions. For all three quality metrics, the performance of our decoder and the spike-sorted decoder decreases as the quality of the sorting decreases. Despite this decrease in performance, performance than the spike-sorted decoder even in the presence of significant motion drift as well as when there is large fraction of missed spikes contaminated units.
method consistently has better
7

---

a)
PO
LP
DG
CA1
VISa
R2
motion energy
0.6
PO
LP
DG
CA1
VISa
0.2
All
Tel
onfe Uhiesin
"Joef 'll K3
duan lhbwad "nnod" al
Cura (fresd
"poo wl
bune driso "good" all KS
durn
Frewah "aooa" cara
KB
KS
K3
K3
KS
KS
KS
KB
wheel speed
+ PO
DG
CA1
VISa
0.4
All
mircsh 'oood" all ours thresh "ppoot" all purs thriss "good"" KS
all KS
ours thresh "pood" all curs throsh "pood" ai ours thrrssh "cood" al ours KS
KS
Ks
NS
Ki
KS
KS KS
KS KS
accuracy
choice
0.9
DG
CA1
VISa
0.8
0.6
All
b)
thresh
"pood" KB
bII KS
thriss "pood"
ours thross "good"
curs thresh 'good"
thresh "oood" KS
all KS
ours
trrssh "pooo" curs
KS K3
KS
KS
KS
KS KS
R2
good KS
ours
0.85
DG
CA1
VISa
0.3
500
200 components
200
100
100
200
R2
components
components
components
components
components VISa
0.85
LP
DG
CA1
0.3
All
1e6
1e6
0.0
1e8
0.2 1e8
0.00
0.25
1e6
spikes
# spikes
spikeste6
# spikes
spikes
# spikes
Figure 3: Decoding comparisons broken down by brain regions. (a) We decoded 20 IBL datasets acquired using NPI probes which inserted into mice performing behavioral task. The locations of the probe insertions the mouse brain and the corresponding brain parcellations along the NP1 probe are shown. We compared the performance of all decoders across different recorded brain regions. For the "All" region, spikes from all brain regions were utilized for decoding. In contrast, for the "PO," "LP," "LP. "DG," "CA1," and "VISa" regions, only spikes from the respective regions were used for decoding. The decoding performance were summarized using box plots showing metrics obtained through a 5-fold CV and averaged across 20 IBL sessions. We observe higher accuracy from PO, LP, and VISa regions when decoding choice; decoding results are more comparable across regions for the continuous behavioral variables. Our proposed decoder consistently achieves higher accuracy in decoding the continuous variables. (b) We use scatter plots to quantify the relationship between decoding quality, measured by R² from decoding motion energy, and the number of components used for decoding. In the case of "all KS" and "good KS", the number of components corresponds to the number KS units. For our method, the number of components refers to the number of MoG components used. For all methods, the decoding performance is higher when using more components (in the regime of a small number components). Our decoding method consistently outperforms spike-sorted decoders based on KS 2.5 while tending to need fewer components.

---

a)
R2
motion
R2
wheel speed
accuracy
choice
0.80
0.8
energy
0.4
0.75
NP2A
0.4
thresh good
all
ours
0.2
thresh good KS
àll ours
0.70
thresh good KS
all
ours
KS
KS
KS
KS
b)
NP1-UUPP
arm force
arm force
observed
KS
0.6
macaria
bran
0.4
5
45
0
0.2
0.658
KS ours
9.85 0 time (s) 9.85
9.85
Figure 4: Decoding performance generalizes across different animals and probe geometry. (a) We compare all decoders on a NP2.4 dataset using box plots showing performance metrics obtained from 5-fold CV. Our method achieves much higher performance than all other decoders on continuous behavior decoding with slightly worse choice decoding than the spike-sorted decoder. (b) We utilize data from single NP1-NHP recording session decode the reaching force of monkey engaged in path-tracking (pacman) behavioral task. The decoders are evaluated through both guantitative analysis (box plots) and qualitative examination of the decoded traces. Each trial within the NPI-NHP recording has a duration of 9.85 seconds. Our method outperforms all other decoders on predicting the arm force.
Different brain regions from 20 datasets. The decoding results across various brain regions for 20 IBL sessions summarized in Figure Overall, our approach consistently achieves higher R2 values compared to other competing methods in decoding both motion energy and wheel speed across the five recorded brain regions. Notably, decoders based on "good" KS units exhibit poor performance across all recorded brain regions when compared to decoders based on all KS units. This observation highlights the importance of utilizing all available information for decoding behaviors rather than solely relying on "good" units based on sorting quality metrics. The scatter plots in Figure (b) indicate general trend where decoding quality tends to increase when more components (i.e., KS units for the spike-sorted decoders and MoG components for our method) are available for decoding. However, our method outperforms spike-sorted decoders even with limited number of components.
Different probe geometry. The decoding results for the NP2.4 and NPI-NHP geometries are illustrated in Figure For NP2.4, our approach significantly outperforms other competing methods when decoding motion energy and wheel speed, while again, the approaches are more comparable when decoding the discrete choice variable (our method performs slightly worse than the spike- sorted decoder). For NP1-NHP, Figure 4 demonstrates that our method achieves better decoding performance (R² z 0.6) compared to the thres--hreshdddd (R² z 0.2) and spike-sorted baselines (R² s 0.3). "Good" KS units are not available in this scenario (the IBL quality criteria were not applied to this primate dataset) and therefore not included in the results.
Comparison to state-of-the-art clusterless decoder. We compare our method to a state-of-the- clusterless point process decoder (Denovellis al. 202T) in Table T Our method has higher decoding performance than the point process decoder on both HD and simulated tetrode datasets for all behavior variables. This performance improvement likely due to the increased flexibility of our decoder compared to state-space models that make stronger assumptions about the dynamics of the decoded signals.
Computation time. In Figure 5 we provide computation time comparison relative to real-time. Our decoding step operates at a sub-real-time pace (0.3 times real-time). The total time after preprocessing for our method is close to real-time.

---

Multiple tetrodes (position)
NPI (wheel speed)
NPI (motion energy)
Denovellis et 2222))
0.91 (= 0.01)
0.50 (# 0.16)
0.55 (# 0.15)
Density-Based
0.97 (= 0.03)
0.63 (#: 0.12)
0.63 (=: 0.14)
Table I: Comparison to a state-of-the-art clusterless decoder. evaluated the performance of both methods using 5-fold cross-validation. We reported the mean correlation between the ground- truth behavior and the decoded behavior along with the standard deviation. All tetrode data was simulated. For the HD datasets, averaged the results across three IBL datasets.
time (real-time)
1.5
0.9
0.6
proprocessing spıi subtraction donolinng local
roglstration
danatty docoding
arcprecess
Figure 5: Computation time measured relative real-time. "Preprocessing" includes destriping, required by all decoders (IBI. et al., 2022). "Total after preprocess" includes spike subtraction, denoising, localization, registration and density-decoding. The computation time of the clusterless point process decoder (Denovellis et al. 2021 is also provided.
5 Discussion
In this work, we introduce a probabilistic model-based neural decoding method that relates spike feature distributions behavioral correlates for more accurate behavior decoding. Our method designed for high-density recording devices such as Neuropixels probes, utilizing novel HD spike features (i.e., spike locations) and maximum peak-to-pesd seap amplitudes of the spikes. We further develop an efficient variational approach to perform inference in this model. We benchmark method across comprehensive set of HD recordings with varying levels of sorting quality, different probe geometries, and distinct brain regions. We demonstrate that our decoding method can consistently outperform spike-thresholded decoders and spike-sorted decoders across wide variety of experimental contexts. This motivates a shift towards a spike-feature based decoding paradigm that avoids the need for spike sorting while also achieving comparable superior decoding performance to approaches relying on well-isolated single-units.
While our method shows promising results, essential to explore avenues for improvement. Two potential improvements
method include utilizing deep learning-based models to capture complex functional associations between firing rates and behaviors and also introducing dependencies among the mixture components account for correlated neuronal firing patterns. An interesting extension of this work would be to apply this spike-feature based paradigm to unsupervised learning of neural dynamics which would enable estimate MoG "firing rates" conditioned on time-varying latent variables. By addressing these challenges and expanding the scope of our method, we can advance our understanding of neural activity and its relationship with behavior in both the supervised and unsupervised settings.
Broader Impact
Neuroscience research aims to uncover how much behavior information can decoded from neural signals in specific brain regions. To do this fairly across regions without bias from spike sorting quality, we need decoding methods that do not rely on spike sorting. However, we should recognize the current decoders' limitations to avoid drawing erroneous conclusions from the decoding results. Our generative model, which improves predictive performance by conditioning on external variables, is also relevant to the machine learning community.
10

---

Acknowledgement
This work supported by grants from the Wellcome Trust (209558 and 216324), National Institutes of Health (IU19NS123716 and K99MH128772) and the Simons Foundation. We thank Matt Whiteway, Tatiana Engel, and Alexandre Pouget for helpful conversations and feedback the manuscript.
References
David M Blei, Alp Kucukelbir, and Jon D McAuliffe. Variational inference: A review for statisticians.
Journal of the American statistical Association, 112(518):59--77, 2017.
Julien Boussard, Erdem Varol, Hyun Dong Lee, Nishchal Dethe, and Liam Paninski. Three- dimensional spike localization and improved motion correction for neuropixels recordings. Ad- vances in Neural Information Processing Systems, 14:22095-221088 2021.
Julien Boussard, Charlie Windolf, Cole Hurwitz, Hyun Dong Lee, Han Yu, Olivier Winter, and Liam Paninski. Dartsort: A modular drift tracking spike sorter for high-density multi-electrode probes. bioRxiv, page 553023, 2023.
Alessio Buccino, Cole L Hurwitz, Samuel Garcia, Jeremy Magland, Joshua H Siegle, Roger Hurwitz, and Matthias H Hennig. Spikeinterface, unified framework for spike sorting. Elife, e61834, 2020.
Alessio Paolo Buccino, Samuel Garcia, and Pierre Yger. Spike sorting: new trends and challenges of the era of high- probes. Progress in Biomedical Engineering, 2022.
Mayo Faulkner Chapuis, Kenneth D Harris, Julia M Huntenburg, Cole Hurwitz, Hyun Dong Lee, Liam Paninski, Cyrille Rossant, Noam Roth, Nicholas A Steinmetz, Charlie Windolf, et al. Spike sorting pipeline for the international brain laboratory. chunnels, 10:6, 2022.
Zhe Chen, Fabian Kloosterman, Stuart Layton, and Matthew A Wilson. Transductive neural decoding for unsorted neuronal spikes of rat hippocampus. pages 1310-1313, 2012.
Xinyi Deng, Daniel F Liu, Kenneth Kay, Loren M Frank, and Uri T Eden. Clusterless decoding of position from multiunit activity using marked point process filter. Neural computation, 27(7): 1438 1460, 2015.
Eric L Denovellis, Anna K Gillespie, Michael E Coulter, Marielena Sosa, Jason E Chung, Uri T Eden, and Loren M Frank. Hippocampal replay of experience at real-world speeds. Elife, 10:e64505, 2021.
George W Fraser, Steven M Chase, Andrew Whitford, and Andrew B Schwartz. Control of brain-computer interface without spike sorting. Journal of neural engineering, 6(5x055004, 2009.
Joshua Glaser, Ari S Benjamin, Raced H Chowdhury, Matthew G Perich, Lee Miller, and Konrad Kording. Machine learning for neural decoding. Eneuro, 7(4), 2020.
Gerrit Hilgen, Martino Sorbaro, Sahar Pirmoradian, Jens-Oliver Muthmann, Ibolya Edit Kepiro, Si- mona Ullo, Cesar Juarez Ramirez, Albert Puente Encinas, Alessandro Maccione, Luca Berdondini. et al. Unsupervised spike sorting for large-scale, high-density multielectrode arrays. Cell reports, 18(10):2511-2532 2017.
Daniel N Hill, Samar B Mehta, and David Kleinfeld. Quality metrics to accompany spike sorting of extracellular signals. Journal of Neuroscience, 31((4):8699 8705, 2011.
Cole Hurwitz, Kai Xu, Akash Srivastava, Alessio Buccino, and Matthias Hennig. Scalable spike source localization in extracellular recordings using amortized variational inference. Advances in Neural Information Processing Systems, 32, 2019.
International Brain Laboratory IBL, Kush Banga, Julius Benson, Niccolò Bonacchi, Sebastian Brui- jns, Rob Campbell, Gaëlle A Chapuis, Anne K Churchland, M Felicia Davatolhagh, Hyun Dong Lee, et al. Reproducibility of in-vivo electrophysiological measurements mice. bioRxiv, pages 2022 05, 2022.
11

---

The International Brain IBL, Valeria Souillon-Rodriguzz, Dora Angelaki, Hannah Bayer, Niccolo Bonacchi, Matteo Carandini, Fanny Cazettes, Gaelle Chapuis, Anne K Churchland, Yang Dan, et al. Standardized and reproducible measurement of decision-making in mice. Elife, 10, 2021.
James J Jun, Nicholas A Steinmetz, Joshua H Siegle, Daniel J Denman, Marius Bauza, Brian Barbarits, Albert K Lee, Costas Anastassiou, Alexandru Andrei, Çağatay Aydm, et Fully integrated silicon probes for high-density recording of neural activity, Nature, 55 (7667)))2332-22 2017.
Fabian Kloosterman, Stuart P Layton, Zhe Chen, and Matthew A Wilson. Bayesian decoding using unsorted spikes in the rat hippocampus. Journal of neurophysiolog, 2014.
Alp Kucukelbir, Dustin Tran, Rajesh Ranganath, Andrew Gelman, and David M Blei. Automatic differentiation variational inference. Journal of machine learning research, 2017.
Jesse A Livezey and Joshua Glaser. Deep learning approaches for neural decoding across architec- tures and recording modalities. Briefings in bioinformatics, 2272):1577-1991, 2021.
Jeremy F Magland and Alex H Barnett. Unimodal clustering using isotonic regression: Iso-split.
arXiv preprint arXiv:15088 04841, 2015.
Najja Marshall, Joshua Glaser, Eric M Trautmann, Elom A Amematsro, Sean M Perkins, Michael N Shadlen, LF Abbott, John Cunningham, and Mark M Churchland. Flexible neural control of motor units. Nature neuroscience, 25(11):4422-150, 2022.
Elon Musk et al. An integrated brain-machine interface platform with thousands of channels. Journal of medical Internet research, 21(10):e16194, 2019.
Marius Pachitariu, Nicholas Steinmetz, Shabnam Kadir, Mattco Carandini, et al. Kilosort: realtime spike-sorting for extracellular electrophysiology with hundreds of channels. BioRxiv, page 061481, 2016.
Marius Pachitariu, Shashwat Sridhar, and Carsen Stringer. Solving the spike sorting problem with kilosort. bioRxiv, pages 2023-01, 2023.
Angelique C Paulk, Yoav Kfir, Arjun Khanna, Martina L Mustroph, Eric M Trautmann, Dan J Soper, Sergey D Stavisky, Marleen Welkenhuysen, Barundeb Dutta, Krishna V Shenoy, et al. Large-scale neural recordings with single neuron resolution using neuropixels probes in human cortex. Nature Neuroscience, 25(2):252-266 223,,222-2221 2022.
Mohammad Reza Rezaei, Kensuke Arai, Loren M Frank, Uri T Eden, and Ali Yousefi. Real-time point process filter for multidimensional decoding problems using mixture models. Journal of neuroscience methwds, 348: 109006, 2021.
Nicholas A Steinmetz, Christof Koch, Kenneth D Harris, and Matteo Carandini. Challenges and opportunities for large-scale electrophysiology with neuropixels probes. Current opinion in neurobiology, 50:92-1000 2018.
Nicholas A Steinmetz, Cagatay Aydin, Anna Lebedeva, Michael Okun, Marius Pachitariu, Marius Bauza, Maxime Beau, Jai Bhagat, Claudia Böhm, Martijn Broux, et al. Neuropixels 2.0: A miniaturized high-density probe for stable, long-term brain recordings. Science, 372(6539): eabf4588, 2021.
Sonia Todorova, Patrick Sadtler, Aaron Batista, Steven Chase, and Valérie Ventura. To sort or not to sort: the impact of spike-sorting on neural decoding performance. Journal of neural engineering, 11(5):056005, 2014.
Eric M Trautmann, Sergey D Stavisky, Subhaneil Lahiri, Katherine C Ames, Matthew T Kaufman, Daniel O'Shea, Saurabh Vyas, Xulu Sun, Stephen Ryu, Surya Ganguli, estimation of neural population dynamics without spike sorting. Neuron, 103(2):292-308 2019.
Accurate
Eric M Trautmann, Janis K Hesse, Gabriel M Stine, Ruobing Xia, Shude Zhu, Daniel J O Shea Bill Karsh, Jennifer Colonell, Frank F Lanfrancchi Saurabh Vyas, et al. Large-scale brain-wide neural recording in nonhuman primates. bioRxiv, pages 2023-02, 2023.
12

---

Anne E Urai, Brent Doiron, Andrew M Leifer, and Anne K Churchland. Large-scale neural recordings call for new insights link brain and behavior. Nature neuroscience, 25(1):11-19,22) 2022.
Valérie Ventura. Spike train decoding without spike sorting. Neural computation, 20(4):923-963, 2008.
Valérie Ventura. Traditional waveform based spike sorting yields biased rate code estimates. Pro- ceedings of the National Academy of Sciences, 106(()):6921-66 6926, 2009.
Charlie Windolf, Angelique Paulk, Yoav Kfir, Eric Trautmann, Samuel Garcia, Domokos Meszéna, William Muñoz, Richard Hardstone, Irene Caprara, Mohsen Jamali, et Robust online multiband drift estimation in electrophysiology data. bioRxiv, pages 2022-12, 2022.
13

---

Sitk Supplementary Material
1 Decoding using automatic differentiation variational inference (ADVI)
In the method section of our paper, we describe the general secoding-decodiigg paradigm. In this section of the supplementary material, delve into specific case that focuses decoding continuous behaviors, 3/k E IR"", that exhibit temporal variations. We introduce the use of ADVI, which allows us model the relationship between firing rates of MoG components and behavior correlates using a generalized linear model (GLM).
Notation
Definition
spike feature at time bin trial k: behavior in trial & MoG assignment of spike feature at time bin t in trial & mean covariance matrix of MoG component MoG mixing proportion at time hin t trial & number of spikes collected at time bin in trial & number of MoG components number time bins within each trial number of trials in a session
Zisk
Ge = (Me, xae)
Tak
Yeck
C
K
Aatk
behavior-depender firing rate of component c at time bin & in trial "baseline firing rate" of component e "behavior-man modulated Thchaviormmouulated firing rate" of component time bin mean and variance of the variational posterior distribution for b mean variance of the variational posterior distribution for Aa mean and variance of the variational posterior distribution for yike
be
ADVI
Bet
Ab. = (1hn. Tbe
nA. (p Bee
Pui (Mora yes
dell, Acto
firing rate of component time bin that switches between two states variational paramter that represents the probability of yı unnormalized posterior probability of assigning Sы to component normalized posterior probability of assigning unnormalized posterior probability of yı and yı normalized posterior probability of yı =
CAVI
Pictk
Tictk
to component
Vhr V&o
Pk
Table 2: Table of notation.
1.1 ADVI-based encoder
Building upon the general model specification outlined in Equations 13 we can describe the generative model of the encoder as follows:
Aeth
exp(b. Bat yrk), b. ~ Nbbei 0, 1), Box N(3ez;
(6)
Notk
Zitk ~ Categorical Zank Tek), Fekk (Teekh )o Terk Ee Aerth"
(7)
Satk ~ (Ssuii Beae), ne = (Mc, Ee).
(8)
where and are the unknown model parameters corresponding to in Equation and are sampled from standard normal prior distribution. The latent spike assignment 2 depends on the mixing proportion of the MoG, which is influenced by the behavior y through the firing rate X. Intuitively, we can interpret Aat as the "firing rate" of component at time t in trial k while b. and 3. describe the "baseline firing rate" and "behavior- moduulae-dddullteed firing rate" of component c, respectively.
To learn the latent variables b and B, we posit mean-field Gaussian variational approximation
g(2, B)
II (b2
for the posterior distribution p(2, b, B ss s.y)- Obtaining exact updates for and 3 is challenging due to the normalization term for π in Equation 7 Therefore, employ ADVI to maximize the
14

---

ELBO and utilize stochastic gradient ascent to update and 3. ADVI requires that the model be differentiable with respect the parameters, and with normal prior, the latent variables b and reside in the real coordinate space and cause issues with differentiability. The variational approximations for and are
(()) II glbbs The = II N (b.) Tho = (Abee ob)c
(10)
((3 II (Bez TB.)
II c,t
(11)
drawback of the parameterization in Equations is that the spike assignment variables are discrete and not compatible with ADVI, An alternative, equivalent parameterization that addresses these problems is to marginalize over 2. The marginalized model is
Autk = exp(be + Ber y(a), be ~ Nbbe 0) Ae ~ NN 0, 1),
(12)
dath
Totk =
Ee Awek"
Tek = [Tock
(13)
Task MassNNsssss De (Mer Ee).
(14)
Under this parameterization, the ELBO for the encoder is
cnc == Re((, (l,,pps b y)] E, [log g(b, [3]]
(15)
E (b,
Σ
log.N(sau: ne) log Tark log (bb,,, ((Bzs 0, 1)]
(16)
A,t,c,i
- -gg(,,3)
Σ
log (bos mo.) + (Bee: no.))
c,t
1.2 ADVI-based decoder
The decoder adopts the same generative model as the encoder described in Equations with two exceptions: 1) The latent variable yıkk is assumed to have a standard normal prior, i.c., you ~ (yra; 0, 1), assuming independence at each time step Alternatively, Gaussian process prior can be chosen to capture temporal correlations between time steps. 2) The parameters and learned from the ADVI-based encoder and kept fixed in the decoder.
The mean-field Gaussian variational approximation for the posterior distribution p(2,3 s) is
q(2, 3) II 9(=a)9(31). e,
(17)
where
((v) II ( 1p.)) R
Il N( Agn.), Myn (Huya o k.t
9t )
(18)
To enable the use of ADVI, we can marginalize out the discrete latent variable 2, thereby transforming the MoG model into differentiable form. Under the marginalized MoG parametrization in Equations 12 14 the ELBO of the decoder is
ADVI dec
: Eg(w) [log (og - By(n 4(v)]
(19)
Ey(Co)
Σ
log W me) + log Firk + logNN (yux; 0, 1)]
(20)
k.t.c.&
[gg(t)
NogN(vrss )
k.
15

---

2 Decoding using coordinate ascent variational inference (CAVI)
We present specific scenario for decoding binary variables, ys E {0,1), 1), where we derive exact updates for the variational variables using CAVI.
2.1 CAVI-based encoder
Extending the general model described in Equations defined as follows:
13 the generative model the encoder can be
Aet1
U2
Aeto
1-ve
cetk Ee
den
Ev derio
(21)
Fitk
CCategoricalFstt mekk; Teek {Teth I (Sitt NENr) = (Mes Me),
(22)
Sick
(23)
where z and À are the latent variables that we aim learn. The behavior-dependen firing rates of each component at time f vary based on the binary variable y, such that the components switch between two behavioral states characterized by firing rates Act and Aao-
The log- lelelikeiiood of the encoder can be written
logp(8, 11) = Σ
Zzetk { ((suii + yn (log Aae - -log Ar1)
(24)
k.t.e,t
+ - yn (log An1o log (ea)
where = Aen and = A To approximate the posterior p(= 8, y), we employ the mean-field variational approximation ((z) = Ne, ((=xt). The ELBO of the encoder is
CCAVI enc
= loggpsss b)] - (E q()))-
(25)
The exact update for g(2) is obtained by maximizing the ELBO with respect to g(z), which leads to the following update equation:
(Fxx) x Exp(E,, (loopp(za 2-0t, y
(26)
where q-ot means
Ie He ((20e). Then.
qlZict = 1) exp { log NogNNsaus ne) + 3r log
Aeel
)
+ (1
ykd log
Pictk (27)
Ar1
Aro
denotes the unnormalized posterior probability of assigning spike collected at time t in trial to component c, while E[zick] Pictk/ Eet Pieltk Fictk represents the normalized posterior probability.
After fixing q(z), the term in the ELBO which depends on X, μ and E can be expressed as
(ottl
(ee20)
L
Σ k.t.o,
Tieth log log (suui Me)
yı log
+ (1 ykk log
}.
An
Awo
16

---

Algorithm CAVI-based encoder
Input: {sakk: i 1 Rek, 1::.t=1:1p 1
1: K, number of components C.
Initialize {Jey Ec)e C while ELBO not converged do K do for all & : do for all i E : Tikk do
for all
Set [Zicck = Pictk-
eq. 27
end for end for end for
for all c E C do
Set Mr. 12 Eg = E.-
eq.
33 39
for all / 1::::: do
Set Aeto New Art
Xan
eq. 32
33
end for
end for
Compute the ELBO CCAVT end while
eq. 25
Return (22). X, He E.
We derive the update for À by setting the gradients V VAL. L and V VAcio L to 0:
V Vetu L- V ren Σ
Tiotk
yn (log doti - log Ae1)
(28)
Enp Fietk yk
Exs Ee Teelth
yk
(29)
Aet
Ae1
Σ Fictk yk (har
Eeer)) = Σ
Σ
Fieltk Us den
(30)
c'scc
d=l
Σ yk Fiett Piett)ad =
Σ
M/
Tietk Σ
de'e))
(31)
c'
Aid
o'sc
Exi Mk rioch ( (Eotso ddres)
Ad1
== Nox1
(32)
Σki yk (Eeffee Tecrek)
Ex,(1 - Mkrrictte EEd dera) Σx,(1 - -ykk (Ee Poeres)
Deto 1
Ada
= Nan-
(33)
Consider the gradient with respect to the no parameter,
V 1o L
Σ Pietk Vn, Iog.N(sass ne)
(34)
Σ
Fietk V Vno 1
(log |3-5 Tr(E, (sick - Mc)(aak Mc)
(35)
17

---

The closed-fone updates for Me and E. are
Preth 2 (Sith Mc) =0
(36)
= He
8atk = Mor Be = Σ
Tictk-
(37)
Re k,t,i
Vx L
Σ
Tieth (Ee (#tth - M5se - [e) =
(38)
k
Tieth (8ith Me) Mc) = De
(39)
The R.,
CAVI-based decoder
The CAVI-based decoder employs the same generative model the CAVI-based encoder, with exception that the behavior-dependent firing rates Aot and Aæo are learned by the encoder and kept fixed, and the behavior is treated as unknown latent variable. We sample from prior distribution, yk (ernoulli@)), where & is a variational parameter that represents the probability that y 1. The log-likelihood can be expressed as follows
log 2, = Σ k,t,o,
Zietk { (oo.. (ssuxit) + yu Aee - log Ac1)
(40)
(1 4k (log Aen log (oo) + log Φ (1 yk) log(1 d).
We use the factorization q(2, y) = (() II ((<)))(V)) to approximate the posterior distribu- tion p(2,y) s). The ELBO of the decoder can be defined as
CAVI dec
E (=.0)) [log 2, y)] E, [log ((2, v)].
(41)
Algorithm 2 CAVI-based decoder
Input: {s }. (Acto: dct1); ): i = 1 Rek, t t= 1: k=
Initialize {pc}. {Ec}. while ELBO not converged do for all & E 1 K do
Set q(vk =)
Ver-
eq.
44
for all t E 11:78 To do
for all i : do
Set (ZZirk
x Pictk-
eq. 43
end for
end for
end for
for all CE 1: do
Set Me 112 Eg = Ed
eq.
(36399 39
end for Set 6²
eq.
46
Compute the ELBO end while
CCAVI dec
D eq.
Return ((2,)), 6, E.
The exact updates for q(x) and qιy) that guarantee an increase in the ELBO are
((2) exp [og
q(v) (og
(42)
where
(Zziot = 1) x Noalli + E[yna] log (2t))) + (1 - E[yz]) log
(e A.
}} = Pietks
(43)
18

---

and E[Zickk] = Psctk/ Ee Pic'tk = Tietk are the unnormalized and normalized posterior probabilities of assigning spike collected at time in trial & component c, respectively. The unnormalized posterior probabilities of and ys 0 are
Qyv: = 1) xc expl
El Eiziol (log Ae1 log Ar1) + log o} := yk1
(44)
duk = 0) x expl) t,c,
Σ
[[ziors](opp X log Aro) + log(1 - -6} :- -yko-
and E[ya] / + vūo) = Vi represents the normalized posterior probability of yı 1.
After fixing (2) and q(v), the term in the ELBO which depends o, μ and E can be written
L
Σ Fieth (log N (ssu; me) + uy (log dat log Ac1)
(45)
K,t,,,i
+1 (log Acto log A + u log δ (1 - vx) log(1 Φ).
Considering the gradient of the ELBO with respect to the o parameter, we obtain its update:
L
Vk-
(46)
The updates for Tee are computed in similar manner described in Equations 36 39
3 MoG initialization
We employ the following procedure to intialize the MoG model used in both the ADVI-base and CAVI-based algorithms: 1) According to Figure (a), the spike feature distribution is highly multimodal. To determine the appropriate number of modes, we utilize isosplit (Magland and Barnett 2015) to cluster the spike features. This step helps in splitting the set of spike features into distinct clusters. 2) For each identified cluster, we compute the mean and variance of the spike features belonging to that cluster, which serve as the parameters for the corresponding Gaussian component. This automatic selection of the number of MoG components and the initialization of means and covariance matrices facilitate the initialization of the ADVI-based and CAVI-based algorithms.
19

---

Data preprocessing
We provide brief overview of our data preprocessing pipeline, which involves the following steps.
Destriping. During the data collection process, we encounter line noise due to voltage leakage on the probe. This translates into large "stripes" of noise spanning the whole probe. To mitigate the impact of these noise artifacts, we apply destriping procedure (Chapuis et al. 2022).
Subtraction-bast spike detection and denoising. We employ the iterative subtraction-bact procedure for spike detection and collision-correction described in Boussard et al. (2023).
Spike localization. each denoised spike.
We employ the method of Boussard al. 2021 estimate the location of
Drift registration. Probe motion (or drift) in the electrophysiology data poses challenge for downstream analyses. Decentralized registration (Windolf et al. (0022) is applied to track and correct for motion drift in the high-density probe recordings.
a)
"good' sorting example
"bad' sorting example
3500
3000
depth (um)
2500
2000
1000
500
50
100
150 time (s)
200
250
50
100
150
200
250
time
Figure 6: Motion drift in "good" and "bad" sorting recordings. (a) The motion-registered spike raster plot of a "good" sorting example that is less affected by drift. (b) The spike raster plot of "bad" sorting example, which is still affected by drift even after registration.
Behavior decoder
To decode binary behaviors, such the mouse's left right choices, we utilize Le-penalized logistic regression. For decoding dynamic behaviors, such as wheel speed, we employ sliding- algorithm aggregate the entries of the weight matrix. Wkats over time. Within the time window [ (t-& 5, f + 5). where 5 is the window size, we stack 28 weight matrix entries, (Wte) ) for time point in trial kt This aggregated weight matrix is then used as input for ridge regression to predict the behavior Vik at time t. The window size 5 and the regularization strength are model hyper-parametern set through cross-validation to achieve the optimal decoding performance.
20

---

Model interpretation
In this section, we provide visualizations to gain insights into the effectiveness our proposed decoder. We quantify the posterior entropy of each spike assignment in Figure (a). Spike assignments with higher entropy correspond
spread of posterior probabilities among multiple MoG components. In contrast, traditional spike sorting or thresholding methods result in deterministic spike assignments, leading to lower entropy and empirically reduced decoding performance. Figure we compare the trial- trial-averaged weight matrices (W) used for decoding between spike-sorted, spike-thresholdd, and our proposed decoders.
b)
W (spikossrilac)
W.' (spko-trreshlleed
W (curs)
War
2.00
10-
3.75
component
20
1.00
0.75
30-
40
0.50
0.25
0.3
Figure 7: Model interpretation. (a) The posterior entropy of the spike assignment is high, when the posterior probability of spike assignment, ((Zsket). is spread out among several MoG components instead being concentrated at a single component. The scatter plot shows that low-ttplt spikes that are difficult to assign have higher posterior entropy than high-amplitude spikes. (b) Visualizations of the averaged weight matrices W's across trials in an example IBL session. For "W (spike-sorted and spike-thresholded)", the W matrix has one-hot rows and each entry Wk is the spike count that belongs to KS unit (channel) and time bin in trial The purple crosses indicate the "good" KS units. For "W (ours)", Wket is the sum of posterior probabilities of spike assignments, and the MoG mixing proportion π depends on the behavior y and changes over time. The arrangement of KS units (channels or MoG components) the heat maps is based on the depth of the NP probe, ensuring comparability across the displayed matrices.
7 Decoding across brain regions
In addition to the previously mentioned five brain regions (PO, LP, DG, CAI, VISa) depicted in Figure 3] we expanded our analysis to include two additional brain regions situated in the cerebellum: the arbor vitae (ARB) and the ansiform cruciform lobule (ANCR). We specifically include the cerebellum in our analysis due to the frequent occurrence of spike sorting issues in this area. In Figure and we present a comparison between our decoder and the spike-sorted decoder that utilizes all KS units across all the brain regions studied. According to Figure g our method consistently outperforms the decoder that relies on all KS units across all brain regions and the majority of IBL. sessions. Furthermore, Figure 9 specifically demonstrates that method consistently achieves superior decoding performance in the recorded regions of the cerebellum, where spike sorting quality issues are commonly encountered. This highlights the robustness and reliability of our method, particularly in challenging recording conditions.
21

---

All
molion energy
wheel speed
choice
0.8
0.6
accuracy
our R2
0.6
our R2
0.4
0.4
0.2
0.2
0.0
0.2
0.4
0.6
0.8
0.0
0.2
0.4
0.6
0.6
0.8
1.0
KS
R²
KS
R²
KS accuracy
PO
accuracy
0.6
our R2
our
0.4
0.4
0.2
0.2
0.6
0.0
0.0
0.4
0.2
0.4
1.0
R²
R²
KS accuracy
LP
0.6
0.5
accuracy
0.9
our R²
our A2
0.4
0.3
0.7
0.2
0.1
0.5
0.0
-0.
0.3
0.0 0.2 0.4 0.6
-0.1
0.1
0.3
0.5
0.3
0.5 0.7 0.9
R²
R²
KS accuracy
DG
accuracy
0.6
0.4
our R2
our R2
0.7
0.4
0.2
0.2
0.5
0.0
0.0
0.0
0.2
0.4 0.6
0.0
0.2
0.4
0.3
0.5
0.7
0.9
KS R²
KS R²
KS accuracy
CA1
accuracy
our
0.6
Rur
0.4
0.4
0.2
0.5
0.0
0.0
0.3
0.2 0.4
0.3
0.5
0.7
R²
R²
KS accuracy
VISa
accuracy
0.5
our R2
0.4
0.9
our
0.4 X
0.2
0.2
0.0
0.5
0.0
0.2
0.4 0.6
0.0
0.2 0.4
0.5
0.7
0.9
1.1
KS R?
R
KS accuracy
Figure 8: Decoding comparisons across brain regions in the thalamus, hippocampus and visual cortex. We compare our decoder to the spike-sorted decoder using all KS units across various brain regions and behavioral tasks. Each point in the scatter plot represents one session from the 20 IBL session previously described in the experiments section. Sessions with "good" sorting quality are depicted in green, while sessions with "bad" sorting quality are shown in purple. The majority of the sessions lie above the gray diagonal line, indicating that our method consistently outperforms the decoder relying on all KS units.
22

---

All
motion energy
wheel speed
choice
0.7
our accuracy
1.1
R²
0.6
our R2
0.5
our
x
x
0.9
0.4
0.3
0.2
0.1
0.7
0.2
0.4
0.6
0.1
0.3
0.5
0.7
0.7
0.9
1.1
KS
R²
KS R²
KS accuracy
ARB
motion energy
wheel speed
choice
0.6
accuracy
OUr R2
0.4
0.6
XXX
++
our R2
0.4
0.7
0.9
0.2
0.2
0.5
0.0
0.0
0.3
0.0 0.2 0.4 0.6
0.0
0.2
0.4
0.6
0.3
0.5 0.7 0.9
KS R²
KS R²
KS accuracy
ANCR
motion energy
wheel speed
choice
0.5
our accuracy
1.1
0.6
Our R²
our R2
X
x
0.4
0.3
0.9
0.2
0.7
0.2
0.4
0.6
0.1
0.3
0.5
0.7
0.9 KS accuracy
1.1
KS R²
KS R
Figure 9: Decoding comparisons across brain regions the cerebellum. We evaluate our decoder against the spike-sorted decoder utilizing all KS units across various brain regions in the cerebellum. The scatter plot visualizes the results for each IBL session in the study, with all sessions depicted in purple indicating "bad" sorting quality. Note that the probes used in the cerebellum sessions were not uniformly implanted in the same set of brain regions. As result, the number of available sessions for decoding in certain cerebellum regions is limited. Notably, the majority of sessions are positioned above the gray diagonal line, indicating that our method consistently achieves better decoding performance compared the decoder relying on all KS units. This outcome highlights the robustness of our method, particularly in the context of the cerebellum where spike sorting quality issues are prevalent.
23

---

8
Ablation study
We conduct an ablation study to investigate the importance of various components in our decoding paradigm. Specifically, we examine how the integration of dynamic beeaavor-deepedder mixing proportions in the MoG and figher-dimensional spike features influence decoding performance. Additionally, we analyze how the inclusion criteria for spike waveforms can affect decoding outcomes.
Effects of dynamic MoG mixing proportion. Table B presents a comparison between the ordinary MoG with fixed mixing proportion (referred to as "fixed T²" and our proposed model with dynamic mixing proportion (referred to as "dynamic r²") The results indicate that this approach leads to improved decoding performance compared to using fixed mixing proportion in the MoG model.
Motion energy (It-)
Wheel speed (R²)
Choice (accuracy)
Fixed x
Dynamic
Fixed
Dynamic IT
Fixed
Dynamic IT
All
0.664 (+ 0.034)
0.742 (+ 0.028)
0.470 (+ 0.062)
0.564 (+ 0.045)
0.948 (+ 0.038)
0.957 (+ 0.036)
PO
0.365 (# 0.031)
0.488 (= 0.046)
0.520 (# 0.019)
0.670 (# 0.015)
0.844 (# 0.016)
0,861 (# 0.035)
LP
0.145 (+: 0.015)
0.464 (= 0.054)
0.114 (+ 0.032)
0.342 (+ 0.027)
0.917 (+ 0.026)
0.931 (+ 0.022)
DG
0.280 (1: 0.033)
0.492 (= 0.042)
0.221 (#: 0.040)
0,381 (#: 0.035)
0.669 (1 0.084)
0.722 (1: 0.050)
CA1
0.407 (#: 0.021)
0.538 (= 0.038)
0,308 (#: 0.051)
0.428 (#: 0.030)
0.621 (# 0.062)
0.626 (# 0.064)
VISa
0,488 (.L 0.046)
0.490 (1 0.047)
0.318 (1 0.056)
0.364 (1 0.074)
0.857 (1 0.065)
0.874 (L 0.066)
Table 3: Effects of dynamic mixing proportion of the MoG on decoding performance.
Effects of higher-dimensiinall spike features. In Table 4 we provide a comparison of decoding performance using two different sets of spike features. The first set includes spike location along the width and depth dimensions of the probe (denoted as x and 2) as well as the maximum peak-to- peak amplitude of the spike (denoted as a). The second set includes the first and second principal components (PCs) of the spike waveforms (denoted as 41 and U2) in addition to x, z and a. The spike features are visually represented using scatter plots in Figure I0 We report the mean and standard deviation of the decoding accuracy (It²) obtained from a 5-fold CV on a single session. Table 4 provides insights regarding the inclusion of additional waveform PC features for decoding. The findings suggest that the incorporation of these additional PC features does not contribute to significant improvements in decoding performance.
Motion energy (IP2)
Wheel speed (RP")
Choice (accuracy)
(2, =.
(,,,,0, 01, Mg]
(2)
[r, 11, U2)
(1)
(2, M1, M1, U2)
All
0.531 0.026)
0.529 0.026)
0.484 (+ 0.042)
0.478 (+ 0.045)
0.917 (+ 0.019)
0.917 (+ 0.019)
PO
0.462 (1: 0.038)
0.457 0.039)
0.469 (1: 0.061)
0.464 (1: 0.051)
0.853 0.019
0.853 (1: 0.025)
LP
0.490 (+: 0.026)
0.489 (= 0.029)
0.479 (# 0.035)
0.473 (# 0.013)
0,864 (# 0.022)
0.849 (#: 0.036)
DG
0.335 (.L 0.028)
0.321 (1 0.037)
0.273 (1 0.030)
0.264 (11 0.025)
0.675 (1 0.038)
0.679 (L 0.049)
CA1
0.449 (1: 0.027)
0.440 (1: 0.046)
0.329 (1: 0.044)
0.328 (1: 0.040)
0.755 (1: 0.045)
0.758 (1: 0.053)
VISa
0,270 (+ 0.021)
0.237 (+ 0.017)
0,225 (+: 0.024)
0.206 (t 0.013)
0.725 (+ 0.048)
0.732 (t: 0.054)
Table 4: Effects of incorporating higher- dimensiaam spike features on decoding performance.
24

---

4000
3500
3000
10
2500
maxpppp amplltude (s.u.)
(um)
2000
1500
70
1000
500
25
96
-100
100
40
50
100
-25
25
x (um)
maxptp amplitude (s.u.)
U1
U2
Figure 10: Visualizations of spike features employed for decoding. Spike localization features, (x,, 2), the locations spikes along the width and depth of the NP1 probe, and waveform features, the maximum peak-to-peak (max ptp) amplitudes of spikes. Amplitude is measured in standard units (s.u.). u and Ug denote the first and second principal components (PCs) of the spike waveforms.
a,
25

---

Effects of inclusion criteria for spike waveforms. To investigate whether the density- bassd decoder is performing better by using additional spikes that a spike sorter would miss, we conducted an experiment using different inclusion criteria for spike waveforms. We fitted our model using only spikes detected by Kilosort 2.5, and compared its performance to decoders using spike-sorted outputs and our subtraction- based spike detection on choice and motion energy decoding. The results are summarized in Table As shown in the table, our decoder can achieve comparable or better decoding performance than the spike-sorted decoder when modeling the same spikes, This suggests that the gain in decoding performance can be attributed to the density-based approach, instead of the spike inclusion criteria.
Choice (accuracy)
Motion energy (R²)
Density-bused (subtraction spikes)
0.876 (= 0.068)
0.589 (+: 0.111)
Density-bused (KS spikes)
0.876 (= 0.079)
0.579 (+: 0.121)
Sorted (KS spikes)
0.887 (+ 0.078)
0.503 (1: 0.117)
Table 5: Effects of inclusion criteria for spike waveforms on decoding performance.
9 Comparison to a state-of-the-art clusterless decoder
In this section, we outline the specific experimental setup for evaluating our density-based decoder in comparison with the clusterless point process decoder (Denovellis et al. 2021 on both multiple tetrodes and high-density (HD) probes data.
Application to tetrodes. We utilized the code provided in the GitHub repositor point process decoder (Denovellis et al.) (02I) to generate simulated neural and behavioral data. This synthetic dataset was designed to mimic recordings from 5 tetrodes, with each tetrode containing 4 channels. Spike amplitudes from each channel of these multiple tetrodes were selected as the spike features for both decoding methods. The objective of the decoding was to estimate the animal's 's position from these simulated spike features. As the original simulated position was too straightforward to decode, we intentionally distorted it by blending with real position data sourced from the GitHub repository Moreover, we introduced random Gaussian noise to the simulated position for added complexity.
! of the clusterless
Application to HD probes. We evaluated both decoders on decoding wheel speed and motion energy from spike features extracted from NP1 probes across three IBL datasets. To preprocess the data, we followed the pipeline outlined above, extracting spike localization features and maximum peak-to-peak amplitudes as common features for both decoders.
We only focus on continuous behaviors for decoding since the clusterless point process decoder is designed exclusively for continuous behaviors. Notably, due to its continuous time nature, the clusterless point process decoder directly decoded behaviors without time binning. In contrast, our density-based model required time binning of both behavioral and spike feature data into equal-sized time intervals. Consequently, we used the time-binned behaviors for decoding with the density- based approach. We employed 5-fold cross-validation to assess the decoding performance of both decoders. We used a random walk as the state transition for the clusterless point process decoder, and used the estimated variance from the animal's 's behavior to set the variance of the random walk. The clusterless point process decoder uses a grid-based approximation of the inferred behavior, discretizing the behavior space into place bins; see Section "Decoding" in Denovellis al. (202)). We determine place bin size based on the square root of the observed behavior variance. Denovellis et al. (202I) use kernel density estimation (KDE) to estimate the distributions of both the behavior variable and the spike features used for decoding. We set the KDE bandwidth that determines the amount of smoothing done for spike features to be 1.0, and the bandwidth for behavior to be the square root of the observed behavior variance; see Section "Encoding clusterless" in Denovellis et af. (2021).
26

---

10 Simulation for model validation
We conducted simulations to illustrate the principles of method. The simulation aimed show that our encoding model can learn the relationship between spike features and behaviors. We performed two tasks, decoding a binary variable, yk- simulated from a Bernoulli distribution, and decoding continuous variable, yk, simulated from Gaussian process. To mimic the data- setta-ggnerating process, we selected Gaussian components with "templates" extracted from a real dataset. The encoding model parameters, and B, were also taken from learned parameters in the same dataset. Given b. B and yk- we simulated the "firing rates" for each Gaussian component the mixture, as described in the Method section of paper. Next, generated spike features based on these simulated "firing Figure [T] displays the
rates," and applied the encoding model to infer the behavior-dependent learned À for each component c, time t, and trial k. The learned "firing rates" closely resembled the simulated ones, indicating the model's s ability to recover the primary associations between spike features and behaviors. With such associations, the decoding model decode behaviors.
a)
simulated
learned
AdA
b)
simulated
learned
X
0
15
component
10
component
20
20
5
5
40
40
40
1.0
1.0
time (second)
time (second)
time (second)
time (second)
Figure 11: Our encoding model recovers the relationship between the simulated spiking activity and the simulated behavior correlate. Panel (a) shows comparison of the simulated firing rates conditioned on the binary behavior variable with the learned firing rates by our encoding model. In Panel (b), compare the simulated firing rates conditioned the continuous behavior variable with the learned firing rates from our encoding model.
27

---

11 Relationship to spike sorting
We conducted experiments investigate the biological interpretation of our MoG units and the correspondence between single cells identified by KS and our MoG units. The agreement matrix between "hard" KS spike assignments and "soft" MoG assignments is shown in Figure 12 We calculated the conditional probability of spikes belonging to each MoG component, given that these spikes belong to the corresponding KS unit. Notably, KS units with large amplitudes are less likely to be split into multiple Gaussian components. This shows a reasonable correspondence between the Gaussian components and the spike-sorted units.
Prob.
10
0.75
20
MoG units
30
0.00
Ampliudde (s.u.)
30
60
20
10
20
40
60
80
100
120
KS units
Figure 12: Correspondence between Kilosort and MoG spike assignment. Units ordered by their depth on the Neuropixel probe. The color bar shows the conditional probability of spikes belonging each MoG component, given that these spikes belong the corresponding KS unit The mean amplitude of each KS unit is shown at the bottom.
28


Reviews:
[{"review_id": "46J0Lv0Bdg", "reviewer": "Reviewer_tDhN", "summary": "The paper develops a decoding method directly on ‘spike features’, without going through spike sorting. The authors model spike assignment uncertainty using a mixture of Gaussians model, and then perform variational inference to model the relationship between the spike features and behavior.", "questions": "It may be a good idea to apply previously developed clusterless methods on this kind of data to show how / why they fail.\n\nSimulated data will be helpful to the reader to better understand the limitations of the modeling approach.", "limitations": "The limitations are adequately discussed. Potential negative societal impact is adequately discussed.", "rating": 7, "confidence": 4, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "The premise and models in this paper are sound, if not particularly novel. The inference methods are not particularly novel either. However, it is clearly shown that the models developed here are empirically quite powerful. The evaluations performed are very extensive, showing very clearly that their methods outperform standard spike sorting and decoding methods on large datasets in multiple species.", "weaknesses": "It is unclear why exactly previous clusterless methods fail on this kind of data. The authors may need to show this concretely to be fair to previous approaches that do not use spike sorting.\n\nSimulation data may help to get intuition on the methodology, and to show validation of the methods.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "**General:**\nWe thank Reviewer tDhN for carefully reviewing our manuscript, and appreciate the opportunity to address the concerns raised. We want to emphasize the novelty of our method, which involves conditioning of the data generating process on external variables, allowing for improved prediction. Our approach is not limited to GMM modeling and GLM decoding; instead, it can be extended to non-Gaussian mixture models and nonlinear models to tackle diverse modeling challenges. Furthermore, compared to previous clusterless decoders based on state-space models, our density-based decoder avoids making explicit assumptions about underlying system dynamics, and is thus more flexible at capturing complex relationships in the data. Our method is also more scalable than previous approaches, which is desirable for large data volumes generated by HD probes.\n\n**Weakness / Question 1:**\nWe have compared our decoder to previous clusterless decoders on both multiple electrodes and HD probes to demonstrate that our method is effective in both scenarios. The global response Table 1 provides a benchmark of clusterless decoders, reaffirming the superiority of our density-based decoder when applied to HD probes.\n\n**Weakness / Question 2:**\nWe appreciate the reviewer's feedback on conducting a simulation study for model validation. The simulation results shown in global response Figure 1 confirm that our encoding model effectively captures the association between the simulated spike features and the simulated behavior of interest. With such learned associations, the decoding model is able to accurately decode the behaviors. To highlight the significance of these learned associations, we have compared our GMM conditioned on the decoding variable with a vanilla GMM that lacks such informative associations. The details of this comparison can be found in Table 2 in Section 8 (ablation study) of the supplementary materials."}, {"review_id": "XbuDYPFWVR", "reviewer": "Reviewer_Fmkq", "summary": "This paper proposed to decode animal behavior with a spike-sorting-free method that is well-suited for high-density recordings, by modeling the distribution of extracted spike features using a mixture of Gaussians (MoG) on uncertainty of spike assignments, and decoding using a generalized linear model (GLM). The authors benchmarked results on Neuropixel recordings from different brain regions, different animals, as well as different probe geometries, and showed that performance outperformed decoders based on multi-unit threshold crossings and single-units sorted by Kilosort.", "questions": "1. I'm not following the authors' arguments on why previous Bayesian decoding methods on unsorted spikes do not apply for HD probes. Specifically, this argument \"the aforementioned approaches are not suitable for HD probes as they rely on sampling-based Bayesian inference and are specifically designed for tetrodes\". It's true that sampling-based Bayesian inference makes certain assumptions on point process, but I don't quite see why those are specific for tetrodes. The original Bayesian inference paper used waveform features such as the maximum amplitudes and principal components for decoding, but that do not mean waveform features cannot be represented in other formats. Can authors elaborate this statement?\n2. What are the baseline decoders used with sorted neurons after sorted by Kilsort? \n", "limitations": "N/A", "rating": 4, "confidence": 4, "soundness": 3, "presentation": 3, "contribution": 2, "strengths": "1. The paper is clearly presented. The writings are easy to follow, and most of methods and results are described clearly. \n2. Benchmarking results from multiple brain regions, animals, and different versions of Neuropixel probes were conducted, with a significant improvement on decoding accuracy compared to baseline methods. \n3. Given that high-density recordings are getting more widely-used in the neuroscience field, this paper brings impacts with a better decoding tool for studies conducted with this kind of high-density probes.", "weaknesses": "1. While I recognize the paper providing extensive analysis results, the novelty of proposed density-based decoding on unsorted spike waveforms is limited from my perspective. First, the proposed algorithm on using MoG and GLM for decoding has been widely used on neural ensembles, or other types of neural signals like LFP or ECoG. For example, the [paper](https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3972894/) also leveraged MoG for unsorted spike decoding, despite that it did leverage Bayesian decoding framework. This [paper](https://ieeexplore.ieee.org/stamp/stamp.jsp?tp=&arnumber=9585446) proposed Gaussian mixture of model (GMM)-assisted PLS (GMMPLS) for decoding in BMI. Second, there are also other methods proposed for decoding behavior from cluster-less neural data. For example this [paper](https://www.biorxiv.org/content/10.1101/760470v1.full.pdf) proposed cluster-less HMM model for decoding. Another [paper](https://www.biorxiv.org/content/biorxiv/early/2021/08/28/2021.08.26.457795.full.pdf) proposed Gaussian-process multi-class regression decoding on neural population data. The concept of modeling spike distributions with MoG has been applied, and GLM decoding also has been commonly used for neural ensemble. I recognized that some of these papers are not specifically targeting unsorted spikes, and acknowledge this paper bring good impacts for the neuroscience community, but I am not sure if NeurIPS is the proper target given novelty of its proposed algorithm. To meet the bar for NeurIPS, I would like to see more originality algorithm improvements.\n2. Unclear whether baseline decoders are fair baselines. First, Kilsort depends on hyperparameters that need to tune properly for each dataset. From some of example decoded curves, I'm not sure if those baseline decoders were calibrated properly to data. In practice, we rarely worked with decoders with poor accuracy. This leads me to suspect the baseline decoders were not trained/tuned properly for comparison. Second, there are more SOTA decoders that have been used, including previous clusterless decoders on unsorted spiking data (see [paper](https://journals.physiology.org/doi/full/10.1152/jn.01046.2012)). Without a benchmark comparison to these decoders, it's hard to be convinced that the proposed algorithm is much superior than those other decoders. \n3. Unclear if the proposed method only worked better for high-density recordings from Neuropixel or not. Despite that the authors provided two different probe versions for different geometries, it's unclear whether this brings similar benefits for other types of neural recordings (e.g., multiple tetrode arrays, etc). Especially if the paper claims that this method would be better than previous clusterless decoders tested on tetrode recordings, it's important for the NeurIPS audience to know if the proposed method works more generally on decoding unsorted spikes, or somehow better tuned for the distribution from Neuropixel recordings. ", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank reviewer Fmkq for their thorough and thoughtful review. We appreciate their valuable feedback and would like to address their main concerns.\n\n**Weakness 1:**\nPlease refer to the global response for the novelty of our method. \n\n**Weakness 2:**\nIt is important to clarify that the baseline decoders are well-tuned and calibrated properly for each dataset, and we used spike-sorted data tuned by IBL, as described in IBL's spike sorting white paper [(IBL et al ‘22)](https://figshare.com/articles/online_resource/Spike_sorting_pipeline_for_the_International_Brain_Laboratory/19705522). The relatively lower decoding accuracy is due to the inherent difficulty of the decoding task with the behavior data from IBL, which are challenging to decode, especially considering that some brain regions contain limited information about the behaviors. In fact, in the IBL's decoding paper titled \"*A Brain-Wide Map of Neural Activity during Complex Behaviour* [(IBL et al ‘23)](https://www.biorxiv.org/content/10.1101/2023.07.04.547681v2),\" their achieved decoding accuracies are comparable to ours. For detailed numerical values, please consult their Figure S6 and S15 in the supplementary materials. Additionally, in the global response Table 1, we showcase high decoding accuracies for multiple tetrodes data, which involves a simpler decoding task compared to the IBL's behavior tasks.\n\nWe acknowledge the importance of conducting a benchmark comparison of clusterless decoders and appreciate the reviewer for bringing up previous clusterless decoding papers. It is worth noting that the Gaussian-process multi-class regression decoding paper [(Greenidge et al ‘21)](https://www.biorxiv.org/content/biorxiv/early/2021/08/28/2021.08.26.457795.full.pdf) and the GMM-PLS paper [(Foodeh et al ‘21)](https://ieeexplore.ieee.org/stamp/stamp.jsp?tp=&arnumber=9585446) are not specifically designed for clusterless decoding, hence we did not consider comparing our method to those. Regarding a comparison with the previous clusterless decoders, please refer to the global response for details.\n\n**Weakness 3:**\nWe acknowledge the reviewer’s point about adapting our decoder to more general probe geometries. The global response Table 1 demonstrates the effectiveness of our decoder across probe geometries beyond Neuropixel probes. Specifically, on both multiple tetrodes and HD probe data, our decoder outperforms the clusterless point process decoder [(Denovellis et al ‘21)](https://elifesciences.org/articles/64505). It is important to note that our method is not exclusively tailored for Neuropixel recordings, and can be applied more generally to probes with different geometries.\n\n**Question 1:**\nThe key point to emphasize is the high data volume associated with HD probes, which generates a significantly larger amount of data compared to tetrodes. As a result, previous clusterless decoders based on sampling-based inference tend to be slow. This computational cost is evident in the clusterless HMM [(Ackermann et al ‘19)](https://www.biorxiv.org/content/10.1101/760470v1.full.pdf) method, as mentioned by the authors in their paper. They admitted the challenges related to high computational cost, with the model taking about 6 hours to fit a relatively small dataset. Although the clusterless point process model [(Denovellis et al ‘21)](https://elifesciences.org/articles/64505) has a lower computational cost, our decoder is faster and outperforms it in high-density settings, as indicated in global response Table 1 and Figure 2. The limitation of their approach lies in its inability to effectively handle the highly overlapping signals of HD probes, where multiple electrodes are closely packed together. Moreover, our density-based decoder is more flexible by avoiding making restrictive assumption about the underlying system dynamics.\n\nThe original method could potentially use features of other formats, but it's important to note that the localization features, which are HD-specific, play a crucial role. HD probes offer exceptional spatial resolution, and utilizing spike localization features is essential for successful decoding. Relying solely on waveform features on each electrode becomes challenging in the context of HD probes without leveraging the spatial advantages offered by localization features. Our ablation study in Section 8 of the supplementary materials supports this statement, showing that the inclusion of waveform features alongside localization features does not contribute significantly to the decoding tasks when compared to decoding solely based on localization features and peak-to-peak amplitudes.\n\n\n**Question 2:**\nWe use ridge regression for decoding continuous behaviors and L2-penalized logistic regression for decoding binary behaviors."}, {"review_id": "sgkbeG7hXl", "reviewer": "Reviewer_C56U", "summary": "This paper presents a way to decode behavior from neural recordings, bypassing an explicitly spike sorting step. They model individual spikes as coming form a mixture of gaussian distribution, and the assignment probability to different mixtures is used instead of a 'hard' assignment of each spike to a cell. Robust performance of this approach is shown across a large dataset of Neuropixels recordings. ", "questions": "1. Could the paper describe the method as \"soft spike-sorting\"? This would help place the method in the right context. \n\n2. Spike sorting is necessary for biological interpretation of neural recordings. Does soft-spike sorting lend such interpretation? \n\n3. If one doesn't do spike sorting, one could use local field potentials, or other bandpass signals to supplement threshold crossings. One should compare with all those features to claim the superiority of the method. \n\n4. How crucial are the preprocessing steps? Were they developed on the same dataset to optimize spike sorting? If so then the results might be confounded. \n\n5. What are the computational costs of this method compared to spike sorting? \n\n6. Using the KS spike assignments, how 'pure' are the gaussians in the MOG? Are single cells split into different gaussians , or do single gaussians correspond to multiple cells? This could shed light on how important is it to deviate from \"hard\" spike sorting. \n\n7. While spike sorting gives stable decoding across time (inspite of template changes over time), does the current method present with similar stability? Can you learn a decoder from the beginning of a recording and apply it towards the end? \n", "limitations": "Limitations have been adequately addressed.\n\nQuestion for Area chair: The paper identified the source of the data as the data from International Brain Lab. It is not clear if the data is open-source. If not, then it clearly violates the double-blind policy of NeurIPS. Please check. ", "rating": 8, "confidence": 4, "soundness": 4, "presentation": 4, "contribution": 4, "strengths": "While spike sorting is thought to incorporate the biological priors (ie., the fact that spikes come from cells, which carry information in a neural circuit), and promote robust decoding, deviating from the prior by developing a \"soft\" spike sorting approach seems to have helped in this paper. This is surprising and very impactful.\n\nWell written paper. ", "weaknesses": "All the results use very high-density neuropixel recordings. However, it is not clear if the method is a general replacement to spike sorting across different arrays, etc. \n\nThe encoding model depends only on behavior. However, there is a large amount of behavior-independent variability (ex. attention states) that are not incorporated in the encoding model. ", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank Review C56U for the thoughtful review and useful comments. We are glad that the reviewer felt that the manuscript was clearly written and appreciate that the reviewer recognizes that our work makes an important contribution to the literature.\n\n**Weakness 1:**\nAlthough our decoder is mainly designed for high-density (HD) probes, we acknowledge the reviewer’s point about adapting our decoder to more general probe geometries. The global response Table 1 demonstrates the effectiveness of our decoder across probe geometries beyond Neuropixel probes. Specifically, on both multiple electrodes and HD probes data, our decoder outperforms the previous clusterless decoder. \n\n**Weakness 2:**\nWhile our model is primarily conditioned on the behavior of interest, we have the flexibility to turn off this conditioning and utilize a vanilla GMM model to capture behavior-independent variability. This approach comes with a slight trade-off in decoding accuracy. For a detailed comparison between our models with and without dependencies on behavior correlates, please refer to Section 8 (ablation study) in the supplementary materials. We agree that including additional latent variables in the model (e.g. attention states) would be an exciting direction for future work.\n\n**Questions:**\n1. Indeed, our method uses a GMM to estimate the probability of a spike belonging to a mixture component, thereby quantifying the uncertainty in spike assignment. This approach allows the model to make errors and retain valuable information, as opposed to \"hard\" spike sorting methods that discard such uncertainty.\n2. While our method may not offer a single-unit interpretation similar to spike sorting, it brings its own advantages. Traditional spike sorting can be stringent and might discard valuable information by eliminating many cells that do not pass the quality control criteria. In contrast, our approach leverages all unsorted spikes within each brain region for decoding behavioral correlates. As a result, our method retains all available information within that region, providing abundant interpretations for specific brain regions without discarding any valuable data.\n3. We agree with the reviewer's suggestion that the inclusion of LFP and other band-passed signals would be very informative. However, incorporating these data modalities would require extensive efforts beyond the scope of our current work. Our focus in this study is to understand what spikes can inform us about the behavior of interest, rather than optimizing decoding accuracies using all available data modalities. Nonetheless, we acknowledge the potential benefits of incorporating additional data modalities and consider it interesting future work.\n4. No, the preprocessing steps were not developed to optimize spike sorting; they are distinct and separate processes. For instance, destriping is specifically employed for data from IBL due to prominent stripes in the recording data, and it is a preprocessing step developed by IBL. We rely on standard preprocessing from IBL to address quality concerns; see details in their white paper [(IBL et al ‘22)](https://figshare.com/articles/online_resource/Spike_sorting_pipeline_for_the_International_Brain_Laboratory/19705522). Additionally, registration, which removes probe motion, is crucial for data quality, and we utilized a separate method proposed in this paper [(Windolf et al ‘22)](https://www.biorxiv.org/content/10.1101/2022.12.04.519043v1). Spike localization is another essential step in our pipeline, and it is introduced separately in this work [(Boussard et al ‘21)](https://openreview.net/pdf?id=ohfi44BZPC4) as well.\n5. In the global response, we provide a quantification of the computational cost for our method. Regarding Kilosort, we didn't run it ourselves and profile its computational cost, as the spike-sorted output is readily available from IBL's public database. It is worth noting that Kilosort is a real-time sorting algorithm.\n6. Please refer to the section about Figure 3 in the global response for details.\n7. While we can train a decoder at the start of a recording and apply it to the end, there could be a slight impact on decoding accuracy. We use spike features with registered locations that remain constant over time to mitigate probe drift or motion artifacts. The GMM components also remain unchanged over time. However, the mixing proportions of the GMM model vary over time. We conducted experiments by training the model on initial recording segments and then decoding subsequent segments for continuous behaviors on three datasets. The table below shows the mean correlation between true and decoded behaviors. In the context of this experiment, *Time shuffled* indicates training the GMM on segments from various time points within the recordings, while *Time ordered* indicates training on earlier segments to decode later ones.\n\n|              | Dataset 1 | Dataset 2 | Dataset 3 |\n|----------|----------|----------|----------|\n| Time shuffled | 0.760 (0.016) | 0.703 (0.015) | 0.802 (0.012) |\n| Time ordered | 0.733 (0.067) | 0.687 (0.008) | 0.792 (0.025) |\n\n**Question for Area Chair:**\nWe use open-source IBL datasets in our paper. The public available datasets can be accessed [here](https://int-brain-lab.github.io/iblenv/public_docs/public_introduction.html). "}, {"review_id": "ylrEf8dy42", "reviewer": "Reviewer_Cpfb", "summary": "The authors demonstrate a novel approach to decoding behavioral variables from neural activity recording using dense silicon probes. Most approaches utilize a sort-then-decode line of attack – first channels are spike sorted, then decoders are trained on unambiguously identified single-units. Here, the authors present a probabilistic model that bypasses explicitly assigning waveforms to units; rather, the authors model spike features using a mixture of Gaussians, which is then leveraged for decoding. The authors go on to show that their method generally outperforms the current state-of-the-art in a variety of experiments spanning recordings in the mouse and in the monkey.", "questions": "1. Did the authors check to see if their technique picks up on noise or motion artifact in these recordings? ", "limitations": "Limitations are adequately addressed.", "rating": 8, "confidence": 3, "soundness": 4, "presentation": 3, "contribution": 3, "strengths": "Strengths:\n\n1. The authors present an elegant solution a highly-challenging problem in neuroscience – decoding behavior in the presence of noisy and highly-overlapping signals. \n1. The results are generally compelling. Clearly, the author's algorithm outperforms Kilosort and multi-unit thresholding in a variety of decoding contexts.\n1. The core idea and implementation is fairly straightforward. Given this, I expect this contribution to be generally useful to the community.", "weaknesses": "Weaknesses:\n\nMajor:\n\n1. One potential issue with \"bypassing\" spike-sorting is noise and motion artifact. Spike-sorting, if carefully curated, can remove unwanted sources of contamination in extracellular electrophysiology. Here, I am somewhat concerned that this method could pick up and amplify contaminating signals that may contain behaviorally relevant information. For instance, Figure 3a shows that \"good\" KS units in CA1 and DG contain barely any information about motion energy; however, multi-unit thresholding, all KS, and the author's method find a surprising amount of information. Do the authors think this is reasonable? I am a bit worried since similar gains are not seen for choice decoding, which may not be as easy to decode from motion artifact. This could absolutely be real information uncovered by the model, but this issue is worth addressing head on. Do the authors think that their method could be picking up on any noise or motion artifact in the original recordings?\n\nMinor:\n\n1. Rhetorical point. Line 47 the authors mention prior attempts to \"bypass\" spike sorting as working with \"limited\" features such as maximum amplitudes and principal components of the spike waveform. Unless I missed something this language seems overly critical since the authors are only using peak-to-peak amplitudes and spike localization features (x, z). \n1. Related to the last point, did the authors consider waveform principal components or aspects of the waveforms other than amplitude? Spike width, e.g., is known to contain information about cell type.\n1. Line 97 typo", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank reviewer Cpfb for their thorough and thoughtful review. We appreciate that the reviewer recognized the contribution of our method to the neuroscience community.\n\n**Major Weakness:**\nWe performed motion correction (registration) during the preprocessing step, ensuring that motion artifacts are minimized in our data. Figure 5 in the supplementary materials provides a comparison of spike raster plots with and without registration, confirming the minimal impact of motion artifacts in the preprocessed data. Additionally, Kilosort has its own motion correction step, reducing the influence of motion artifacts on the spike-sorted decoder.\n\nOur decoder exhibits strong performance not only in decoding motion energy but also in decoding wheel speed and other continuous behaviors that are unrelated to motion artifacts. Decoding continuous behavioral variables is generally easier to improve compared to decoding binary choices. We acknowledge the importance of improving choice decoding, and we are actively exploring this direction in our ongoing research.\n\nThe subpar performance of \"good KS\" is due to IBL's stringent quality control criteria for spike sorting. Only those single-units that pass the quality control procedure are categorized as \"good\" KS. In some brain regions like CA1 and DG, there are often less than 10 units that meet these criteria, explaining the poor decoding performance of \"good\" KS in these areas. This is precisely why decoding using all Kilosort units (including probable multi-unit activity) yields significantly better results.\n\n**Minor Weakness:**\n1. We apologize for being overly critical with our choice of words. Our intention was to emphasize the importance of spike localization features, which are uniquely applicable to high-density probes with high spatial resolution. This novel technique, as described in this paper [(Boussard et al ‘21)](https://openreview.net/pdf?id=ohfi44BZPC4), is of great importance as it extracts valuable information that cannot be achieved solely through waveform features. Our ablation study in Section 8 of the supplementary materials supports this statement, showing that the inclusion of waveform features alongside localization features does not contribute significantly to the decoding tasks when compared to decoding solely based on localization features and peak-to-peak amplitudes.\n2. We appreciate the insightful suggestion from the reviewer. We did explore the possibility of incorporating more waveform features to enhance decoding performance. As shown in our ablation study in Section 8 of the supplementary materials, we included waveform principal components in addition to localization features and peak-to-peak amplitudes. However, supplementary Table 3 clearly demonstrates that the addition of waveform features did not result in a significant improvement in the decoding tasks compared to using only localization features and peak-to-peak amplitudes.\n3. We thank the reviewer for bringing this typo to our attention, and we will correct it in the future version.\n\n**Question 1:**\nWe examined our method to ensure that it did not capture motion artifacts in the recordings. We manually inspected the spike raster plots after registration, and we found no evidence of motion artifacts affecting the decoding results. For a more comprehensive discussion on this topic, please refer to our response in the Major Weakness section. "}]

Author rebuttal:
We appreciate the reviewer's feedback on our manuscript. We were encouraged that the reviewers recognized our paper's empirical robustness and its solution to the challenging task of decoding behaviors from highly overlapping neural signals. To address the concerns, we conducted the following experiments:
- Benchmark comparison to other clusterless decoders.
- Simulation for model validation.
- Computational time cost comparison.
- Correspondence between spike sorting and GMM assignment. 

One important point to emphasize is the novelty of our method, which involves conditioning the generative model for the spike feature distribution on the behavior of interest, making the inference problem non-trivial. See supplementary materials for detailed derivations. **The probabilistic model and inference method presented have broader relevance within the NeurIPS community. The conditioning of the data generating process on external variables allows for improved prediction in various scenarios.  Our approach is not limited to GMM and GLM; instead, it can accommodate non-Gaussian mixture models and nonlinear models.** While our paper serves as a proof-of-concept, the proposed method can be extended to tackle diverse modeling challenges. 

**Our novelty also lies in the first demonstration of a clusterless decoding method designed for high-density (HD) probes (which are now used in hundreds of labs), and our use of localization features which are only meaningful for HD probes.** By combining these factors, we overcome previous challenges and improve decoding performance. Furthermore, compared to previous decoders based on state-space models, our density-based decoder offers increased flexibility and scalability. Our method avoids explicit assumptions about underlying system dynamics, making it more flexible at capturing complex relationships in the data. The increased scalability is critical for large data volumes generated by HD probes.

We appreciate the suggestion to benchmark clusterless decoders beyond HD probes. However, limited available code in clusterless decoding papers makes comparisons challenging within a short timeframe. While the suggested Bayesian decoding paper [(Chen et al ‘12)](https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3972894/) lacks functionality for the main task (due to incomplete code), we re-implemented the clusterless HMM model [(Ackermann et al ‘19)](https://www.biorxiv.org/content/10.1101/760470v1.full.pdf), but the training time was prohibitive. Despite these limitations, we made efforts to compare our method to clusterless point process decoders [(Denovellis et al ‘21)](https://elifesciences.org/articles/64505) on both tetrodes and HD probes. This clusterless point process decoder is similar to the Bayesian decoding paper, as both use a marked point process to link spike features and behaviors. We used the same spike features and calibrated each decoder to the data. For multiple-tetrode data, both decoders used spike amplitudes from 4 channels of 5 tetrodes. Similarly, for HD probes, both decoders relied on spike localization features and peak-to-peak amplitudes. **Our method's advantage over the clusterless point process decoder is evident in Table 1,  owing to the increased flexibility of our decoder compared to earlier clusterless state-space models.**

We conducted simulations to illustrate the principles of our method. The simulation aimed to show that our encoding model can learn the relationship between spike features and behaviors. We performed two tasks, decoding a binary variable, $y_k$, simulated from a Bernoulli distribution, and decoding a continuous variable, $y_k$, simulated from a Gaussian process. To mimic the data-generating process, we selected Gaussian components with "templates" extracted from a real dataset. The encoding model parameters, $b$ and $\beta$, were also taken from learned parameters in the same dataset. Given $b, \beta$, and $y_k$, we simulated the "firing rates" $\lambda$ for each Gaussian component in the mixture, as described in the Method section of our paper. Next, we generated spike features based on these simulated ''firing rates'', and applied the encoding model to infer the behavior-dependent $\lambda$. Figure 1 displays the learned $\lambda$ for each component $c$, time $t$, and trial $k$. **The learned "firing rates" $\lambda$ closely resembled the simulated ones, indicating the model's ability to recover the primary associations between spike features and behaviors. With such associations, the decoding model can decode behaviors.**

We appreciate the question about the computational cost of our method. For spike sorting, we did not personally run Kilosort (KS), given the accessible spike-sorted output from IBL's public database. Notably, Kilosort operates in close to real-time, implying that it takes 1000 seconds to sort a 1000-second recording. In Figure 2, we provided a computational time comparison relative to real-time. **Our decoding step operates at a sub-real-time pace (0.3 times real-time), which is 4 times faster than the point process decoder (1.2 times real-time). The total time after preprocessing for our method is close to real-time.**  While we didn't measure the time cost of the clusterless HMM due to the extensive model fitting time, the authors acknowledged that their model takes around 6 hours to fit on a small dataset in the paper.

We appreciate the inquiry about the insights our method provides into "hard" spike sorting. In response, we computed an agreement matrix (Figure 3) between "hard" KS assignments and "soft" GMM assignments. We calculated the conditional probability of spikes belonging to each GMM component, given that these spikes belong to the corresponding KS unit. Notably, KS units with large amplitudes are less likely to be split into multiple Gaussian components. In conclusion, **Figure 3 shows a reasonable correspondence between the Gaussian components and the spike-sorted units.** 

Meta-review: The reviewers largely agree that this is an interesting and useful piece of work. Their opinions have been strengthened by the author responses, and the additional results/clarification provided. Please be sure to include these in the final version of the manuscript.
Decision: Accept (spotlight)