Bubblewrap: Online tiling and real-time flow prediction on neural manifolds
Authors: Anne Draelos
Pranjal Gupta
Na Young Jun
Chaichontat Sriworarat
John Pearson
Venue: NeurIPS 2021 Poster (2021)
Keywords: neuroscience
neural populations
online
real-time

Abstract: While most classic studies of function in experimental neuroscience have focused on the coding properties of individual neurons, recent developments in recording technologies have resulted in an increasing emphasis on the dynamics of neural populations. This has given rise to a wide variety of models for analyzing population activity in relation to experimental variables, but direct testing of many neural population hypotheses requires intervening in the system based on current neural state, necessitating models capable of inferring neural state online. Existing approaches, primarily based on dynamical systems, require strong parametric assumptions that are easily violated in the noise-dominated regime and do not scale well to the thousands of data channels in modern experiments. To address this problem, we propose a method that combines fast, stable dimensionality reduction with a soft tiling of the resulting neural manifold, allowing dynamics to be approximated as a probability flow between tiles. This method can be fit efficiently using online expectation maximization, scales to tens of thousands of tiles, and outperforms existing methods when dynamics are noise-dominated or feature multi-modal transition probabilities. The resulting model can be trained at kiloHertz data rates, produces accurate approximations of neural dynamics within minutes, and generates predictions on submillisecond time scales. It retains predictive performance throughout many time steps into the future and is fast enough to serve as a component of closed-loop causal experiments.

Bubblewrap: Online tiling and real-time flow prediction on neural manifolds
Anne Draelos Biostatistics & Bioinformatics
Pranjal Gupta
Psychology & Neuroscience
Duke University anne ddraelos@duke edu
Duke University
pranjal. guptaduuke.ee..eee. edu
Na Young Jun Neurobiology Duke University nayoung. jundduke.eduee.o. edu
Chaichontat Sriworarat
Biomedical Engineering
Duke University chaichontat. sdduke. edu
John Pearson
Biostatistics & Bioinformaties
Electrical & Computer Engineering
Neurobiology
Psychology & Neuroscience
Duke University
john. pearsooarood eduu
Abstract
While most classic studies of function in experimental neuroscience have focused the coding properties of individual neurons, recent developments in recording technologics have resulted in an increasing emphasis on the dynamics of neural populations. This has given rise to wide variety of models for analyzing pop- ulation activity in relation to experimental variables, but direct testing of many neural population hypotheses requires intervening in the system based current neural state, necessitating models capable of inferring neural state online. Exist- ing approaches, primarily based dynamical systems, require strong parametric assumptions that are easily violated in the noise-dominated regime and do not scale well to the thousands of data channels in modern experiments. To address this problem, we propose a method that combines fast, stable dimensionality re- duction with soft tiling of the resulting neural manifold, allowing dynamics to be approximated as a probability flow between tiles. This method can be fit efficiently using online expectation maximization, scales tens of thousands tiles, and outperforms existing methods when dynamics are noise-dominated feature multi-modal transition probabilities. The resulting model can trained kiloHertz data rates, produces accurate approximations of neural dynamics within minutes, and generates predictions on submillisecond time scales. retains predic- tive performance throughout many time steps into the future and is fast enough serve as a component of closed-loop causal experiments.
Introduction
Systems neuroscience is in the midst of a data explosion. Advances in microscopy and probe technology 4 have made it possible to record thousands of neurons simultaneously in behaving
35th Conference on Neural Information Processing Systems (NcurIPS 2021).

---

animals. At the same time, growing interest in naturalistic behaviors has increased both the volume and complexity of jointly recorded behavioral data, On the neural side, this has resulted in a host of new modeling and analysis approaches that aim to match the complexity of these data, typically using artificial neural network models as proxies for biological neural computation 6
At the same time, this increase in data volume has resulted in increasing emphasis on methods for dimensionality reduction p and focus on neural populations in preference to the coding properties of individual neurons [0 However, given the complexity of neural dynamics, it remains difficult to anticipate what experimental conditions will be needed to test population hypotheses in post hoc analyses, complicating experimental design and reducing power. Conversely, adaptive experiments, those in which the conditions tested change in response to incoming data, have been used in neuroscience to optimize stimuli for experimental testing [T] 22 [ [4] in closed- -loo
designs 15 16. 17,, and even to scale up holographic photostimulation for inferring functional connectivity in large circuits [8]
Yet, despite their promise, adaptive methods are rarely applied in practice for two reasons: First, although efficient online methods for dimensionality reduction exist 19 PC 21 22 [3
these methods do not typically identify stable dimensions to allow low-dimensional representations of data to be compared across time points. That is, when the spectral properties of the data are changing in time, methods like incremental SVD may be projecting the data into an unstable basis, rendering these projections unsuitable as inputs to further modeling. Second, while many predictive models based on the dynamical systems approach exist
29 including online approaches 30 311 16 R2), they typically assume system with lawful dynamics perturbed by Gaussian noise. However, many neural systems of interest are noise-dominated, with multimodal transition kernels between states.
In this work, we are specifically interested in closed loop experiments in which predictions of future neural state are needed in order to time and trigger interventions like optogenetic stimulation or a change in visual stimulus. Thus, our focus is on predictive accuracy, preferably far enough into the future to compensate for feedback latencies. To address these goals, we propose an alternative to the linear systems approach that combines a fast, stable, online dimensionality reduction with a semiparametric tiling of the low-dimensional neural manifold. This tiling introduces a discretization of the neural state space that allows dynamics to be modeled as a Hidden Markov Model defined by a sparse transition graph. The entire model, which we call "Bubblewrap." can be learned online using a simple EM algorithm and handles tilings and graphs of up to thousands of nodes at kiloHertz data rates. Most importantly, this model outperforms methods based on dynamical systems in high-noise regimes when the dynamics are more diffusion-like. Training can be performed at a low, fixed latency x 10ms using GPU, while cached copy of the model in main memory is capable of predicting upcoming states at Ims latency. As a result, Bubblewrap offers a method performant and flexible enough to serve as neural prediction engine for causal feedback experiments.
2 Stable subspaces from streaming SVD
As detailed above, one of the most pressing issues in online neural modeling is dealing with the increasingly large dimensionality of collected data - hundreds of channels per Neuropixels probe
4 15, tens of thousands of pixels for calcium imaging. However, as theoretical work has shown
36 k true neural dynamics often lie on a low-dimensional manifold, so that population activity
can be accurately captured by analyzing only a few variables.
Here, we combine two approaches to data reduction: In the first stage, we use sparse random projections to reduce dimensionality from an initial d dimensions (thousands) to n (a few hundred) By simple scaling, for a fixed budget of N cells in our manifold tiling, we expect density (and thus predictive accuracy) to scale as N à in dimension n. and so we desire n to be as small as possible. However, by the Johnson-Lindenseee Lindenstrauss Lemma [37 36), when reducing from n dimensions, the distance between vectors u, and v, in the reduced space is related to the distance between their original versions u and v by
B55
36
(1) with probability 1 - 5 if > 0log(1/6)/e3) Unfortunately, even for & ~ 0.1 (10% relative error), the required n may be quite large, making this inappropriate for reducing to the very small numbers of effective dimensions characterizing neural datasets.
s (1 e) |Iu - -0|2 < ||u. 1. = (1+ x||
2

---

Algorithm 1 Procrustean SVD (proSVD)
Given: Initial data Xo. decay parameter (0, 1]
2: Initialize: QR Factorization: Xo Qo
3:
4: for = 1...... do
5:
Fetch new columns of data, X+
6:
C X+. X1 < X+ Q--1C,
Q1, R1 + QR({1))
Gram-Schmidt
Q - Q1). R+ R1-1 C
7:
QR of augmented data
R
8:
- SVD(R))
9:
I+
D Discount old data > U1 contains the first & columns of U Orthogonal Procrustes: minr||QU,T Q--ill
10:
M M4-Q-IQU = []x 0kxt] U1 T-UVT
11:
12:
Qi QUTT
Update left subspace basis > QR right subspace, update inner block
13:
R, + QR(V).
R + TE2Q0
14:
end for
Thus, in the second stage, we reduce from n ~ ((100) to k O(10) dimensions using streaming singular value decomposition. This method based on the incremental block update method of
120 22 with an important difference: While the block update method aims to return the top-k SVD every time point, the directions of the singular vectors may be quite variable during the course of an experiment (Figure Id) which implies unstable representation of the neural manifold. However, as we show below, the top-k subspace spanned by these vectors stabilizes in seconds on typical neural datasets and remains so throughout the experiment. Therefore, by selecting stable basis (instead of the singular vector basis) for the top-k subspace, we preserve the same information while ensuring a stable representation of the data for subsequent model fitting.
More specifically, let x, I be vector of input data after random projections. In our streaming setup, these are processed samples at time, with b = reasonable for slower methods like calcium imaging and 40 more appropriate for electrophysiological sampling rates of ~20kHz. Then, if the data matrix X has dimension x T, adding columns over time, the incremental method of produces at each time step a factorization X Q QRW where the columns of the orthogonal matrices Q and W span the left and right top-k singular subspaces, respectively. If the matrix R were diagonal, this would equivalent to the SVD. In the incremental algorithm, R is augmented each timestep based on new data to form Rt which block diagonalized via an orthogonal matrix and truncated the top-k subspace, allowing for exact reduced-ram SVD (Appendix A).
20 22
However, as reviewed in subspace, there are likewise multiple choices of block diagonalization for R. In (0 22], the form of this operation is chosen for computational efficiency. But an equally valid option is to select the orthogonal matrix that minimizes the change in the singular subspace basis Q from one timestep to the next:
since there are multiple choices of basis Q for for the top-k singular
min | Q.- || quinl|QU,TT' - Q1-1
where Q is augmented basis for the top-(k: + b) singular subspace, U, contains the first k left singular vectors of A. and 7 is an orthogonal matrix (Appendix A This minimization is known as Orthogonal Procrustes problem and has a well-known solution 38 T ÜV T, where Ü and V are the left and right singular vectors, respectively, of M = QUIQQU (See 39 for a recent application of similar ideas in brain-computer interfaces). This Procrustean SVD (proSVD) procedure is summarized in Algorithm ! There, lines 1-8 follow 20, 22), while lines 10 and 11 perform the Orthogonal Procrustes procedure. Line 9 serves as a leak term that discounts past data in 40
Figures - illustrates the performance the two-stage dimension reduction algorithm for case of d = I randomly generated Gaussian data. While ProSVD yields minimal distortion (due to truncation of the spectrum to A: = 6), random projections require & ~ O(100) to achieve the same result (Figure Da). By contrast, random projections are much faster (Figure Thus, we can trade off distortion against time by adjusting n, the number of intermediate dimensions. As Figure T: shows, the optimal tradcoff occurs somewhere around 200 for this example.
3

---

distortion of reducing N 1e 04 dims to dims
time to reduce
pareto front of distortion timing
a
N 1e+04 dims to k dims
0.6
SVD to
timee (mr)
10
distortion
RP to
0.10
distertion
0.4
RP
5V0 n
SVD M to
107
0.05
RP (1
0.2
0.0
10c
2000
4000
6000
8000
10000
2000
4000
6000
8000
10000
0.00
0.05
0.10
0.15
time
Projection single trial onto singular vectors (top) or proSVD basis vectors (bottom)
10 bins (0.15 s) seen
100 bins (1.5 s) seen
1000 bins (15.0 s) seen
streaming SVD
singular vector
singulas vector 2
singular vecter
whole data SVD
singular vector
singular vector
singular
Grssod basee e Sectorr
proSvD
prosvo
proSVD
prosVD
basis vector
basis
prosvD basis vector
Figure 1: Timing and stability of two-stage dimension reduction. a) Distortion (e) as a function of number of dimensions retained (n) for both sparse random projections and proSVD on random Gaussian data with batch size b = 1000. b) Time required for the dimensionality reduction in (a), amortized for batch size. While random projections are extremely efficient, proSVD time costs grow with the number of dimensions retained. c) Pareto front for the time-distortion tradeoff of random projections followed by pro Color indicates n, the number of dimensions retained by random projections. Black arrow indicates the particular tradeoff we chose of n 200. d-f) Embedding of single trial (green line) into the basis defined by streaming SVD for different amounts of data seen. Dotted line indicates the same trial embedded using SVD on the full data set. Rapid changes
estimates of singular vectors early on lead to an unstable representation. g-i) Same trial and conventions as (d-f) for the proSVD embedding. Dotted lines in the two rows are the same curve in different projections.
Figures - show results for neural data from recorded from monkey motor cortex 2 in cued reach task. While projection the data into the basis defined by streaming SVD remains unstable early in data collection (top), the proSVD representation is nearly equivalent to the full offline result after only a few trials 155,,, middle). This due to the fact that, all data sets we examined, the top-k SVD subypace identified extremely quickly; proSVD simply ensures the choice of stable basis for that subspace.
3 Bubblewrap: a soft manifold tiling for online modeling
As reviewed above, most neural population modeling approaches are based the dynamical systems framework, assuming lawful equation of motion corrupted by noise. However, for animals engaged task-free natural behavior 41 42 4x), trajectories are likely to sufficiently complex that simple dynamical models fail. For instance, dynamical systems models with Gaussian noise necessarily pro- duce unimodal transition probabilities centered around the mean prediction, while neural trajectories may exhibit multimodal distributions beginning at the same system state. By contrast, pursue an

---

alternative method that trades some accuracy in estimating instantaneous system state for flexibility modeling the manifold describing neural activity.
Our approach is to produce a soft tiling of the neural manifold in the form of Gaussian mixture model (GMM), each component of which corresponds between tiles via a Hidden Markov Model (HMM), which allows us to capture multimodal probability flows. As the number of tiles increases, the model produces increasingly finer-grained description of dynamics that assumes neither an underlying dynamical system nor particular distribution of noise.
single tile. We then approximate the transitions
More specifically, let zy be the low-dimensional system state and let 2 E 1...NN index the tile to which the system is assigned time t. Then we have for the dynamics
p(z j|z011 = Asj p(a))2)) .V (M,,, Ex,) plu,, NIW >g. Vy, (3)
where we have assumed Normal- cetmmal-nvveree--ee -iishart priors the parameters of the Gaussians. Given its exponential family form and the conjugacy of the priors, online expectation maximization updates
available in closed form [4] 45 46 for each datum, though we opt, as in 45 for gradient-based optimization of an estimate the evidence lower bound
C(A, uu, -
(T) + B.g - log As + ((i,,(T) + AsHOO) " Ps
(4)
Σ tr((v, + (T) + AyMotHosj + +A, + (n,(1)))55 11251)
2
j
Σ (vs + (T) + d+ 2) log det Σ
2
with accumulating (estimated) sufficient statistics
ag(t) Σ as(t (t)
(t(() (1 e,) taatt (t)
(5)
ht(t) Σ (t)
(t,(t) = (1 - - 1) + at)
(2) (1 (1- e.) (1- + a (1)
where (1( = p(2; = =j) is the filtered posterior, I (1) is the update matrix from the forward algorithm [4], and E, is forgetting term that discounts previous data. Note that even for = = 0, L is only an estimate of the evidence lower bound because the sufficient statistics are calculated using att and not the posterior over all observed data.
In setting Normal-Inversee. Inverse-Wishart priors over the Gaussian mixture components, take empirical Bayes approach by setting prior means May to the current estimate of the data center of mass and prior covariance parameters V to N- times the current estimate of the data covariance (Appendix BB. For initializing the model we use small data buffer M O(10). We chose effective observation numbers (A, = 10-3 and trained this model to maximize C(A, μu, 2)) using Adam 47,, enforcing parameter constraints by replacing them with unconstrained variables a.j and lower triangular L; with positive diagonal: A.y = exp(a.y)) Σ, exp(a.g), Σ, L,
Finally, in order to prevent the model from becoming stuck in local minima and encourage more effective tilings, we implemented two additional heuristics as part of Bubblewrap: First, whenever a new observation was highly unlikely to be in any existing mixture component (log dog,(z; (24) O,, for all 22), we teleported node at this data point by assigning (t) = 1 for an unused index J. For initial learning this results "breadcrumbing" approach where nodes are placed at the locations of each new observed datum. Second, when the number of active nodes was equal to our total node budget N, we chose to reclaim the node with the lowest value of ii()) and zeroed out its existing sufficient statistics before teleporting it to a new location. In practice, these heuristics substantially improved performance, especially early in training (Appendix D). The full algorithm is summarized in Algorithm 2
5

---

Algorithm 2 Bubblewrap
1: Given: Hyperparameters Xj, Vj, 3t, forgetting rate Et, teleport threshold 8, step size initial
data buffer M
2: Initialize with {xI... 2M M4 nr. Ey E, aij a. d2
3:
4: for / 1...dd do
5:
Observe new data point Zt.
6:
if' log p(z2 [x) < 8 42z then
b Teleport
7:
end J4J 2t, (t) = for J = arg min, (t)
8:
9:
Calculate i; (t) via forward filtering 44
Update sufficient statistics via 5
step
I Σ, Su Σ E. Suz
11:
Global mean and covariance update
9 ~ NV,N 1²) Poj t apoj + (1-a)p cj. V3
Update priors (Appendix
N
Perform gradient-based update of L(A,
M step
14: end for
Experiments
We demonstrated the performance of Bubblewrap both simulated non non-linear dynamical systems and experimental neural data. We compared these results to two existing online learning models for neural data, both of which based dynamical systems systems, we generated noisy trajectories from a two-dimensional Van der Pol oscillator and a three- dimensional Lorenz attractor. For experimental data, used four publicly available datasets from range of applications: 1) trial-based spiking data recorded from primary motor cortex monkeys
3T 32
To simulate low-dimensional
performing reach task 4 40
preprocessed by performing online jPCA 4 2) continuous video data and 3) trial-based wide-field calcium imaging from a rodent decision-making task (50 51: 4)
high-throu Neuropixels data 52 53
For each data set, we gave each model the same data as reduced by random projections and ProSVD. For comparisons across models, quantified overall model performance by taking the mean log predictive probability over the last half each data set (Table Tb. For Bubblewrap, prediction T' steps into the future gives
log (og log E
= =
51 i)p(z; = i(p1.e)
= log
(6)
where AT the 7'-th power of the transition matrix. Conveniently, these forward predictions can be efficiently computed due to the closed form 66. while similar predictions in comparison models 5 must be approximated by sampling (Appendix C In addition, for Bubblewrap. which is focused on coarser transitions between tiles, we also report the entropy of predicted transitions:
BC
H(t,
Σ pleur= f(æ1a) log p(244 f(æ1a)
(A4)) log
Additional detailed experimental results and benchmarking of our GPU implementation in JAX 154
are Appendix [ We compared performance of our algorithm against both [30) (using our own implementation in JAX) and Variational Joint Filtering [32]] (using the authors' implementation) Our implementation of Bubblewrap. as well as code to reproduce our experiments, is open-source and available online at http:://githu.. wttcll..comppeeee com/pearsonlab. paaaarssonnaa/uubbberrrggg /Bbbblewrap..
When tested on low-dimensional dynamical systems, Bubblewrap successfully learned tilings of both neural manifolds, outperforming VJF 332 on both datasets (Figure 2b.b) while it was comparable to algorithm of 30 on one of the 2D (but neither of the 3D) cases (Figure 2). This is surprising,
since both comparison methods assume an underlying dynamical system and attempt to predict

---

Table I: Model comparison results as mean + standard deviation of the log predictive probability
over the last half of the dataset. Asterisks (*) indicate models that degenerated a random walk.
Log predictive probability
Dataset
Bubblewrap
VJF 32
ZP (2016) 30
2D Van der Pol, 0.05
0.965 + 1.123
-..38 + 0.427
0.121 + 0.857
2D Van der Pol, 0.20
-1.088 = 1.184
-1.140 = 0.879
-0.506 = 0.964
3D Lorenz, 0.05
-7.338 + 1.289
-16.98 + 1.9.2
-12.39 + 1.723*
3D Lorenz, 0.20
+
-17.30+22.112
-12.42 + 1.708
Monkey reach Wide-field calcium
3.046 + 4.959
-5.159 + 0.987
3.818 9.118
5.974
3.768 6.204
1.613 + 4.083
Mouse video
-00993 + 2.386
-15.66 = 1.084
-0065 + 4.145*
Neuropixels
12.84 + 6.017
-12.06
5.244
+ 5.244
-12.2 + 4.567
Lerenz
Van der Pol
-1000
200
-3-
-1400
122
2P 125281
-2
0
Une shop
Figure 2: Modeling of low-dimensional dynamical systems. a) Bubblewrap end tiling of a 2D Van der Pol oscillator (data in gray; 5% noise case corresponding to line of Table D Tile center locations are in black with covariance 'bubbles' for 3 sigma in orange. b) Bubblewrap end tiling of a 3D Lorenz attractor (5% noise), where tiles are plotted similarly (a). e) Log predictive probability across all timepoints for each comparative model for the 2D Van der Pol, 0.05 case (top) and for the 3D Lorenz, 0.05 case (bottom).
differences between individual data points, while Bubblewrap only attempts to localize data to within a coarse area of the manifold.
We next tested each algorithm on more complex data collected from neuroscience experiments. These data exhibited a variety of structure, from organized rotations (Figure 3 to rapid transitions between noise clusters (Figure Bb)) to slow dynamics (Figure 3b) In each case, Bubblewrap learned a tiling of the data that allowed to equal or outperform state predictions from the comparison algorithms (Figure 3 f. blue). In some cases, as with the mouse dataset, the algorithm of (30 produced predictions for x, by degenerating to a random walk model [ with *; Appendix D) Regardless, Bubblewrap's 's tiling generated transition predictions with entropies far below those of a random walk (Figure Bd-f, green), indicating successfully identified coarse structure, even in challenging datasets. Thus, even though these data are noise-doming and lack much of the typical structure identified by neural population models, noarse-graining identifies some reliable patterns.
We additionally considered the capability of our algorithm to scale to high-dimensional or high- sampling rate data, As a case study, we considered real-time processing (including random projections, ProSVD, and Bubblewrap learning) of Neuropixels data comprising 2688 units with 74,330 timepoints from 30 ms bins. As Figure 4 shows. Bubblewrap once again learns a tiling of the data manifold (a). capturing structure in the probability flow within the space (b) with predictive performance comparable to finer-grained methods (Table MD. More importantly, all these steps can performed well within the 30ms per sample time of the data (c). In fact, when testing on representative examples of d 10a dimensions, kHz sampling rates, or N 20, 000 tiles, our algorithm was able to maintain amortized par-sample processing times below those of data acquisition. In practice,
7

---

25.0
150
15.0
15.0
12.
-10
-20-
5.0
-30-
3000 35.00
2000
20000
40000 09300
10000
Figure 3: Bubblewrap results on experimental datasets. a) Bubblewrap results for example trials (blue) from the monkey reach dataset 48 49 projected onto the first jPCA plane. All trials are shown in gray. The tile center locations which were closest to the trajectories are plotted along with their covariance "bubbles." Additionally, large transition probabilities from each tile center are plotted as black lines connecting the nodes. Bubblewrap learns both within-trin -trial and across-trre trross--rial transitions, as shown by the probability weights. b) Bubblewrap results widefield calcium imaging from 50 51, visualized with UMAP. A single trajectory comprising 1.5. data shown blue. Covariance "bubbles" and transition probabilities omitted for clarity. c) Bubblewrap results when applied to videos of mouse behavior (50) 51 visualized by projection onto the first. SVD plane. Blue line: 3.3s of data. d, f) Log predictive probability (blue) and entropy (green) over time for the respective datasets in (a,b.c). Black lines are exponential weighted moving averages of the data. Dashed green line indicates maximum entropy (log, (N)).
found that even in higher-dimensional datasets (as in the Neuropixels case), only 1-2 thousand tiles used by the model, making easy round trip GPU latencies to perform gradient updates, online predictions using slightly stale estimates of Bubblewrap parameters could be performed far faster, in tens of microseconds.
kHz data rates. What's more, while learning involved
Just as importantly, when used for closed loop experiments, algorithms must be able to produce predictions far enough into the future for interventions be feasible. Thus examined the performance of our algorithm and comparison models for predicting T steps ahead into the future. Bubblewrap allows to efficiently calculate predictions even many time steps into the future using shows the
(6) whereas the comparison models require much costlier sampling approaches. Figure 5 mean log predictive probabilities for all models many steps into the future for each experimental dataset (top row), and the entropy of the predicted transitions using Bubblewrap (bottow row). Our algorithm consistently maintains performance even when predicting 10 steps ahead, providing crucial lead time to enable interventions at specific points along the learned trajectory. In comparison, predictive performance of (30]. which initially matches or exceeds Bubblewrap for two datasets, rapidly declines, while Variational Joint Filtering 32 with lower log likelihood, also exhibits a slow decay in accuracy.
Discussion
While increasing attention has been paid in neuroscience to population hypotheses of neural function 100 and while many methods for modeling these data offline exist, surprisingly few methods function online, though presumably online methods will be needed to test some population dynamics hypotheses [17] While the neural engincering literature has long used online methods based on

---

proj
bro6M02 Bunblewrap
log peotictie probzbilty
(aaa) Adaau Time (ins)
(hhi
30
10-1
10'
20000
40000
60000
timestaps
NP
NPb
10' dim
20k
1k
Prediction
Figure 4: High-throughput data & benchmarking. a) Bubblewrap results for example trajectories (blue) in the Neuropixels dataset [52 53 (data in gray) visualized with UMAP. b) Log predictive probability (blue) and entropy (green) over time. Black lines are exponential weighted moving averages of the data. Dashed green line indicates maximum entropy. e) Average cycle time (log scale) during learning prediction (last bar) for each timepoint. Neuropixels (NP) is run as in (a.b) with no optimization and all heuristics, and Bubblewrap easily able to learn at rates much faster than acquisition (30 ms). By turning off the global mean and covariance and priors updates and only taking a gradient step for C every 30 timepoints, we are able to run at close to 1 kHz (NPb). All other bars show example timings from Van der Pol synthetic datasets optimized for speed: 10 dim, where we randomly project down to 200 dimensions and used proSVD to project to 10 dimensions for subsequent Bubblewrap modeling learning: N 20k, 10k, and 1 nodes, showing how algorithm scales with the number of tiles; and Prediction, showing the time cost to predict one step ahead for the N 1k case.
Kalman filtering, (e.g.. [6).. and these methods are known to work well many practical cases, they also imply strong assumptions about the evolution of activity within these systems. Thus, many studies that employ less constrained behavior or study neural activity with less robust dynamics may benefit from more flexible models that can be trained while the experiment is running.
Here, address this need, we have introduced both new dimension reduction method that rapidly produces stable estimates of features and a method for rapidly mapping and charting transitions on neural manifolds. Rather than focus on moment-by-moment prediction, we focus on estimating a coarse tiling and probability flow among these tiles. Thus, Bubblewrap may be less accurate than methods based on dynamical systems when state trajectories are accurately described by smooth vector fields with Gaussian noise. Conversely, when noise dominates, is multimodal, or only large- scale probability flow discernible over longer timescales, Bubblewrap is better poised to capture these features. We saw this in our experiments, where the model of 30 exhibited better overall performance in the mouse video dataset (Figure 5 when it did not learn to predict and degenerated to random walk. Indeed, the most relevant comparison the two approaches is the duality between stochastic differential equations and Fokker- Planck equations, where ours is a (softly) discretized analog of the latter. Nonetheless, in many of the cases consider, Bubblewrap produces superior results even for state prediction. Nonetheless, like many similar models, ours includes multiple hyperparameters that require setting. While we did not experience catastrophic failure or sensitive dependence on parameters in our testing, and while our methods adapt to the scale and drift of the data, some tuning required in practice.
As detailed above, while many methods target population dynamics, and few target closed-loop settings 31 16. 55), very few models are capable of being trained online. Thus, the most closely related approaches are those in 30 [2] to which we provide extensive comparisons. However, these comparisons are somewhat strained by the fact that provided all models with the same proSVD-reduced ciw-dimensional data, while [22] is capable of modeling high-dimensional data in its own right and 30 was targeted at inferring neural computations from dynamical systems. We thus view this work as complementary to the dynamical systems approach, one that may be preferred when small distinctions among population dynamics are less important than characterizing highly noisy, non-repcating neural behavior.
Finally, we showed that online training of Bubblewrap can be performed fast enough for even kiloHertz. data acquisition rates if small latencies are tolerable and gradient steps can performed for small numbers of samples at a time. Yet, for real-time applications, it is not training time but the

---

Monkey reach
Widefield
Mouse video
Neuropixels
Log pred prob
-20
-12
-5-
40
-15
Bubblewrap
-14
-25
VJF
60
ZP (2016)
10
Entropy
4
5
4
0
0
1
10
1
10
5 10
1
5
10
Steps ahead
Steps ahead
Steps ahead
Steps ahead
Figure 5: Multi-step ahead predictive performance. (top) Mean log predictive probability as function of the number of steps ahead used for prediction for each of the four experimental datasets studied. Colors indicate model. (hottom) Bubblewrap entropy as function of the number of steps ahead used for prediction. Higher entropy indicates more uncertainty about future states. Dashed lines denote maximum entropy for each dataset (log of the number of tiles).
time required to make predictions that is relevant, and demonstrate prediction times of tens of microseconds. Moreover, Bubblewrap is capable of producing effective predictions multiple time steps into the future, providing ample lead time for closed-lood interventions. Thus, coarse- graniiiggiinng methods like ours open the door to online manipulation and steering of neural systems.
Acknowledgments and Disclosure of Funding
Research reported in this publication was supported by NIH BRAIN Initiative Planning Grant (R34NS116738; JP), and a Swartz Foundation Postdoctoral Fellowship for Theory in Neuroscience (AD). AD also holds Career Award at the Scientific Interface from the Burroughs Wellcome Fund.
References
[1] Misha В Ahrens, Michael B Orger, Drew N Robson, Jennifer M Li, and Philipp J Keller. Whole-brain functional imaging at cellular resolution using light-sheet microscopy. Nature methods, 10(5):413 420, 2013.
[2] Valentina Emiliani, Adam E Cohen, Karl Deisseroth, and Michael Häusser. All-optical interro- gation of neural circuits. Journal of Neuroxcience, 35(41):13117-1396 2015.
[3] Ian Stevenson and Konrad Kording. How advances in neural recording affect data analysis. Nature neuroscience, 14(2):199-142, 2011.
[4] Nicholas A Steinmetz, Christof Koch, Kenneth D Harris, and Matteo Carandini. Challenges and opportunities for large-scale clectrophysiology with neuropixels probes. Current opinion in neurobiology, 50:92-100, 2018.
[5] Nicholas 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: miniaturized high-density probe for stable, long-term brain recordings. bioRxiv, 2020.
[6] Valerio Mante, David Sussillo, Krishna V Shenoy, and William T Newsome. Context-dependent computation by recurrent dynamics in prefrontal cortex. nature, 503(7474):78 84,888847788 2013.
Kanaka Rajan, Christopher Harvey, and David Tank. Recurrent network models of
sequence generation and memory. Neuron, 90(1): 2016.
10

---

[8] Francis Song, Guangyu R Yang, and Xiao-Jing Wang. Training inbiiiiio recurrent neural networks for cognitive tasks: a simple and flexible framework. PLoS computational biology, 1202)ee0 2016.
[9] John P Cumningham and M Yu Byron. Dimensionality reduction for large-scale neural record- 2014.
ings. Nature neuroscience,
[10] Becket Ebitz and Benjamin Y Hayden. The population doctrine revolution in cognitive neurophysiology. arXiv preprint arXiv:2/04000144 2021.
[11] Eric T Carlson, Russell Rasquinha, Kechen Zhang, and Charles E Connor. A sparse object coding scheme in area v4. Current Biology, 21(4):288-299 2011.
[12] Christopher DiMattina and Kechen Zhang. Adaptive stimulus optimization for sensory systems neuroscience. Frontiers in neural circuits, 7:101, 2013.
[13] Benjamin Cowley, Ryan Williamson, Katerina Clemens, Matthew Smith, and Byron M Yu. Adaptive stimulus selection for optimizing neural population responses. In Advances Neural Information Processing Systems, volume 30, 2017.
[14) Reza Abbasi-Asl, Yuansi Chen, Adam Bloniarz, Michael Oliver, Ben DB Willmore, Jack I. Gallant, and Bin Yu. The deeptune framework for modeling and characterizing neurons in visual cortex area v4. bioRxiv, page 465534, 2018.
[15] Zihui Zhang, Lloyd E Russell, Adam M Packer, Oliver M Gauld, and Michael Häusser. Closed- loop all-optical interrogation of neural circuits in vivo. Nature methods, 15(12):1037-104) 2018.
[16] Michael F Bolus, Adam A Willats, Christopher J Rozell, and Garrett B Stanley. State- space op- timal feedback control of optogenetically driven neural activity. Journal of Neural Engineering, 2020.
[17] Diogo Peixoto, Jessica R Verhein, Roozbeh Kiani, Jonathan C Kao, Paul Nuyujukian, Chan- dramouli Chandrasekaran, Julian Brown, Sania Fong, Stephen I Ryu, Krishna V Shenoy, et al. Decoding and perturbing decision states in real time. Nature, pages 1-6, 2021.
[18) Anne Draelos and John Pearson. Online neural connectivity estimation with noisy group testing. Advances in Neural Information Processing Systems, 33, 2020.
[19] Matthew Brand. Incremental singular value decomposition of uncertain data with missing values. Computer Vision? 2002, pages 707-720, 2002.
[20] Christopher G Baker. A block incremental algorithm for computing dominant singular subspaces. Master's thesis, Florida State University, 2004.
[21] Matthew Brand. Fast low-rank modifications of the thin singular value decomposition. Linear algebra and its applications, 415(1):20 30, 2006.
[22] Christopher G Baker, Kyle A Gallivan, and Paul Van Dooren. Low-rank incremental methods for computing dominant singular subspaces. Linear Algebra and its Applications, 436(8):2866 2888, 2012.
[23] Julien Mairal, Francis Bach, Jean Ponce, and Guillermo Sapiro. Online learning for matrix factorization and sparse coding. Journal of Machine Learning Research, 11(1), 2010.
[24] Evan Archer, Il Memming Park, Lars Buesing, John Cunningham, and Liam Paninski. Black 2015.
box variational inference for state space models. arXiv preprint
[25] Yuanjun Gao, Evan Archer, Liam Paninski, and John P Cunningham. Linear dynamical neural population models through nonlinear embeddings. arXiv preprint arXiv: 1605008444 2016.
[26] Chethan Pandarinath, Daniel J O' Shea, Jasmine Collins, Rafal Jozefowicz, Sergey D Stavisky, Jonathan C Kao, Eric M Trautmann, Matthew T Kaufman, Stephen Ryu, Leigh R Hochberg, et al. Inferring single-trial neural population dynamics using sequential auto-encoders. Nature methods, 15(10))805-815 2018.
11

---

[27] Scott Linderman, Matthew Johnson, Andrew Miller, Ryan Adams, David Blei, and Liam Paninski. Bayesian learning and inference recurrent switching linear dynamical systems. In Artificial Intelligence and Statistics, pages 914-922, 2017.
[28] Scott W. Linderman Matthew J. Johnson*, Andrew C. Miller, Ryan P. Adams, David M. Blei, and Liam Paninski. Bayesian learning and inference in recurrent switching linear dynamical systems. In Proceedings of the 20th International Conference on Artificial Intelligence and Statistics (AISTATS), 2017.
[29] Josue Nassar, Scott W Linderman, Monica Bugallo, and Memming Park. Tree-structured recurrent switching linear dynamical systems for multi-scale modeling. arXiv preprint arXiv: 12386, 2018.
[30] Yuan Zhao and II Memming Park. Interpretable nonlinear dynamic modeling of neural trajec- tories. In D. Lee, M. Sugiyama, U. Luxburg, I. Guyon, and R. Garnett, editors, Advances in Neural Information Processing Systems, volume 29. Curran Associates, Inc., 2016.
[31] Yuxiao Yang, Allison T Connolly, and Maryam M Shanechi. A control-theoree system identification framework and real-time closed-loop elinical simulation testbed for electrical brain stimulation. Journal of neural engineering, 15(6):066007, 2018.
[32] Yuan Zhao and II Memming Park. Variational online learning of neural dynamics. Frontiers in computational neuroscience, 14, 2020.
[33] Peiran Gao, Eric Trautmann, Byron Yu, Gopal Santhanam, Stephen Ryu, Krishna Shenoy, and Surya Ganguli. A theory of multineuronal dimensionality, dynamics and measurement. BioRxiv, page 214262.
[34] Eric M Trautmann, Sergey D Stavisky, Subhaneil Lahiri, Katherine Ames, Matthew Kaufman, Daniel J O' Shea, Saurabh Vyas, Xulu Sun, Stephen Ryu, Surya Ganguli. al. Accurate estimation of neural population dynamics without spike sorting. Neuron, 103(2):292- 308, 2019.
[35] Dimitris Achlioptas. Database-friendly random projections: Johnson--lndeentrrass with binary coins. Journal of computer and System Sciences, 66(4):671-687 2003.
[36] Ping Li, Trevor J Hastie, and Kenneth W Church. Very sparse random projections. In Proceed- ings of the 12th ACM SIGKDD international conference on Knowledge discovery and data mining, pages 287 296, 2006.
[37] William B Johnson and Joram Lindenstrauss. Extensions of lipschitz mappings into hilbert space. Contemporary mathematics, 26( (61189-206):,, 1984.
[38] Peter H Schönemann. generalized solution of the orthogonal procrustes problem. Psychome- Irika, 31(1):1-10, 1966.
[39] Alan D Degenhart, William E Bishop, Emily R Oby, Elizabeth C Tyler-Kabara, Steven M Chase, Aaron P Batista, and Yu Byron. Stabilization of brain compucr interface via the alignment of dimedinnlli spaces of neural activity. Nature biomedical engineering, 6(7):672-685. 2020.
[40] David A Ross, Jongwoo Lim, Ruei-Sung Lin, and Ming-Hsuan Yang. Incremental learning for robust visual tracking. International journal of computer vision, 77(1):125-14) 2008.
[41] Gordon J Berman, Daniel M Choi, William Bialek, and Joshua W Shaevitz. Mapping the stereotyped behaviour of freely moving fruit flies. Journal of The Royal Society Interface, 11099)20140672, 2014.
[42] Gordon J Berman, William Bialek, and Joshua W Shaevitz. Predictability and hierarchy in drosophila behavior. Proceedings of the National Academy of Sciences, (9(42):11443-1999 2016.
[43] Talmo Pereira, Joshua Shaevitz, and Mala Murthy, Quantifying behavior to understand the brain. Nature Neuroscience, pages 1-13, 2020.
12

---

[44] Gianluigi Mongillo and Sophie Deneve. Online learning with hidden markov models. Neural computation, 2007):1706-1716, 2008.
[45] Olivier Cappé and Eric Moulines, On-line expectation-maximizaii algorithm for latent data models. Journal of the Royal Statistical Society: Series B (Statistical Methodology), 71(3):593 613, 2009.
[46] Sylvain Le Corff, Gersende Fort, et al. Online expectation maximization based algorithms for inference in hidden markov models. Electronic Journal of Statistics, 7:763-792, 2013.
[47] Diederik P Kingma and Jimmy Ba. Adam: A method for stochastic optimization. arXiv preprint arXiv:1412.6980, 2014.
[48] Mark Churchland. Churchland lab code.
columbia. wwwmmbii..edd/ottett////ooo
[49] Mark M Churchland, John P Cunningham, Matthew Kaufman, Justin Foster, Paul Nuyu- jukian, Stephen I Ryu, and Krishna V Shenoy. Neural population dynamics during reaching. Nature, 4 2012.
[50] Simon Musall, Matthew T Kaufman, Ashley L Juavinett, Steven Gluf, and Anne Churchland. Single-trial neural dynamics are dominated by richly varied movements. Nature neuroscience, 22(10):1777-166, 2019.
[51] Simon Musall, Matthew T. Kaufman, Ashley L. Juavinett, Steven Gluf, and Anne K. Churchland, Single-trial neural dynamics are dominated by richly varied movements: dataset. Technical report, October 2019.
[52] Nick Steinmetz, Marius Pachitariu, Carsen Stringer, Matteo Carandini, and Kenneth Harris. Eight-probe neuropixels recordings during spontaneous behaviors, Mar 2019.
[53] Carsen Stringer, Marius Pachitariu, Nicholas Steinmetz, Charu Bai Reddy, Matteo Carandini, and Kenneth D Harris. Spontaneous behaviors drive multidimensional, brainwide activity. Science, 364(6437), 2019.
[54] James Bradbury, Roy Frostig, Peter Hawkins, Matthew James Johnson, Chris Leary, Dougal Maclaurin, George Necula, Adam Paszke, Jake VanderPlas, Skye Wanderman-Minn and Qiao Zhang. JAX: composable transformations of Python+Nummy programs, 2018.
[55] Omid G Sani, Hamidreza Abbaspourazad, Yan T Wong, Bijan Pesaran, and Maryam M Shanechi. Modeling behaviorally relevant neural dynamics enabled by preferential subspace identification. Technical report, Nature Publishing Group, 2020.
13


Reviews:
[{"review_id": "jxkgorValQ-", "reviewer": "Reviewer_Xhub", "summary": "This paper presents a method to efficiently characterize high-dimensional neural dynamics using a soft \"tiling\" of a lower-dimensionality manifold. Specifically, the method is a combination of two algorithms, (i) \"Stable Streaming SVD\", a dimensionality reduction method that aims to find the top-k-dimensional subspace instead of the top-k bases; (ii) \"BubbleWrap\", a Gaussian mixture model of the probability density with Hidden Markov Model transitions, that works on the low-dimensional manifold generated in (i). This is a model-free approach, in the sense that it is not based on a model of the dynamical system, and the authors claim that the flexibility of the method makes it outperform other (model-based) approaches in high-noise regimes. The authors also highlight that the proposed inference is fast, possibly allowing online applications.", "questions": "", "limitations": "", "rating": 6, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "This work addresses an important and interesting problem in the analysis of neural data: a first-pass, low-dimensional characterization of the neural dynamics, especially when the underlying dynamical system is not well understood and/or when it is important to have a rough picture of the dynamics with a limited time budget. With the flexibility of the model formulation and the fast inference, the method has a potential of being a useful tool in the field.\n\nI am enthusiastic about the first part, ssSVD; this could be a very useful method for characterizing low-dimensional dynamics. The idea of finding a stable subspace, instead of the specific basis set, is simple but effective. The demonstration in Fig 1 is also impressive.\n\nI think the \"BubbleWrap\" algorithm is an interesting take to the problem, with a novel combination of ideas in the context of neural dynamics. However, the presented results are less convincing. In particular, it is not clear what are lessons from the analyses (and what is the potential use case for the method), or how the presented results support the usefulness of the method. My major concerns are below.\n\nMajor concerns:\n\n- Log predictive probability is a key quantity in the analysis in the Experiments section, but it is not sufficiently described in the main paper. (It is understandable to keep the details in the supplements, but the main text should be self consistent.) In particular, it is not clear how one should interpret the log predictive probability. What is the message that should be read out of Table 1? Is the absolute value (mean) meaningful, or should one only compare values across different methods within each row? Also in Fig 3d-f and Fig 4b, the values are not really meaningful without having a sense of \"how good is good\".\n\n- The authors already note in the paper that the proposed method takes flexibility at the cost of reduced accuracy. But in this case, it should be discussed how the (possibly inaccurate) outcome of this method should be used towards a final understanding of the neural data. For example, in Fig 3a, although the highlighted tiles (red) roughly overlap with the trajectory (blue), there are considerable discrepancy, especially for the lower trajectory; if you do any downstream analysis based on this representation, it is inevitable that these analyses will be flawed. If BubbleWrap is meant to be a preliminary analysis, what follow-up analyses could it guide?\n\n- Related to the accuracy issue: the \"bubble\" structure, due to the Gaussian-mixture nature of the model, introduces an artificial heterogeneity to the probability density. It would be necessary to discuss the effect of these \"bubbles\" on further analysis based on the resulting probability, and/or ways to remove this artifact by adding a postprocessing step.\n\n- Fig 3b: Text description (L180) says that this dataset was chosen to demonstration a case with rapid transitions between noise clusters. Because the model explicitly infers transition probabilities, it is important that they are shown and discussed in this case, not \"omitted for clarity\". Speaking of the transition probabilities, it seems that the black lines in Fig 3a do not always correctly capture the trajectory.\n\nMinor comments:\n\n- L107 and Fig 1c: it is not clear how to read out the tradeoff value around n=200 from this figure. In fact, the color-coded points looks like n is transitioning from somewhere around 10 (blue, first point) to around 500 (red, second point) with no intermediate points in between. In general, please indicate the actual values of n in the figure if the value is relevant.\n\n\n**EDIT**: The author response to my comments, as well as to the comments from the other reviewers, helped clarify my original concerns, which were:\n\n> it is not clear what are lessons from the analyses (and what is the potential use case for the method), or how the presented results support the usefulness of the method.\n\nAs for the first point (use case), it was crucial for me to understand that the target use case of this set of methods was close-loop experiments, which I missed in my initial reading of the paper. Now I better appreciate the potential usefulness of the method in real time application.\n\nAs for the second point (results), the authors addressed my specific questions in the response. Also importantly, as I was reading the responses to the other reviewers, I realized that I was misunderstanding the \"entropy\" in the BubbleWrap results as the entropy of the state-space (x) distribution, where it was the entropy of the predictive conditional (z|x). This makes a lot of sense and it allows me to re-evaluate the results in the paper. Clarification on the HMM/GMM terminology was also helpful.\n\nGiven this, I am happy to increase my rating from 4 to 6, and now vote for acceptance, under the assumption that the authors will describe these potentially confusing points more clearly in the final version of the manuscript.", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "The authors discussed limitations of their approach, including the decreased accuracy as a tradeoff to gain flexibility. \n\nI do not see any potential negative societal impact of this work.", "review_text": "", "extra_scores": {"ethical_concerns": "None noted.", "needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "5"}, "extra_text": {}}, {"review_id": "ZKZFk7eatlF", "reviewer": "Reviewer_oLEL", "summary": "The authors present an online method for estimating low-dimensional neural dynamics that scales well to large observation dimensions and high sampling rates. The first step of their method is an online version of SVD, designed so that it converges to a stable subspace with very little training data. The second step of their method is an online version of an HMM with a large number of states that tile the low-dimensional space provided by SVD. They demonstrate their method on both simulated and real data, demonstrating both its speed and accuracy.", "questions": "", "limitations": "", "rating": 7, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "originality: this work is based on a growing literature of online algorithms designed to capture neural dynamics; the contribution here is a semiparametric model that is not based on a dynamical system (such as LDS), and the authors show their model works better in the high noise regime. \n\nquality: the model is well-motivated and sources methods from a variety of disciplines (random matrix theory, linear algebra, graphical models). The main feature of the model is its online ability to fit data, and with this in mind the authors do a good job of benchmarking the model, as well as comparing it to relevant baselines. However, I am still left wondering what a potential use-case for this method is - perhaps the authors could give an example in the introduction?\n\nclarity: the manuscript is well written and the model is well-motivated. One point of confusion for me was L123-127: the authors mention both a GMM and an HMM, and I first I thought we were dealing with an HMM with GMM emissions; however, it appears that the emissions are standard Gaussians. Maybe the GMM terminology could be dropped in this case?\n\nsignificance: online estimation of neural dynamics/state is crucial for causal experiments, and as such successful methods will likely be widely adopted. This method seems well suited to such estimation, though the authors have not shown if their method is accurate enough to enable the targeted interventions required for a specific experiment.", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "The authors have addressed some of the limitations in the discussion. One thing I'm curious about that did not seem to be addressed was the sheer number of states in the HMM - >1k!! In L49 there is mention of a \"sparse transition graph\", but I found no further mention of this in the main text or the supplemental. It would be useful if the authors could briefly mention how they handle such a large transition matrix, both during training and inference.", "review_text": "", "extra_scores": {"needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "1.5"}, "extra_text": {}}, {"review_id": "Wwz9Dq6q4BT", "reviewer": "Reviewer_KUqF", "summary": "The authors present Bubblewrap which is a method for dimensionality reduction and discrete approximation of dynamics for high dimensional timeseries datasets. Bubblewrap stabilizes incremental basis estimates into an incremental subspace estimate and the authors show that this converges quickly in practice. The authors then use an HMM with Gaussian emissions to model the projected data along with several heuristics to achieve better performance. The method in benchmarked in terms of likelihood and latency/data rate on synthetic and real neural and behavioral data.", "questions": "", "limitations": "", "rating": 6, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "The manuscript is clearly written and methods that can be run in online during high dimensional neural recordings are important for closed-loop experiments. The reported results are clear and although I have not run the code, it looks complete and uses publicly available datasets.\n\nI'm not an expert in closed-loop neuroscience experiments, but it seems like methods might potentially trade-off\n1. latency,\n2. the ability to predict the future states of the data Bubblewrap is modeling (e.g., neural states for optogenetic manipulation), or\n3. the ability to predict paired timeseries like behavior (e.g., for robotic or task manipulations during neural recordings).\n\nRun-time comparisons for 1. versus refs [30, 32] could be made. 2. is analyzed clearly for this method.\n\nIn my opinion, the main area where the manuscript could be strengthened is by including some evaluation of the utility of the inferred discrete states (i.e., 3.). Can the behavioral state of the animal be predicted? The increased log predictive probability compared to other online methods is promising, but some analysis of the ability to predict behavior would make the method more compelling.\n\nGeneral comments/questions\n- Lines 126-127: I think many people would consider a HMM with Gaussian emissions to be a state-space model with a noise distribution. Can this sentence be made more precise? The related sentence on lines 214-215 could also be made more clear, I believe the continuous versus discrete dynamics distinction is more salient than moment-to-moment prediction versus tiling. Ultimately, Bubblewrap also learns the tiling through maximizing the data likelihood, which is another way to say prediction. I think the sentence on lines 120-122 has a more clear statement of the comparison.\n- Could you provide some interpretation of the plots of entropy in Figs 3 and 4? Is it meaningful that Fig 3 d and f drop closer to zero while Fig 3 e stays closer to the max? Does it tell us something about the neural dynamics?\n- In Fig S2, what is the grey in the background of the right column mean?\n- Lines 225-226: can you expand on what the most important parameters to tune were, if any?\n- The caption for the plots of log predictive probability and entropy don't describe what the black points are? Smoothed average?\n\nMinor comments/questions\n- Is there a reason that sparse random projections are used for the initial projection rather than dense random projections? Is there a computational benefit?\n- For Fig 4 c, having stacked bars and a log-scale y axis makes it difficult to compare times if they don't happen to be anchored to the same lower point. Having 4 bars (3 parts + total time) grouped side-by-side for each case may help. I would have also expected the \"Prediction\" time to scale with N, am I misunderstanding something?", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "I would like to see more analysis or discussion about how the choice of discrete inferred states may impact potential uses for this method. Can the sequence of inferred centroids be used? Are the centroids or networks of centroids related to behavior?", "review_text": "", "extra_scores": {"needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "5"}, "extra_text": {}}]

Author rebuttal:

Meta-review: After discussion, the reviewers all converged on accept or weak accept ratings for this paper, once the stated changes are implemented. The methodological novelty is somewhat limited, but the problem is important, the approach is sensible, and the results are good. 
Decision: Accept (Poster)