Bubblewrap: Online tiling and real-time flow prediction on neural manifolds
Year: 2022

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.

heading: 1 Introduction
text: Systems neuroscience is in the midst of a data explosion. Advances in microscopy [1, 2] and probe technology [3, 4, 5] have made it possible to record thousands of neurons simultaneously in behaving
35th Conference on Neural Information Processing Systems (NeurIPS 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, 7, 8].
At the same time, this increase in data volume has resulted in increasing emphasis on methods for dimensionality reduction [9] and a focus on neural populations in preference to the coding properties of individual neurons [10]. 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 [11, 12, 13, 14], in closed-loop designs [15, 16, 17], and even to scale up holographic photostimulation for inferring functional connectivity in large circuits [18].
Yet, despite their promise, adaptive methods are rarely applied in practice for two reasons: First, although efficient online methods for dimensionality reduction exist [19, 20, 21, 22, 23], 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 [6, 24, 25, 26, 27, 28, 29], including online approaches [30, 31, 16, 32], they typically assume a 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 ⇡10ms using a GPU, while a cached copy of the model in main memory is capable of predicting upcoming states at <1ms latency. As a result, Bubblewrap offers a method performant and flexible enough to serve as a neural prediction engine for causal feedback experiments.

heading: 2 Stable subspaces from streaming SVD
text: 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, 5], tens of thousands of pixels for calcium imaging. However, as theoretical work has shown [33, 34], 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) [35, 36]. 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 1n in dimension n, and so we desire n to be as small as possible. However, by the Johnson-Lindenstrauss Lemma [37, 36], when reducing from d to 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
(1 ")ku vk2  ku⇤ v⇤k 2  (1 + ")ku vk2 (1)
with probability 1 if n > O(log(1/ )/"2). 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.
Algorithm 1 Procrustean SVD (proSVD) 1: Given: Initial data X0, decay parameter ↵ 2 (0, 1] 2: Initialize: QR Factorization: X0 = Q0R0 3: 4: for t = 1 . . . do 5: Fetch b new columns of data, X+ 6: C Q>t 1X+, X? X+ Qt 1C, Q?, R? QR(X?) . Gram-Schmidt
7: Q̂ [Qt 1 Q?], R̂  Rt 1 C 0 R?
. QR of augmented data
8: U,⌃, V SVD(R̂) 9: ⌃ ↵⌃ . Discount old data 10: M Q>t 1Q̂U1 = [ k⇥k 0k⇥b]U1 . U1 contains the first k columns of U 11: Ũ , ⌃̃, Ṽ SVD(M), T Ũ Ṽ > . Orthogonal Procrustes: minT kQ̂U1T> Qt 1kF 12: Qt Q̂U1T> . Update left subspace basis 13: Qv, Rv QR(V ), Rt T⌃1Q>v . QR right subspace, update inner block 14: end for
Thus, in the second stage, we reduce from n ⇠ O(100) to k ⇠ O(10) dimensions using a streaming singular value decomposition. This method is based on the incremental block update method of [20, 22] with an important difference: While the block update method aims to return the top-k SVD at every time point, the directions of the singular vectors may be quite variable during the course of an experiment (Figure 1d–h), which implies an 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 a 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 xt 2 Rn be a vector of input data after random projections. In our streaming setup, these are processed b samples at a time, with b = 1 reasonable for slower methods like calcium imaging and b = 40 more appropriate for electrophysiological sampling rates of ⇠20kHz. Then, if the data matrix X has dimension n ⇥ T , adding columns over time, the incremental method of [20, 22] produces at each time step a factorization X = 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 be equivalent to the SVD. In the incremental algorithm, R is augmented at each timestep based on new data to form R̂, which is block diagonalized via an orthogonal matrix and truncated to the top-k subspace, allowing for an exact reduced-rank SVD (Appendix A).
However, as reviewed in [20, 22], since there are multiple choices of basis Q for for the top-k singular subspace, there are likewise multiple choices of block diagonalization for R̂. In [20, 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:
minkQt Qt 1kF = min T kQ̂U1T
> Qt 1kF , (2)
where Q̂ is an augmented basis for the top-(k + b) singular subspace, U1 contains the first k left singular vectors of R̂, and T is an orthogonal matrix (Appendix A). This minimization is known as the Orthogonal Procrustes problem and has a well-known solution [38]: T = Ũ Ṽ >, where Ũ and Ṽ are the left and right singular vectors, respectively, of M ⌘ Q>t 1Q̂U1. (See [39] for a recent application of similar ideas in brain-computer interfaces). This Procrustean SVD (proSVD) procedure is summarized in Algorithm 1. 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 as in [40].
Figures 1a-c illustrates the performance of the two-stage dimension reduction algorithm for a case of d = 104 randomly generated Gaussian data. While proSVD yields minimal distortion (due to truncation of the spectrum to k = 6), random projections require k ⇠ O(100) to achieve the same result (Figure 1a). By contrast, random projections are much faster (Figure 1b). Thus, we can trade off distortion against time by adjusting n, the number of intermediate dimensions. As Figure 1c shows, the optimal tradeoff occurs somewhere around n = 200 for this example.
Figures 1d-i show results for neural data from recorded from monkey motor cortex [26] in a cued reach task. While projection of 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 (⇡15s, middle). This is due to the fact that, in all data sets we examined, the top-k SVD subspace was identified extremely quickly; proSVD simply ensures the choice of a stable basis for that subspace.

heading: 3 Bubblewrap: a soft manifold tiling for online modeling
text: As reviewed above, most neural population modeling approaches are based on the dynamical systems framework, assuming a lawful equation of motion corrupted by noise. However, for animals engaged in task-free natural behavior [41, 42, 43], trajectories are likely to be sufficiently complex that simple dynamical models fail. For instance, dynamical systems models with Gaussian noise necessarily produce unimodal transition probabilities centered around the mean prediction, while neural trajectories may exhibit multimodal distributions beginning at the same system state. By contrast, we pursue an
alternative method that trades some accuracy in estimating instantaneous system state for flexibility in modeling the manifold describing neural activity.
Our approach is to produce a soft tiling of the neural manifold in the form of a Gaussian mixture model (GMM), each component of which corresponds to a single tile. We then approximate the transitions 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 an increasingly finer-grained description of dynamics that assumes neither an underlying dynamical system nor a particular distribution of noise.
More specifically, let xt be the low-dimensional system state and let zt 2 1 . . . N index the tile to which the system is assigned at time t. Then we have for the dynamics
p(zt = j|zt 1 = i) = Aij p(xt|zt) = N (µzt ,⌃zt) p(µj ,⌃j) = NIW(µ0j , j , j , ⌫j), (3)
where we have assumed Normal-inverse-Wishart priors on the parameters of the Gaussians. Given its exponential family form and the conjugacy of the priors, online expectation maximization updates are available in closed form [44, 45, 46] for each new datum, though we opt, as in [45] for a gradient-based optimization of an estimate of the evidence lower bound
L(A, µ,⌃) = X
ij
(N̂ij(T ) + ij 1) logAij + X
j
(Ŝ1j(T ) + jµ0j) >⌃ 1j µj (4)
1
2
X
j
tr(( j + Ŝ2j(T ) + jµ0jµ > 0j + ( j + n̂j(T ))µjµ > j )⌃ 1 j )
1
2
X
j
(⌫j + n̂j(T ) + d+ 2) log det⌃j
with accumulating (estimated) sufficient statistics
↵j(t) = X
i
↵i(t 1) ij(t) N̂ij(t) = (1 "t)N̂ij(t 1) + ↵i(t 1) ij(t) (5)
n̂j(t) = X
i
N̂ij(t) Ŝ1j(t) = (1 "t)Ŝ1j(t 1) + ↵j(t)xt
Ŝ2j(t) = (1 "t)Ŝ2j(t 1) + ↵j(t)xtx > t
where ↵j(t) = p(zt = j|x1:t) is the filtered posterior, ij(t) is the update matrix from the forward algorithm [44], and "t is a 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 ↵(t) and not the posterior over all observed data.
In setting Normal-Inverse-Wishart priors over the Gaussian mixture components, we take an empirical Bayes approach by setting prior means µ0j to the current estimate of the data center of mass and prior covariance parameters j to N 2 k times the current estimate of the data covariance (Appendix B). For initializing the model we use a small data buffer M ⇠ O(10). We chose effective observation numbers ( , ⌫) = 10 3 and trained this model to maximize L(A, µ,⌃) using Adam [47], enforcing parameter constraints by replacing them with unconstrained variables aij and lower triangular Lj with positive diagonal: Aij = exp(aij)/ P j exp(aij), ⌃ 1 j = LjL > j .
Finally, in order to prevent the model from becoming stuck in local minima and to 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 p(xt|zt) < ✓n for all zt), we teleported a node at this data point by assigning ↵J(t) = 1 for an unused index J . For initial learning this results in a “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 n̂(t) 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.
Algorithm 2 Bubblewrap 1: Given: Hyperparameters j , ⌫j , t, forgetting rate "t, teleport threshold ✓, step size , initial
data buffer M 2: Initialize with {x1 . . . xM}: µj µ̄, ⌃j ⌃̄, aij 1N , ↵j j 1 N . 3: 4: for t = 1 . . . do 5: Observe new data point xt. 6: if log p(xt|zt) < ✓ 8zt then . Teleport 7: µJ = xt,↵J(t) = 1 for J = argminj n̂j(t) 8: end if 9: Calculate ij(t) via forward filtering [44]. 10: Update sufficient statistics via (5). . E step 11: µ̄ P
j Ŝ1jP j n̂j
, ⌃ P
j Ŝ2jP j n̂j µ̄µ̄ T . Global mean and covariance update
12: ✏j ⇠ N (0, ⌘2), µ0j aµ0j + (1 a)µ̄+ ✏j , j ⌃ N 2 k . Update priors (Appendix B) 13: Perform gradient-based update of L(A, µ,⌃) (4) . M step 14: end for

heading: 4 Experiments
text: We demonstrated the performance of Bubblewrap on both simulated non-linear dynamical systems and experimental neural data. We compared these results to two existing online learning models for neural data, both of which are based on dynamical systems [30, 32]. To simulate low-dimensional systems, we generated noisy trajectories from a two-dimensional Van der Pol oscillator and a threedimensional Lorenz attractor. For experimental data, we used four publicly available datasets from a range of applications: 1) trial-based spiking data recorded from primary motor cortex in monkeys performing a reach task [48, 49] preprocessed by performing online jPCA [49]; 2) continuous video data and 3) trial-based wide-field calcium imaging from a rodent decision-making task [50, 51]; 4) high-throughput 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, we quantified overall model performance by taking the mean log predictive probability over the last half of each data set (Table 1). For Bubblewrap, prediction T steps into the future gives
log p(xt+T |x1:t) = log X
i,j
p(xt+T |zt+T = j)p(zt+T = j|zt = i)p(zt = i|x1:t)
= log X
i,j
N (xt+1;µj ,⌃j)(A T )ij↵i(t), (6)
where AT is the T -th power of the transition matrix. Conveniently, these forward predictions can be efficiently computed due to the closed form (6), while similar predictions in comparison models [30, 32] 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: H(t, T ) = X
j
p(zt+T = j|x1:t) log p(zt+T = j|x1:t) = X
ij
(AT )ij↵i(t) log X
k
(AT )kj↵k(t).
(7)
Additional detailed experimental results and benchmarking of our GPU implementation in JAX [54] are in Appendix D. 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://github.com/pearsonlab/Bubblewrap.
When tested on low-dimensional dynamical systems, Bubblewrap successfully learned tilings of both neural manifolds, outperforming VJF [32] on both datasets (Figure 2a,b) while it was comparable to the 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
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 3a) to rapid transitions between noise clusters (Figure 3b) to slow dynamics (Figure 3c). In each case, Bubblewrap learned a tiling of the data that allowed it to equal or outperform state predictions from the comparison algorithms (Figure 3d–f, blue). In some cases, as with the mouse dataset, the algorithm of [30] produced predictions for xt by degenerating to a random walk model (Table 1 marked with *; Appendix D). Regardless, Bubblewrap’s tiling generated transition predictions with entropies far below those of a random walk (Figure 3d–f, green), indicating it successfully identified coarse structure, even in challenging datasets. Thus, even though these data are noise-dominated and lack much of the typical structure identified by neural population models, coarse-graining identifies some reliable patterns.
We additionally considered the capability of our algorithm to scale to high-dimensional or highsampling 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 1). More importantly, all these steps can be performed well within the 30ms per sample time of the data (c). In fact, when testing on representative examples of d = 104 dimensions, 1 kHz sampling rates, or N = 20, 000 tiles, our algorithm was able to maintain amortized per-sample processing times below those of data acquisition. In practice, we
found that even in higher-dimensional datasets (as in the Neuropixels case), only 1-2 thousand tiles were used by the model, making it easy to run at kHz data rates. What’s more, while learning involved 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.
Just as importantly, when used for closed loop experiments, algorithms must be able to produce predictions far enough into the future for interventions to be feasible. Thus we examined the performance of our algorithm and comparison models for predicting T steps ahead into the future. Bubblewrap allows us to efficiently calculate predictions even many time steps into the future using (6), whereas the comparison models require much costlier sampling approaches. Figure 5 shows the 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.

heading: 5 Discussion
text: While increasing attention has been paid in neuroscience to population hypotheses of neural function [10], 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 engineering literature has long used online methods based on
Kalman filtering, (e.g., [16]), and these methods are known to work well in 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, to address this need, we have introduced both a 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 largescale probability flow is 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 3e) when it did not learn to predict and degenerated to a random walk. Indeed, the most relevant comparison to 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 we 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 was required in practice.
As detailed above, while many methods target population dynamics, and a 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, 32], to which we provide extensive comparisons. However, these comparisons are somewhat strained by the fact that we provided all models with the same proSVD-reduced low-dimensional data, while [32] 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-repeating 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 be performed for small numbers of samples at a time. Yet, for real-time applications, it is not training time but the
time required to make predictions that is relevant, and we 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-loop interventions. Thus, coarse-graining 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 a NIH BRAIN Initiative Planning Grant (R34NS116738; JP), and a Swartz Foundation Postdoctoral Fellowship for Theory in Neuroscience (AD). AD also holds a Career Award at the Scientific Interface from the Burroughs Wellcome Fund.

Review summaries:
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.

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.

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.