Optical calcium imaging is a powerful tool for recording neural activity across a wide range of spatial scales, from dendrites and spines to whole-brain imaging through two-photon and widefield microscopy. Traditional methods for analyzing functional calcium imaging data rely heavily on spatial features, such as the compact shapes of somas, to extract regions of interest and their associated temporal traces. This spatial dependency can introduce biases in time trace estimation and limit the applicability of these methods across different neuronal morphologies and imaging scales. To address these limitations, the Graph Filtered Temporal Dictionary Learning (GraFT) uses a graph-based approach to identify neural components based on shared temporal activity rather than spatial proximity, enhancing generalizability across diverse datasets. Here we present significant advancements to the GraFT algorithm, including the integration of a more efficient solver for the L1 least absolute shrinkage and selection operator (LASSO) problem and the application of compressive sensing techniques to reduce computational complexity. By employing random projections to reduce data dimensionality, we achieve substantial speedups while maintaining analytical accuracy. These advancements significantly accelerate the GraFT algorithm, making it more scalable for larger and more complex datasets. Moreover, to increase accessibility, we developed a graphical user interface to facilitate running and analyzing the outputs of GraFT. Finally, we demonstrate the utility of GraFT to imaging data beyond meso-scale imaging, including vascular and axonal imaging.

Modern recordings of neural activity provide diverse observations of neurons across brain areas, behavioral conditions, and subjects; presenting an exciting opportunity to reveal the fundamentals of brain-wide dynamics. Current analysis methods, however, often fail to fully harness the richness of such data, as they provide either uninterpretable representations (e.g., via deep networks) or oversimplify models (e.g., by assuming stationary dynamics or analyzing each session independently). Here, instead of regarding asynchronous neural recordings that lack alignment in neural identity or brain areas as a limitation, we leverage these diverse views into the brain to learn a unified model of neural dynamics. Specifically, we assume that brain activity is driven by multiple hidden global sub-circuits. These sub-circuits represent global basis interactions between neural ensembles—functional groups of neurons—such that the time-varying decomposition of these sub-circuits defines how the ensembles' interactions evolve over time non-stationarily and non-linearly. We discover the neural ensembles underlying non-simultaneous observations, along with their non-stationary evolving interactions, with our new model, **CREIMBO** (**C**ross-**R**egional **E**nsemble **I**nteractions in **M**ulti-view **B**rain **O**bservations). CREIMBO identifies the hidden composition of per-session neural ensembles through novel graph-driven dictionary learning and models the ensemble dynamics on a low-dimensional manifold spanned by a sparse time-varying composition of the global sub-circuits. Thus, CREIMBO disentangles overlapping temporal neural processes while preserving interpretability due to the use of a shared underlying sub-circuit basis. Moreover, CREIMBO distinguishes session-specific computations from global (session-invariant) ones by identifying session covariates and variations in sub-circuit activations. We demonstrate CREIMBO's ability to recover true components in synthetic data, and uncover meaningful brain dynamics in human high-density electrode recordings, including cross-subject neural mechanisms as well as inter- vs. intra-region dynamical motifs. Furthermore, using mouse whole-brain recordings, we show CREIMBO's ability to discover dynamical interactions that capture task and behavioral variables and meaningfully align with the biological importance of the brain areas they represent.

The brain is a highly complex organ consisting of a myriad of subsystems that flexibly interact and adapt over time and context to enable perception, cognition, and behavior. Understanding the multi-scale nature of the brain, i.e., how circuit- and moleclular-level interactions build up the fundamental components of brain function, holds incredible potential for developing interventions for neurodegenerative and psychiatric diseases, as well as open new understanding into our very nature. Historically technological limitations have forced systems neuroscience to be local in anatomy (localized, small neural populations in single brain areas), in behavior (studying single tasks), in time (focusing on specific stages of learning or development), and in modality (focusing on imaging single biological quantities). New developments in neural recording technology and behavioral monitoring now provide the data needed to break free of local neuroscience to global neuroscience: i.e., understanding how the brain’s many subsystem interact, adapt, and change across the multitude of behaviors animals and humans must perform to thrive. Specifically, while we have much knowledge of the anatomical architecture of the brain (i.e., the hardware), we finally are approaching the data needed to find the functional architecture and discover the fundamental properties of the software that runs on the hardware. We must take this opportunity to bridge between the vast amounts of data to discover this functional architecture which will face numerous challenges from low-level data alignment up to high level questions of interpretable mathematical models of behavior that can synthesize the myriad of datasets.

Time series data across scientific domains are often collected under distinct states (e.g., tasks), wherein latent processes (e.g., biological factors) create complex inter- and intra-state variability. A key approach to capture this complexity is to uncover fundamental interpretable units within the data, Building Blocks (BBs), which modulate their activity and adjust their structure across observations. Existing methods for identifying BBs in multi-way data often overlook inter- vs. intra-state variability, produce uninterpretable components, or do not align with properties of real-world data, such as missing samples and sessions of different duration. Here, we present a framework for Similarity-driven Building Block Inference using Graphs across States (SiBBlInGS). SiBBlInGS offers a graph-based dictionary learning approach for discovering sparse BBs along with their temporal traces, based on co-activity patterns and inter- vs. intra-state relationships. Moreover, SiBBlInGS captures per-trial temporal variability and controlled cross-state structural BB adaptations, identifies state-specific vs. state-invariant components, and accommodates variability in the number and duration of observed sessions across states. We demonstrate SiBBlInGS's ability to reveal insights into complex phenomena as well as its robustness to noise and missing samples through several synthetic and real-world examples, including web search and neural data.

Time-varying linear state-space models are powerful tools for obtaining mathematically interpretable representations of neural signals. For example, switching and decomposed models describe complex systems using latent variables that evolve according to simple locally linear dynamics. However, existing methods for latent variable estimation are not robust to dynamical noise and system nonlinearity due to noise-sensitive inference procedures and limited model formulations. This can lead to inconsistent results on signals with similar dynamics, limiting the model's ability to provide scientific insight. In this work, we address these limitations and propose a probabilistic approach to latent variable estimation in decomposed models that improves robustness against dynamical noise. Additionally, we introduce an extended latent dynamics model to improve robustness against system nonlinearities. We evaluate our approach on several synthetic dynamical systems, including an empirically-derived brain-computer interface experiment, and demonstrate more accurate latent variable inference in nonlinear systems with diverse noise conditions. Furthermore, we apply our method to a real-world clinical neurophysiology dataset, illustrating the ability to identify interpretable and coherent structure where previous models cannot.

Closed-loop neuroscience experimentation, where recorded neural activity is used to modify the experiment on-the-fly, is critical for deducing causal connections and optimizing experimental time. Thus while new optical methods permit on-line recording (via Multi-photon calcium imaging) and stimulation (via holographic stimulation) of large neural populations, a critical barrier in creating closed-loop experiments that can target and modulate single neurons is the real-time inference of neural activity from streaming recordings. In particular, while multi-photon calcium imaging (CI) is crucial in monitoring neural populations, extracting a single neuron's activity from the fluorescence videos often requires batch processing of the video data. Without batch processing, dimmer neurons and events are harder to identify and unrecognized neurons can create false positives when computing the activity of known neurons. We solve these issues by adapting a recently proposed robust time-trace estimator---Sparse Emulation of Unused Dictionary Objects (SEUDO) algorithm---as a basis for a new on-line processing algorithm that simultaneously identifies neurons in the fluorescence video and infers their time traces in a way that is robust to as-yet unidentified neurons. To achieve real-time SEUDO (realSEUDO), we introduce a combination of new algorithmic improvements, a fast C-based implementation, and a new cell finding loop to enable realSEUDO to identify new cells on-the-fly with no "warm-up" period. We demonstrate comparable performance to offline algorithms (e.g., CNMF), and improved performance over the current on-line approach (OnACID) at speeds of 120 Hz on average. This speed is faster than the typical 30 Hz framerate, leaving critical computation time for the computation of feedback in a closed-loop setting.

Learning interpretable representations of neural dynamics at a population level is a crucial first step to understanding how observed neural activity relates to perception and behavior. Models of neural dynamics often focus on either low-dimensional projections of neural activity or on learning dynamical systems that explicitly relate to the neural state over time. We discuss how these two approaches are interrelated by considering dynamical systems as representative of flows on a low-dimensional manifold. Building on this concept, we propose a new decomposed dynamical system model that represents complex non-stationary and nonlinear dynamics of time series data as a sparse combination of simpler, more interpretable components. Our model is trained through a dictionary learning procedure, where we leverage recent results in tracking sparse vectors over time. The decomposed nature of the dynamics is more expressive than previous switched approaches for a given number of parameters and enables modeling of overlapping and non-stationary dynamics. In both continuous-time and discrete-time instructional examples, we demonstrate that our model effectively approximates the original system, learns efficient representations, and captures smooth transitions between dynamical modes. Furthermore, we highlight our model’s ability to efficiently capture and demix population dynamics generated from multiple independent subnetworks, a task that is computationally impractical for switched models. Finally, we apply our model to neural “full brain” recordings of C. elegans data, illustrating a diversity of dynamics that is obscured when classified into discrete states.

The sparse coding model posits that the visual system has evolved to efficiently code natural stimuli using a sparse set of features from an overcomplete dictionary. The original sparse coding model suffered from two key limitations; however: (1) computing the neural response to an image patch required minimizing a nonlinear objective function via recurrent dynamics and (2) fitting relied on approximate inference methods that ignored uncertainty. Although subsequent work has developed several methods to overcome these obstacles, we propose a novel solution inspired by the variational autoencoder (VAE) framework. We introduce the sparse coding variational autoencoder (SVAE), which augments the sparse coding model with a probabilistic recognition model parameterized by a deep neural network. This recognition model provides a neurally plausible feedforward implementation for the mapping from image patches to neural activities and enables a principled method for fitting the sparse coding model to data via maximization of the evidence lower bound (ELBO). The SVAE differs from standard VAEs in three key respects: the latent representation is overcomplete (there are more latent dimensions than image pixels), the prior is sparse or heavy-tailed instead of gaussian, and the decoder network is a linear projection instead of a deep network. We fit the SVAE to natural image data under different assumed prior distributions and show that it obtains higher test performance than previous fitting methods. Finally, we examine the response properties of the recognition network and show that it captures important nonlinear properties of neurons in the early visual pathway.

Optical imaging of calcium signals in the brain has enabled researchers to observe the activity of hundreds-to-thousands of individual neurons simultaneously. Current methods predominantly use morphological information, typically focusing on expected shapes of cell bodies, to better identify neurons in the field-of-view. The explicit shape constraints limit the applicability of automated cell identification to other important imaging scales with more complex morphologies, e.g., dendritic or widefield imaging. Specifically, fluorescing components may be broken up, incompletely found, or merged in ways that do not accurately describe the underlying neural activity. Here we present Graph Filtered Temporal Dictionary (GraFT), a new approach that frames the problem of isolating independent fluorescing components as a dictionary learning problem. Specifically, we focus on the time-traces-the main quantity used in scientific discovery-and learn a time trace dictionary with the spatial maps acting as the presence coefficients encoding which pixels the time-traces are active in. Furthermore, we present a novel graph filtering model which redefines connectivity between pixels in terms of their shared temporal activity, rather than spatial proximity. This model greatly eases the ability of our method to handle data with complex non-local spatial structure. We demonstrate important properties of our method, such as robustness to morphology, simultaneously detecting different neuronal types, and implicitly inferring number of neurons, on both synthetic data and real data examples. Specifically, we demonstrate applications of our method to calcium imaging both at the dendritic, somatic, and widefield scales.

Neural-Matrix style, high-density electrode arrays for brain-machine interfaces (BMIs) and neuroscientific research require the use of multiplexing: Each recording channel can be routed to one of several electrode sites on the array. This capability allows the user to flexibly distribute recording channels to the locations where the most desirable neural signals can be resolved. For example, in the Neuropixel probe, 960 electrodes can be addressed by 384 recording channels. However, currently no adaptive methods exist to use recorded neural data to optimize/customize the electrode selections per recording context. Here, we present an algorithm called classification-based selection (CBS) that optimizes the joint electrode selections for all recording channels so as to maximize isolation quality of detected neurons. We show, in experiments using Neuropixels in non-human primates, that this algorithm yields a similar number of isolated neurons as would be obtained if all electrodes were recorded simultaneously. Neuron counts were 41-85% improved over previously published electrode selection strategies. The neurons isolated from electrodes selected by CBS were a 73% match, by spike timing, to the complete set of recordable neurons around the probe. The electrodes selected by CBS exhibited higher average per-recording-channel signal-to-noise ratio. CBS, and selection optimization in general, could play an important role in development of neurotechnologies for BMI, as signal bandwidth becomes an increasingly limiting factor. Code and experimental data have been made available.

Tracking algorithms such as the Kalman filter aim to improve inference performance by leveraging the temporal dynamics in streaming observations. However, the tracking regularizers are often based on the ℓp-norm which cannot account for important geometrical relationships between neighboring signal elements. We propose a practical approach to using the earth mover's distance (EMD) via the earth mover's distance dynamic filtering (EMD-DF) algorithm for causally tracking time-varying sparse signals when there is a natural geometry to the coefficient space that should be respected (e.g., meaningful ordering). Specifically, this letter presents a new Beckmann formulation that dramatically reduces computational complexity, as well as an evaluation of the performance and complexity of the proposed approach in imaging and frequency tracking applications with real and simulated neurophysiology data.

Calcium imaging has become a fundamental neural imaging technique, aiming to recover the individual activity of hundreds of neurons in a cortical region. Current methods (mostly matrix factorization) are aimed at detecting neurons in the field-of-view and then inferring the corresponding time-traces. In this paper, we reverse the modeling and instead aim to minimize the spatial inference, while focusing on finding the set of temporal traces present in the data. We reframe the problem in a dictionary learning setting, where the dictionary contains the time-traces and the sparse coefficient are spatial maps. We adapt dictionary learning to calcium imaging by introducing constraints on the norms and correlations of the time-traces, and incorporating a hierarchical spatial filtering model that correlates the time-trace usage over the field-of-view. We demonstrate on synthetic and real data that our solution has advantages regarding initialization, implicitly inferring number of neurons and simultaneously detecting different neuronal types.

Understanding why and how certain neural networks outperform others is key to guiding future development of network architectures and optimization methods. To this end, we introduce a novel visualization algorithm that reveals the internal geometry of such networks: Multislice PHATE (M-PHATE), the first method designed explicitly to visualize how a neural network's hidden representations of data evolve throughout the course of training. We demonstrate that our visualization provides intuitive, detailed summaries of the learning dynamics beyond simple global measures (i.e., validation loss and accuracy), without the need to access validation data. Furthermore, M-PHATE better captures both the dynamics and community structure of the hidden units as compared to visualization based on standard dimensionality reduction methods (e.g., ISOMAP, t-SNE). We demonstrate M-PHATE with two vignettes: continual learning and generalization. In the former, the M-PHATE visualizations display the mechanism of "catastrophic forgetting" which is a major challenge for learning in task-switching contexts. In the latter, our visualizations reveal how increased heterogeneity among hidden units correlates with improved generalization performance. An implementation of M-PHATE, along with scripts to reproduce the figures in this paper, is available at https://github.com/scottgigante/M-PHATE.

Tracking time-varying signals is an important task for practical systems working with large discretized domains. Under such settings, sparsity-based approaches improve tracking accuracy since typically few targets appear in the scene (i.e. few locations in the discretized space are occupied). Discretization introduces a unique challenge: the traditional lp -norm dynamic constraints produce significant errors when there is even a small spatial mismatch between the predicted and true state. To overcome this, we present a tracking algorithm leveraging concepts from optimal transport, namely utilizing the earth-movers distance (EMD) as a dynamic regularizer to the l1 -regularized inference problem (i.e., LASSO [1], or BPDN [2]). We extend the problem formulation to complex valued signals and modify the optimization program to reduce the computational burden. We demonstrate the efficacy of our approach in imaging and frequency tracking applications.

Neurons in many brain areas exhibit high trial-to-trial variability, with spike counts that are overdispersed relative to a Poisson distribution. Recent work (Goris, Movshon, & Simoncelli, 2014 ) has proposed to explain this variability in terms of a multiplicative interaction between a stochastic gain variable and a stimulus-dependent Poisson firing rate, which produces quadratic relationships between spike count mean and variance. Here we examine this quadratic assumption and propose a more flexible family of models that can account for a more diverse set of mean-variance relationships. Our model contains additive gaussian noise that is transformed nonlinearly to produce a Poisson spike rate. Different choices of the nonlinear function can give rise to qualitatively different mean-variance relationships, ranging from sublinear to linear to quadratic. Intriguingly, a rectified squaring nonlinearity produces a linear mean-variance function, corresponding to responses with a constant Fano factor. We describe a computationally efficient method for fitting this model to data and demonstrate that a majority of neurons in a V1 population are better described by a model with a nonquadratic relationship between mean and variance. Finally, we demonstrate a practical use of our model via an application to Bayesian adaptive stimulus selection in closed-loop neurophysiology experiments, which shows that accounting for overdispersion can lead to dramatic improvements in adaptive tuning curve estimation.

Two-photon (TP) calcium imaging is an important imaging modality in neuroscience, allowing for large-scale recording of neural activity in awake, behaving animals at behavior-relevant timescales. Interpretation of TP data requires the accurate extraction of temporal neural activity traces, which can be accomplished via manual or automated methods. In this work we seek to improve the accuracy of both manual and automated TP microscopy demixing methods by introducing a denoising algorithm based on a statistical model of TP data which includes spatial contiguity, sparse activity and Poisson observations. Our method leverages recent developments in stochastic filtering of structured signals based on Laplacian-scale mixture models (LSMs) to model the neural activity in TP data as a set of spatially correlated sparse variables. We apply our method on TP images taken from the visual cortex of an awake, behaving mouse, and demonstrate improved neural activity demixing over current pre-processing techniques.

Recurrent neural networks (RNNs) have drawn interest from machine learning researchers because of their effectiveness at preserving past inputs for time-varying data processing tasks. To understand the success and limitations of RNNs, it is critical that we advance our analysis of their fundamental memory properties. We focus on echo state networks (ESNs), which are RNNs with simple memoryless nodes and random connectivity. In most existing analyses, the short-term memory (STM) capacity results conclude that the ESN network size must scale linearly with the input size for unstructured inputs. The main contribution of this paper is to provide general results characterizing the STM capacity for linear ESNs with multidimensional input streams when the inputs have common low- dimensional structure: sparsity in a basis or significant statistical dependence between inputs. In both cases, we show that the number of nodes in the network must scale linearly with the information rate and poly-logarithmically with the input dimension. The analysis relies on advanced applications of random matrix theory and results in explicit non-asymptotic bounds on the recovery error. Taken together, this analysis provides a significant step forward in our understanding of the STM properties in RNNs.

Sparsity-based models have enabled significant advances in many image processing tasks. Hyperspectral imagery (HSI) in particular has benefited from these approaches due to the significant low-dimensional structure in both spatial and spectral dimensions. Specifically, previous work has shown that sparsity models can be used for spectral superresolution, where spectral signatures with HSI-level resolution are recovered from measurements with multispectral-level resolution (i.e., an order of magnitude fewer spectral bands). In this letter, we expand on those results by introducing a new inference approach known as reweighted l1spatial filtering (RWL1-SF). RWL1-SF incorporates a more sophisticated signal model that allows for variations in the SNR at each pixel as well as spatial dependences between neighboring pixels. The results demonstrate that the proposed approach leverages signal structure beyond simple sparsity to achieve significant improvements in spectral superresolution.

Cortical networks are hypothesized to rely on transient network activity to support short-term memory (STM). In this letter, we study the capacity of randomly connected recurrent linear networks for performing STM when the input signals are approximately sparse in some basis. We leverage results from compressed sensing to provide rigorous nonasymptotic recovery guarantees, quantifying the impact of the input sparsity level, the input sparsity basis, and the network characteristics on the system capacity. Our analysis demonstrates that network memory capacities can scale superlinearly with the number of nodes and in some situations can achieve STM capacities that are much larger than the network size. We provide perfect recovery guarantees for finite sequences and recovery bounds for infinite sequences. The latter analysis predicts that network STM systems may have an optimal recovery length that balances errors due to omission and recall mistakes. Furthermore, we show that the conditions yielding optimal STM capacity can be embodied in several network topologies, including networks with sparse or dense connectivities.

Accurate estimation of undersampled time-varying signals improves as stronger signal models provide more information to aid the estimator. In class Kalman filter-type algorithms, dynamic models of signal evolution are highly leveraged but there is little exploitation of structure within a signal at a given time. In contrast, standard sparse approximation schemes (e.g., L1 minimization) utilize strong structural models for a single signal, but do not admit obvious ways to incorporate dynamic models for data streams. In this work we introduce a causal estimation algorithm to estimate time-varying sparse signals. This algorithm is based on a hierarchical probabilistic model that uses re-weighted L1 minimization as its core computation, and propagates second order statistics through time similar to classic Kalman filtering. The resulting algorithm achieves very good performance, and appears to be particularly robust to errors in the dynamic signal model.

The sparse coding hypothesis has generated significant interest in the computational and theoretical neuroscience communities, but there remain open questions about the exact quantitative form of the sparsity penalty and the implementation of such a coding rule in neurally plausible architectures. The main contribution of this work is to show that a wide variety of sparsity-based probabilistic inference problems proposed in the signal processing and statistics literatures can be implemented exactly in the common network architecture known as the locally competitive algorithm (LCA). Among the cost functions we examine are approximate l(p) norms (0 ≤ p ≤ 2), modified l(p)-norms, block-l1 norms, and reweighted algorithms. Of particular interest is that we show significantly increased performance in reweighted l1 algorithms by inferring all parameters jointly in a dynamical system rather than using an iterative approach native to digital computational architectures.

This paper presents an algorithm for an ℓ1-regularized Kalman filter. Given observations of a discrete-time linear dynamical system with sparse errors in the state evolution, we estimate the state sequence by solving an optimization algorithm that balances fidelity to the measurements (measured by the standard ℓ2norm) against the sparsity of the innovations (measured using the ℓ1norm). We also derive an efficient algorithm for updating the estimate as the system evolves. This dynamic updating algorithm uses a homotopy scheme that tracks the solution as new measurements are slowly worked into the system and old measurements are slowly removed. The effective cost of adding new measurements is a number of low-rank updates to the solution of a linear system of equations that is roughly proportional to the joint sparsity of all the innovations in the time interval of interest.