CREIMBO: Cross-Regional Ensemble Interactions in Multi-view Brain Observations
Authors: Noga Mudrik
Ryan Ly
Oliver Ruebel
Adam Shabti Charles
Venue: ICLR 2025 Spotlight (2025)
Keywords: computational neuroscience
multi-regional brain interactions
sparsity
cross-session variability
dynamical systems modeling
neural dynamics
non-simultaneous neural recordings
Primary area: applications to neuroscience & cognitive science
TL;DR: We propose a framework for uncovering latent multi-regional neural sub-circuits by leveraging the richness and variability of multi-session neural data.

Abstract: 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.

Published as a conference paper at ICLR 2025
CREIMBO: CROSS-REGIONAL
ENSEMBLE
INTERACTIONS IN MULTI-VIEW BRAIN OBSERVATIONS
Noga Mudrik
Ryan Ly Scientific Data Division Lawrence Berkeley National Laboratory Berkeley, CA, USA. rlyłlbl.
Biomedical Engineering, Kavli NDI, CIS
The Johns Hopkins University Baltimore, MD, USA. nmudrikl k1@ jhu. edu
Oliver Rübel
Adam S. Charles
Scientific Data Division
Biomedical Engineering, Kavli NDI, CIS
Lawrence Berkeley National Laboratory Berkeley, CA. USA. oruebel @eooelloblggo
The Johns Hopkins University Baltimore, MD, USA. adamso jhu. edu
ABSTRACT
Modern recordings of neural activity provide diverse observations neurons across brain areas, behavioral conditions, and subjects; presenting exciting opportunity to reveal the fundamentals of brain-wide dynamics. Current analysis methods, however, often fail to fully hamness the richness of such data, they provide either uninterpretable representations (e.g., via deep networks) or oversimplify models (e.g., by assuming stationary dynamics analyzing each session independently). Here, instead of regarding asynchronous neural recordings that lack alignment in neural identity or brain areas 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 ceseemble-fuucttinn groups of neurons- such that the time-varys 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-simultancous observations, along with their non-stationary evolving interactions, with our new model, CREIMBO (Cross Pross-Regional Ensemble Interactions in Multi-view Brain Observations). CRE- IMBO 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 the use of shared underlying sub-circuit basis. Moreover, CREIMBO distinguishes session-specific computa- tions from global (session-invariang ones by identifying session covariates and variations in sub-circuit activations. We demonstrate CREIMBO's 's ability to re- cover true components synthetic data, and uncover meaningful brain dynamics in human high-density electrode recordings, including cross-subject neural mech- anisms as well as inter- vs. intra-region dynamical motifs. Furthermore, using mouse mmole-brain recordings, interactions that capture task and behavioral variables and meaningfully align with the biological importance of the brain areas they represent.
show CREIMBO's 's ability discover dynamical
INTRODUCTION
Identifying the interactions between and within brain areas is fundamental to advancing our under- standing of how the brain gives rise to behavior. Recent advances in neural recording technologies present an exciting opportunity to study brain-wide interactions by enabling simultancous recording of neural activity across many brain areas through multiple high-density electrodes. Such experiments,

---

Published as a conference paper at ICLR 2025
A observation space
N² Nº
nauron observed asaaion
challenges
time
(266
partial
view
learn
Iaaantt space
81
B
'reaa space
U...
ensembles
brain areas
compose ansamblee
Tibutons
shared underlying system
global dynamics prior
wwssion-spacific time changinnnnnngggggnnn dynamica' reweighting
Figure I: CREIMBO's Illustration. A: Real- Real-world multi-regional neural datasets consist of multiple (D) non-simultaneous recording sessions ({ (=1))=1) that cannot matched in terms of individual neuron identities, quantity, or cross-region distribution, resulting in neurons appearing only in certain sessions (4-pointed stars) or not at all (gray circles). This variability hinder our ability to draw unified conclusions of whole-brain neural computation, while analyzing sessions individually could lead to session-specifi biased results. This challenge further complicated by variability in session (or trial) duration, the non-stationarity of brain dynamics (even within session), and the presence of multiple trials within each session. B ~>: Instead of considering these asynchronous sessions frame them "views" onto the shared underlying brain system, thereby facilitating the learning of unified brain dynamics model. We assume that brain computations are mediated by multiple neural ensembles- * groups of same-region neurons with shared functionalites--Whir neural representations. However, the specific neurons and their membership degrees within each ensemble are unknown, and a neuron may belong to multiple ensembles with varying membership levels based on its diverse functionalities. C: CREIMBO leverages partial brain views to co-learn ensemble compositions per session ({ ((A") their temporal activity ((x" J (x 1), and cross-regional
challenge,
opportunity to obtain distinct, potentially complementary
interactions yield meaningful
interactions 2 Fdaze. It posits that these non-stationary interactions stem from session-shared dictionary prior of up to K global interactions, here termed sub-cireuits ( whose sparse time-changing decomposition = , shapes the overall ensemble dynamics at each time point and session d. The sub-circuit temporal coefficients enable distinguishing session-speciffe interactions from session-shared interactions in different time periods (pink curves in subplot C right). Under certain assumptions (Sec. D). the universal sub-circuit dictionary (f H ). shared across sessions, allows CREIMBO to auto-sort ensembles across sessions by functionality, aligning the j-th ensemble of session with the j-th ensemble of session d', despite differences in observed neurons.
repeated over many sessions with different settings, offer multiple asynchronous measurements of the brain system, with each session encompassing hundreds of distinet neurons across regions.
Individual recording sessions, however, provide singular perspective of the brain dynamics, as the reinsertion of recording devices results in the capture of different subsets of neurons in each session. Thus, fitting population-level models often reverts per-session analysis (i.e., fitting each session independently). which is inherently constrained by the inability to incorporate the full set of recorded activity across all recorded brain areas across all sessions. Moreover, analyzing the data dession-b--sesee may hinder our ability to distinguish between computations that are session- or subject-specific and those that are task-related or globally session-invariant Compounding the difficulty of merging data from different sessions is the non-linear and non-stationary nature of neural activity, along with noise and variability across trials-both both of which require additional model complexity account for. Moreover, to extract scientifically meaningful insights, models must be interpretable, with model parameters directly related to task variables or connections between recorded units. This gap between the opportunities offered by modern neural data and the limited
ressioniivariannoommuundiii
2

---

Published as a conference paper at ICLR 2025
capabilities of current methods necessitates new approaches to leverage the richness in modern brain data and discover the fundamental neural sub-cireuits governing brain activity.
One current approach used to combine multi-session data into unified model is deep learning. By training to predict all sessions' data, deep networks implicitly merge the datasets into the fit network weights. While these powerful models can combine information across sessions, their uninterpretable nature and complexity typically hinders their ability to reveal the fundamental building blocks of neural computation. Other approaches, including tensor factorization and its variations (e.g.. Chen et al. )2015); Harshman (1970g) may be limited in their ability to merge sessions of varying durations or model the interactions between the identified neural components.
Here, we re- frame cross-session recording variability not as a drawback, but as a valuable advantage, viewing this data as providing multiple, complementary "views" into a single, shared brain system (Fig. B > A). This approach allows comprehensive discovery of the underlying system by extracting joint information from the entire dataset, despite individual sessions differing in the identity and quantity of the neurons they cover. We hypothesize that the population dynamics at each session is driven by a universal set of global time-invariant sub-circuits shared across sessions. To effectively capture the non-linear and non-stationary dynamics of the brain, we enable these sub- circuits to present varying activation levels over time and across sessions. This approach allows us to identify evolving neural patterns by analyzing the joint activity of the sub-circuits each time point, captured through their time-varying decomposition. We assume that these sub-circuits govern the interactions between neural ensembles, which in turn drive the neural dynamics within a latent space shared across sessions, despite variations in recorded neuron identities across sessions. To account for this cross-session variability in neuron identities, we propose that the joint low-dimensional space in which the latent dynamics evolve transforms to each session's 's observation space via per-session mapping, while ensuring cross- session alignment of the neural groups' functionalities.
In this work, we lay out this new model, which we term CREIMBO Reoos--egiiona Ensemble Interactions Multi-view Brain Observations), and demonstrate its ability to capture latent dynamics in multi-view neural data. Specifically, our contributions include:
We discover multi-regional brain dynamics through leveraging the richness of modern neural datasets while ensuring interpretability.
We distinguish between intra- and inter-region interactions by identifying sparse, localized neural ensembles and disentangling co-active global circuits that govern their interactions. We accurately recover ground truth components in synthetic data and discover meaningful cross-regional brain interactions underlying human and mouse brain data.
2 RELATED WORK
Over the past decade, remarkable advances in neural data acquisition technologies have enabled the measurement of hundreds of neurons across brain regions, subjects, and experimental settings. This has spurred follow-up work to leverage this extensive data for novel neural discoveries, includ- ing merging non-simultaneous sessions and examining cross-regional brain interactions, all while accounting for the complex, non-linear, and non-stationary nature of brain activity.
Merging Non-simultaneous Neural Recording Sessions. Existing work for combining non- simultancous neuronal population recordings, here termed "sessions", include Turaga et al. (2013) who introduced the term "stitching" for this problem and modeled neural activity with a shared latent dynamical systems model and (Soudry et al. (013)) who proposed to stitch asynchronous recordings for inferring functional connectivity by assuming shared cross-session latent input. However, these methods are not designed to capture more than two simultaneously interacting neural sub-circuits across multiple populations. They also do not distinguish between within- and between-area interac- tions, nor do they specifically address temporal non-stationarity Bishop (2015) and [Bihhop & Yu (2014) proposed integrating the neural activity structure into inference and expanding stitching to study the communication between two neuronal populations from potentially different brain areas. Nonnenmacher et al. (2017) suggested extracting oow-dimensional dynamics from multiple sessions by learning temporal co-variations across neurons and time. Both of these methods, though, assume at least a small overlap in the identity of cross-session observed neurons. LFADS (Keshtkaran et al.) 2021: Pandarinath et al. 2018) enables stitching while coping with limited neural overlap by finding shared non-linear latent dynamical model. However, its latent dynamics may be difficult to interpret in relation to the specific neural sub-circuits that compose it, and because it is based on a variational autoencoder, tuning it can be challenging.
3

---

Published as a conference paper at ICLR 2025
Maintaining Interpretability in Non-Stationary Neural Dynamics Models. Real-world neural activity follows non-stationary nature due to changing environmental settings. Classical dynamical models often simplify neural dynamics to linear (c.g..jPCA) (hhurchCldl et al. (2012) or non-linear
but stationary (e.g.. (e.g.. [Gallego et al. (2000) processes. Models that capture non-stationarity often rely on "black box" deep learning models re.g..(Scnneiie
al. (2023): Pandarinath al. (2018): Zhu
iet al.] (202)). which are powerful but often present limited interpretability with respect to neural interactions. Other models are built on switching piece-wise linear models Linderman et al. (2016)
2017); Murphy (199), which do not enable the identification of multiple co-active dynamic processes. dLDS (decomposed Linear Dynamical Systems) (Mudrik et al.) [024a] CCen et al. (2024b) and its multi-step extension (Mudrik et al. (024c)), addresses this gap by generalizing the switching assumption to sparse time-changing decomposition of dynamical systems, thus capturing co-active sub-circuits. However, di dLDS is not inherently designed to distinguish multi-session within- and between-regional interactions with session-to-session variability in neuron identity.
Identification of Functional Neural Ensembles. A potential way to find interpretable loading matrix that links low to high dimensional space is through the identification of functional neural ensembles here referred to as groups of neurons with similar functionalities- that determine the axes in the latent space. For instance, the mDLAG model (Gokcen et al., 2024p used Automatic Relevance Determination (ARD) to promote population-wise sparsity patterns of each latent variable. Recently. Mudrik et al.] (2024b) presented an approach to identify sparse hidden building blocks underlying multi-way data while accounting for trial and condition variability. Their method includes a graph-driven regularization technique that groups co-active units into the same building blocks while pushing apart units with different functionalities. also enables capturing structural adjustments in block configurations across task conditions and per-trial variability in the temporal activity of the building blocks. However, their approach does not model the dynamic interactions between units nor provide a closed-form prior for the evolution of the functional groups' temporal traces. While other methods, e.g., clustering approaches like (Grossberger et al.) 2018) or Sparse Principal Component Analysis (SPCA) (Zou et al. (206)) enable the recognition of underlying sparse groups, they either do not allow a neuron to belong to multiple groups with varying degrees of membership or do not enable the ensembles to present structural adjustments to specific conditions.
Multi-regional Brain Interactions. With the ability to record from multiple brain areas and the
recognition of the importance of multi-area communication (Pesaran et al. (02T), recent works have explored communication within and between brain areas. In (Gokcen et al. 2022; 2024). the authors proposed an approach to identify time delays in cross-regional brain interactions and to distinguish indirect, concurrent, and bidirectional interactions in two or more populations. However, their model focuses on analyzing sessions individually and is not intended to uncover the full underlying set of transition matrices, but rather to recognize meaningful time delays that can imply interactions. Other models addressed multi-regional interactions through a communication subspace (Smeeoo et al., (2019) using dimensionality reduction, or by Generalized Linear Models (GLMs). with either Poisson Yates) cet al. 2017) or Gaussian (Yates et al./ (017) statistics to identify functional coupling between areas. These approaches, though, are not tailored to capture non-simultaneous cross-sessions variability in neuron identity and subject processing. Other recent methods (Karniol- TTmbbur et al. (0024);
Glaser et al. (2020); Li et al. (2024)) were tailored to model multi-regionil interactions. Particularly, they have proposed studying brain-wide dynamics using switching dynamical systems approaches, where the model switches between different states over time capture distinct periods. However, this approach cannot disentangle interactions among multiple simultaneously active sub-circuits, nor can it resolve interactions at different temporal resolutions occurring concurrently.
PROBLEM INTRODUCTION AND APPROACH
Let (Y be of estimated neural firing rates (e.g. from Neuropixels data) over D asynchronous sessions, indexed by d....... The neural recordings of each d-th session, yd E x capture T time observations of the activity of N' neurons from up to J dis- tinct brain areas, thus offering a partial view of the brain system. Note that sessions may refer to the same different subjects, can vary in duration or number of neurons (i.c., N' and ddd need not be equal to N' and Tdd for different sessions d, d'), and can differ in the subset of brain areas (out of J) they capture. Hence, these cross-session data matrices cannot be aligned into a single data array (e.g.. via a tensor), which hinders the direct application of existing analyses such as tensor factorization.
4

---

Published as a conference paper at ICLR 2025
We assume that these brain observations reflect the hidden activity X E RPXT" of p << N functional neural groups that evolve and interact over time. Each of these groups encompasses a sparse set of same-area neurons, where each neuron belong to more than one group with varying degrees of cross-group membership. The low-dimensional latent space in which the groups interact is shared across sessions and projected to each session's 's observed neurons via an unknown per-session projection g" R RNN with Y = (X
Since brain dynamics have been shown to be temporally changing (Petersen & Sporns 2015). we assume that the interactions between these functional groups follow a per-session non non-linear and non-stationary dynamics F IPP + Ip where X (x((1)) Finding such Fý, however, is nontrivial (as solving arg maxFe F. p(FF xe x4-1) each time point is intractable due having more unknowns than equations). Therefore, we follow dLDD (Mudrik et al. 2024a) and assume that these dynamics can be described by set [f of K global time- and session-invariant "sub-cireuits". Each sub-cireuit fx RRPP represents a core interaction between ensembles and is reused by the system at different time-points throughout its trajectory. Each of these latent sub-circuits may capture either global session-invariant brain interactions and/or session- and subject-specific interactions. Moreover, these sub-circuits may be active simultaneously or intermittently, yielding a model that can flexibly fit neural trajectories that differ between sessions while drawing upon single underlying mechanism. CREIMBO thus addresses the problem of leveraging the joint information from asynchronously collected observations { to 1) identify the K latent multi-regional sub-circuits (fi ) and their non -tattitionnay activation levels, and 2) recognize how they jointly and individually drive the interactions between neural ensembles, and ultimately define their trajectories in a shared-session latent space.
CREIMBO
Our approach simultaneously fits the global sub-circuits that underlie the inter-ensemble interactions across multiple sessions. It is predicated on three key assumptions: 1) ensembles of neurons, rather than individual neurons, form the basic units that interact in brain dynamies, 2) it is feasible to identify functionally analogous neuronal ensembles across sessions, and 3) the full repertoire of interactions between the neural ensembles is spanned by a global set of linear dynamical systems, with their linear combinations evolving over time to capture the non-stationary dynamics.
Let A be a sparse matrix where each column, Ad encodes the neuronal composition of the i-th ensemble in the d-th session, such that non-zero values represent the membership magnitudes of neurons to this i-th ensemble. Since we aim to capture multi-regional interactions, we design the ensemble matrix a A block diagonal matrix where the j-th block, a' E xps, contains neural ensembles from the j-th area only. We denote P2 and n, = No. as the number of ensembles and the number of observed neurons in that area for session d, respectively (p E/- - M)
Let := Id indicate the number of time-points recorded in each m-th trial of session d, such that the firing rate observations in that trial are captured by y == Y. E RRN"XT, and arise from the joint activities of these hidden neural ensembles, up to normally-distributed error e -N r). Namely, ax + E, where x := X E xgxT are the temporal trajectories of the ensembles. Please refer to Section G for Poisson statisties inference in low spiking rate regimes.
The evolution of the ensemble temporal trajectories ({X ) is determined via the latent interac- tions between them, which we model with a non-stationary linear dynamical system xt+1 F for t = 1.... T, where F; E RPxp, A key desired property of our model is to capture the overlapping activity of multiple sub-cireuits of interacting ensembles. We thus model the interactions (F,) through a time-varying sparse decomposition F, = =k= C where Cы = Ci., are the sub-circuits' time-varying coefficients in session d that capture the modulation of each circuit interactions.
m mm,a)
Fitting this model requires identifying the sub- circuits, their per-session coefficients, and the ensemble compositions. Hence, we employ alternating approach that iterates until convergence, between updating the circuits {
K and inferring, for each session d, 1) the ensemble compositions (a = ) 2) their trajectories (x = x²), and 3) the sub-circuit activations (Ck: = =CH)
Ensemble Update: We update a per-row (neuron) (a,)) Following the work of Mudrik al. (2024b), we infer the sparse ensemble structures with graph-driven re-weighted & regularization that groups together neurons with shared activity patterss, while pushing apart neurons that do not
5

---

Published as a conference paper at ICLR 2025
behave similarly. The neuron- neuron-similarity graph for each area j. h/ R is calculated based on a data-driven Gaussian kernel that measures the temporal similarity between neurons, such that 12 = expl [lv. v ), where on a hyperparameter. Specifically, we infer a row-wise, where the n-th row (an, E og TR?) indicates which ensembles neuron n belongs to and to what degree. Due to the block-diagonal structure imposed on a, in practice, we only need to update the membership magnitudes for ensembles within the region to which neuron n belongs (region j), while nullifying the membership to the ensembles of other areas. This is done via:
31
an aarg min yo. an,.: 55.11. + Mngland where An. 32 + |ān + 33 1 a.
(1)
on.
where 3y the activity of the ensembles of the j-th region and f1, B2. and B3 are scalars that control the effect of the graph on the regularization. The multiplication hn.. a. tests the correspondence be- tween 1) the temporal similarity of all neurons to neuron n (hn,.). and 2) the membership magnitudes of these neurons to the current estimate of the j-th ensemble (a..). This way, if hn, a..j j is high (i.e., the neighbors of neuron n on the graph are part of the current estimate of the j-th ensemble), Xn.j is small, resulting in lower sparsity regularization for including neuron in that ensemble, which ultimately promotes its inclusion. Alternatively, if O. is small, the regularization weight Xn.5 is large, and neuron will be less likely to be included in ensemble j.
The ensemble trajectories (æ;) and the sub-circuits coefficients ((cu,))e 1/ i) are updated iteratively for each time-point via the LASSO optimization
Latent State and Dynamic Coefficients Update:
Ge C. arg min 11 21,6r
ax, 112
1+41
Σ Co
Aellcelli
(2)
where dg E R is the sparsity regulatyrrgullaiizattin weight on the sub-circuit activations.
Dynamical System Update: The sub-circuits ((fx) which capture global core ensemble interactions, are assumed to be sparse (i.e., not fully interconnected), and are identified directly by:
Fall arg min |x= pill
Full T [
+ Xp
Σ In (k, #k2)
p(fki,
(3)
where Fall E RexpK is an horizontal concatenation of all {f)s, and xt E (Rx T is 16 the horizontal concatenation of all { where in practice this operation is taken on random subset of sessions D C D for computational complexity considerations. Above, L E Kp I" is the horizontal concatenation of all {{ (( 15 with &ेd E Kpx defined & x where 3 is the Kronecker product, and o is element- wise multiplication. The operator vec(-) flattens matrix a vector and Ap promotes sparsity within the sub-circuits. p(fk, fk.)) is de-correlation term with weight Xp used to ensure that distinct fs are not too similar. See Algorithm for a method summary and Section for complexity analysis.
as y - [c & [1]1xp) o [[]]Ix
5 EXPERIMENTS
CREIMBO Recovers Ground Truth Components in Synthetic Data: We first assess CREIMBO's ability to recover ground-truth components from synthetic data. We generated K = 3 synthetic sub- circuits represented by set of rotational matrices (Fig. ().. These sub-cireuits capture interactions among p = 3 ensembles from J = 3 regions (py 1 for all j 1,..., J). Each region consists of a random number of neurons drawn from a uniform distribution between 4 9. We generated the observation data {Yd)5_ 1==1 by simulating 5 sessions with different neuron identities and varying numbers of neurons per region (Fig. 5/ and D). CREIMBO was able to recover the ground truth components with high accuracy, as measured by correlation with the ground truth, across 312 random noisy initialization repeats with different parameters (Tab. []. The accurate recovery included all parameters-] 1) the reconstructed observations y (Fig. AA), 2) ensemble compositions A (Fig. 2B) and activity x (Fig. (C). 3) the sub-circuits F (Fig. 2D. and 4) circuit coefficients c (Fig. 2 When comparing CREIMBO's is results to components identified under diverse ablations (Sec. D). found that approaches that either separate the ensemble identification from the dynamic identification, analyze each session separately, or do not consider the localized structure of the ensembles- fail to

---

Published as a conference paper at ICLR 2025
A
D
30200073
6459
IND diCuRa
Without Dynamic or Regione (spame)
All Regions Onet Sparse)
Single
Vithout Dynamic
All Regions Isparsel
Rugians onoc
PCA Regions
Reglons
UAA CREMBO
Single
Ground
CREIBOO
Ablation Experiment
1-1
1-1
1-2
1
sesSions
witsou AlS
Aegist
RCA
Deckon
CHUMIO
00
sesacro
-44
62
-12
Comparison Others
SLDS grand SLDS Gierorad
FSLDS Odaa) FSLDS miag
FLLES- Cauns tra ISUS Coss (nis)
ISLS Foasso (pentred nq. SLDS Pooson Brabi
CAEIBO
Ground Trth
-logi Ri1
M
sessions
5ets/ons
525sions
Sessionz
Figure 2: A-E: Testing CREIMBO over 312 repeats with varying random initializations and random
seeds (Sec. B reveals high correlations between ground truth and fit components for all unknowns: Y.
A. x. F. c (subplots A-E respectively). F-G:Ablation Experiment: Comparison of the ensemble trajectory (ad subplot F) and the ensemble compositions (A4, subplot G) for the 1-st session as identified by CREIMBO vs. various ablations (details in Sec. D see Figure 7 for more sessions). H-J: Comparing CREIMBO ablation variants in terms of (- log (MSE)) between ground truth and identified dynamics, ensembles, and sub-circuits. "Single Session" refers to training CREIMBO on sessions individually. Details in Sec. [ K: Baseline Comparisons: The latent dynamics identified
by the baselines (SDLS, rSL.DS (Linderman et al., 2016), mp-rSI.S (Glaser et al.) (020)) for the first two sessions (details in Sec. D merged-session basclines provide a single trajectory). L-O: Comparison of the baselines' latent dynamics and ensemble matrices with the ground truth.
recover the ground-truth components (i.e., their resulting fit model parameters are less correlated with the ground-truth observations, Fig. - 7). We further evaluated CREIMBO's ability recover hidden components compared to other methods, including SLDS, rSLDS, and mp-rSLDS, and found that CREIMBO more accurately recovers the ground truth hidden dynamics and ensembles (Fig. K-O, Fig. g baseline details in Sec. D). Importantly, these existing models are not designed to identity multiple co-occurring processes or sparse ensembles, in contrast to CREIMBO, and they are not suitable for multi-sessions with varying neural identities, which limits the comparison from the outset.
Next, we tested CREIMBO on a richer synthetic example with = 15 sessions. Each session had a varying number of neurons (a random integer between 14 and 19 per region from uniform distribu- tion) and distinct sub-circuit coefficients (see Fig. II for the generated ground truth components and Fig. for regional distributions and masking). The sub-circuits were set rotational matrices, similar to the previous experiment. Across 204 random initializations with different parameters (see Sec. B we again found that CREIMBO consistently identifies the ground truth latent dynamics, ensemble compositions, and sub-circuits. Moreover, CREIMBO accurately reconstructed the observations with high correlation to the ground truth across all initializations (Fig. 1T with low relative error
7

---

Published as a conference paper at ICLR 2025
A
Uanaa snteri cortoe right
krft
G
f
6o
f1
6z
5,60
5.60
6/51
5160
$160
Figure 3: CREIMBO identifies cross-regional neural sub-cireuits underlying multi multi-session human brain recordings. A: Two exemplary observations (Y, right) compared to their reconstruction by CREIMBO (Y, left). B: Reconstruction performance in terms of R² and relative error across all sessions. C: The sub-circuit coefficients for 3 exemplary subjects. D: The standard deviation (std) of the coefficients over time across the two exemplary sessions. E: Exemplary identified sparse ensemble matrices (f (144)2)) for the observations from A. F: Relative std of the coefficients across four subjects reveals usage of similar sub-circuits. G: The identified ensembles' activity (Xd) for the two exemplary sessions. H: The identified sub-cireuits. Edge width indicates effect magnitude. Nodes represent ensembles, with colors indicating the ensemble area. Black arrows near fs and fu highlight within regiin interactions. I: Pairwise correlations between within-session sub-circuit coefficients (left, middle) and average within-session correlations across all sessions (size of '+' markers indicates std). J: Mean and std of fs values. K: Percentage of non-zero connections per fx-
(Fig. T... We further tested CREIMBO on more advanced data that included D = 40 sessions with a maximum of J = regions and py = 3 ensembles per region (Fig. [2]. We again found that CREIMBO recovers the components with high correlation with the ground-truth (Fig. [3].
CREIMBO Discovers Multi- Melti-Regional Dynamics in Electrophysiology Data: We then tested CREIMBO on human neural recordings from high-density electrode array provided by Kyzar ct al. (024b)) The data consists of neural activity from overall J = 10 brain areas with limited cross-session regional overlap (Fig. 14 15p, recorded while subjects performed a screening task (details in Kyzar et al. (2024b)). The data encompasses 21 subjects across D = 41 non-simultaneous sessions. For each subject (except Subject 19), the data offers 1) "Screening" session, and "Sternberg test" session (details in Kyzar et al. (0244)). We first converted the spike-sorted units to firing rate matrices by convolving the spike trains with 30ms Gaussian kernel (Fig. I0), and then tested CREIMBO on all sessions together with maximum of K = 12 sub-circuits and at most pj- = ensembles per region (full parameter list in Tab. P First, we note that CREIMBO can accurately reconstruct the data with high accuracy across all sessions (Fig. BA. B). Interestingly, the ensemble compositions identified by CREIMB (Fig. (3] present sparse patterns with one dense ensemble for most regions that we hypothesize capture the "background" mean field activity of that region (right-most column of diagonal blocks in Fig. 33.. Meanwhile, the other ensembles (the sparse ones) likely capture more specialized, nuanced functionality. The ensemble trajectories (X"d Fig. 3G) reveal some constantly active ensembles, implying their importance for neural processing.
The universal sub-circuits Hf 1) exhibit distinct localized motifs (Fig. 3D, i.e., most sub-circuits present clear trends of either sourcing or targeting the same area. This is in contrast to other approaches (Sec. D). that yield overlapping uninterpretable sub-circuits (Fig. 22 that emphasize only limited circuitry (Fig. 24). Interestingly, sub-circuin f10. f12 ended up almost empty (mainly zero-ish values, Fig. [9). This highlights CREIMBO's 's ability to automatically nullify redundant sub-circuits- i.e., form of model selection-vie the sparsity regularization. The number of active connections and
8

---

Published as a conference paper at ICLR 2025
average connection strength (Fig. 10J, K) show modest variability between sub-circuite highlighting the model's 's ability to identify distinct sub-cireuits with varying interaction styles, including those with varying numbers of connections, rather than being limited to sub-circuits with fixed sparsity. From a neuroscience perspective, this suggests that underlying neural sub-circuits can vary in density and the number of participating areas and ensembles. Interestingly, most identified interactions occur between distinct areas (i.e., inter-regional interactions) with either multiple ensembles from the same area affecting together ensembles of another area (e.g., fa). or multiple ensembles from diverse areas converging to the same target area (e.g., fs). These source or target ensemble groups vary across sub-circuits, meaning a single sub-circuit cannot capture the full repertoire of activity. This highlights the importance of CREIMBO's s abilitty to disentangle the combined activity of multiple circuits. We also observe, but to lesser degree, within-area interactions (e.g., black arrows in Ss, f11), which emphasize that within-area amplification or regulation exist in addition to cross-region computations.
When exploring the activation patterns of the sub-cireuits (Fig. P. 20), CREIMBO reveals cross- subject and cross-session variability. Notably, some sub-circuits are consistently used with high variability over time (e.g., Ss. in light-green), while others, such as f12, are not utilized at all in any exemplary session (Fig. BF, F). Moreover, when further exploring the sub-circuie that exhibit high activity variability across sessions (Fig. 3D.F), we observe that some are used within individual subjects in both tasks (e.g., f2 in subject 14, Fig. 3D), while others are used both across subjects and tasks (e.g.. f11, used by subjects 10 and 14 in both tasks). Within-session pairwise correlations between sub-cireuit coefficients (Fig. I01) reveals low correlations (<0.1) for most circuit pairs, which indicates that distinct sub-circuits differ in activity and reflect different cognitive processes.
Since ground-truth components are unavailable in real data, direct validation of such components is impossible. Hence, we sought a proxy to evaluate CREIMBO's ' performance on real-data by testing its consistency under increasing observation noise levels. Particularly, we re-ran CREIMBO on the human recordings after adding increasing variance i.i.d Gaussian noise (with K = 8, = 6, Fig. 23 CREIMBO remains robust to increasing noise, and experiences a phase transition at a specific noise level ( (Zuise ~ 0.2, Fig. A3A B, F, F. H, D) which align with the dictionary-leannig
Can
literature (Studer& & Baraniuk 2012). We found that the structure of correlations between within-area ensemble compositions (Fig. 23C, D) remains consistent under increasing noise levels, until around Onoise R 10, after which these correlations weaken as noise increases, This may imply that high noise in the data can obscure meaningful relationships between ensembles. Additionally, CREIMBO identifies that the sub-circuit coefficients remain robust until nnose 1.8, while exhibiting an increase in internal frequency as noise rises, potentially to account for noise that causes the observations to deviate from the system's ' dynamical rules (Fig. 23G 25). This suggests that rapid transitions in CREIMBO's circuit coefficients may indicate noise that disrupts the system's 's dynamics.
CREIMBO Discovers Regional Interactions Predictive of Task Variables: We tested CREIMBO's
ability to infer task-related variables from mice whole-brain Neuropixels multi-session data during a memory-guided movement task (data from Chen et al. (2024a), Fig. 26 27). CREIMBO identifies intra- and inter-area brain interactions via the sub-circuits (Sec. Fig. 2),.D), including cross- regional flows into or from key areas associated with skills needed for the task, such as memory (the hippocampus, serves as a source in fs and as target in fs ), planning (flows from frontal cortex in fs), and movement (flows from primary motor cortex fa). Within-area sub-circuits (Fig. (9)) further show within- secondary motor cortex in f2, fr), along with seff-activation//niiitooo of ensembles other sub- circuits (e.g.. f4). Moreover, the sub-circuit coefficients ((Cur) K k=n)) 1) capture task-related patterns across trials and sessions (Fig. [A.)). When training a regularized logistic regression model with the sub-circuit activations (fcR)) 1 as the only inputs, we were able predict various task variables, including outcome, early lick, and lick side, well above chance levels (Fig. [4C-,, p-values < x 10-11). Furthermore, sub-circuit activations across different task periods highlight how specifie multi-regional interactions capture different aspects of the task. For instance, towards the trial end (t3 window in Fig. 4 which includes the lick movement, C7 shows an increased importance (Fig 4B). Notably, e7 corresponds to the activity of the fr sub-circuit that captures flows to the secondary motor cortex 19g. Another example is cs (activity of fs that includes flows to hippocampus, Fig. 29g) which shows an increased feature importance in the first time window to- This carly-in-tas, importance of this hippocamppal-iiress
ensemble interactions (e.g., the basal ganglia in fi and the
aligns with the expectation that memory required at the beginning of this memory-guided task. CREIMBO thus demonstrates the ability to capture interpretable neural interactions and to predict task variables that reveal complex regional interactions (Sec. F).

---

Published as a conference paper at ICLR 2025
reearly
Sgnore
Ag1
O1
becores
D
E
Feature Importance
L0
H
wath and
MI
www.mua
Figure 4: Task-variabe prediction using CREIMBO's ' dynamic coefficients as input, based on coefficients from 19 sessions (40-60 trials each, Sec. F. A: Sub-circuit coefficients from an example session, shade marks within trials, bordered by starting and ending points. B: Sub-circuit coefficients for hit vs. miss trials. C: Accuracy score of predicting task variables based on coefficients vs. chance. D, F,G.H.IJJ: Confusion matrices of predicting varying task-variable based on sub-cireuit coefficients. E: importance of the different coefficients and time points for prediction.
DISCUSSION, LIMITATIONS, AND FUTURE WORK
Here, we introduced CREIMBO-a novel approach for uncovering multi-regional dynamics in neural data collected across multiple sessions. CREIMBO addresses the challenge of integrating non- simultaneous neural recordings by joint dynamical inference and sparsity regularization to capture the underlying neural sub-circuits governing brain activity. We further demonstrated the efficacy of CREIMBO through multiple synthetic and neural data, and found that CREIMBO recovered ground truth components and is robust to noise in identifying cross-regional motifs that span cross-session interactions. CREIMBO offers several advantages over existing methods. Chiefly, it identifies a cross-session shared latent space where non-stationary ensemble interactions are governed by time-varying decomposition of universal basis dynamics. By structuring these dynamics in terms of global sub-cireuits, CREIMBO allows the discovery of sub-cireuits meaningful for various cognitive processes, enabling the identification of variability in neural activity across subjects and sessions.
An important feature of CREIMBO is its ability to unify sessions with different neuron subsets through the universal dictionary of dynamical interactions prior that aligns ensembles in terms of functionality. This can enable the robust inference of ensemble activities even if different neurons are observed across sessions. This ability, however, depends on the extent to which the same sub-circuits are used across sessions, the distinctiveness of different {f Js, and the premise that sub-circuits containing interactions from or to ensembles that are entirely absent in certain sessions do also include ensembles that are observed in those same sessions (Sec. D). Another limitation is CREIMBO's reliance on linear projections from the latent space to per-session observation space, restricting flexibility. Extending non-linear projections offers potential for development, though introduces computational and interpretability challenges. Finally, CREIMBO uses dictionary learning, which computationally demanding, and future iterations will include parallel processing. Extending CREIMBO to additional applications, including non-neural data (e.g., immune-cell counts), is an exciting future step. Moreover, integrating CREIMBO with mDLAG (Gokcen et al., 2024), with the latter identifying optimal communication delays and the former leveraging these to further identify the underlying set of co-active LDSs, presents an exciting future work.
10

---

Published as a conference paper at ICLR 2025
ACKNOWLEDGMENTS
N.M. was supported by the Kavli NetroData Discovery Award of the Kavli Foundation, and as a Kavli Fellow by the Kavli Neuroscience Discovery Institute Johns Hopkins University. R.I., O.R, and N.M. were supported by the National Institute of Neurological Disorders and Stroke under Award Number 1224NS10005 (PI: O. Ruebel). A.S.C supported by NSF CAREER Award 2340338 and Johns Hopkins Bridge Grant.
CODE AND DATA SHARING
The human neural data used in this study Kyzar et al. (2024baa) is publicly available via DANDI
Archive at https: / /daadiiiii.mmmm rrchive org hrng//andd set veet/000 00070000499/0 240123 23. 1806 The mice neural data (Chen et al., 2024a 2023) is publicly available via DANDI Archive at https //dandiarchiv available on GitHub at this link
The code is
REFERENCES
Kathryne M Allen, Angeles Salles, Sangwook Park, Mounya Elhilali, and Cynthia F Moss. Ef- feet of background clutter neural discrimination in the bat auditory midbrain. Journal of neurophysiology, 126(5):172--1782, 2021.
Mikio Aoi and Jonathan W Pillow. Model-based targeted dimensionality reduction for neuronal population data. Advances in Neural Information Processing Systems, 31, 2018.
William E Bishop. Combining Neural Population Recordings: Theory and Application. PhD thesis,
Carnegie Mellon University, USA, 2015.
William E Bishop and Byron M Yu. Deterministic symmetric positive semidefinite matrix completion. Advances in Neural Information Processing Systems, 27, 2014.
Po- oo-Hauan (Cameron) Chen, Janice Chen, Yaara Yeshurun, Uri Hasson, James Haxby, and Peter Ramadge. A reduced-dimension fmri shared response model. In C. Cortes, N. Lawrence, D. Lee, M. Sugiyama, and R. Garnett (eds.), Advances in Neural Information Processing Systems, vol- ume 28. Curran Associates, Inc., 2015.
Susu Chen, Thinh Nguyen, Nuo Li, and Karel Svoboda. Mesoscale activity map dataset (ver- 00033/
sion 0.231012.2129, 2023. URL 231012. 2129,
[Data set].
Susu Chen, Yi Liu, Ziyue Aiden Wang, Jennifer Colonell, Liu D Liu, Han Hou, Nai-Wen Tien, Tim Wang, Timothy Harris, Shaul Druckmann, et al. Brain-wide neural activity underlying memory-guided movement. Cell, 187(3):67--611, 2024a.
Yenho Chen, Noga Mudrik, Kyle A Johnsen, Sankaraleengam Alagapan, Adam S Charles, and Christopher Rozell. Probabilistie decomposed linear dynamical systems for robust discovery of latent neural dynamics. arXiv preprint htt.:2408.116666 2024b.
Mark M Churchland, John P Cunningham, Matthew T Kaufman, Justin D Foster, Paul Nuyujukian, Stephen Ryu, and Krishna V Shenoy. Neural population dynamics during reaching. Nature, 487 (7405):51-56, 2012.
David F Crouse.
On implementing 2d rectangular assignment algorithms. IEEE Transactions on Aerospace and Electronic Systems, 52(9)1167--1666 2016.
Juan A Gallego, Matthew G Perich, Raeed H Chowdhury, Sara A Solla, and Lee E. Miller. Long-term stability of cortical population dynamics underlying consistent hehavior. Nature neuroscience, 23 (2):260-270, 2020.
Joshua Glaser, Matthew Whiteway. John P Cunningham, Liam Paninski, and Scott Linderman. Re- current switching dynamical systems models for multiple interacting neural populations. Advances in neural information processing systems, 33:14867-144 14878, 2020.
11

---

Published as a conference paper at ICLR 2025
Evren Gokcen, Anna Jasper, João D Semedo, Amin Zandvakili, Adam Kohn, Christian K Machens, and Byron M Yu. Disentangling the flow of signals between populations of neurons. Nature Computational Science, ((8):512-525, 2022.
Evren Gokcen, Anna Jasper, Alison Xu, Adam Kohn, Christian K Machens, and Byron M Yu. Uncovering motifs of concurrent signaling across multiple neuronal populations. Advances in Neural Information Processing Systems, 36, 2024.
Lukas Grossberger, Francesco P Battaglia, and Martin Vinck. Unsupervised clustering of temporal patterns high-dimensional neuronal ensembles using novel dissimilarity measure. PLoS computational biology, 14(7)xe1002283, 2018.
Z. T. Harmany, R. F. Marcia, and R. M. Willett. This is spiral-tap: Sparse poisson intensity reconstruction algorithms theory and practice. IEEE Transuctions on Image Processing, 21: 1084-1096, Mar 2012.
Richard A Harshman. Foundations of the parafac procedure: Models and conditions for an "explana- tory" multimodal factor analysis. UCLA Working Papers in Phonetics, 16(1):1-8: 1970.
Orren
David M. Zoltowski, E. Mika Diamanti, Lucas Pinto, Carlos D Brody, David W. Tank, and Jonathan W. Pillow. Modeling state-dependent communication between brain regions with switching nonlinear dynamical systems. In The Twelfth International Confer- ence on Learning Representations, 2024. URL https: http:::///ppeeeeeeeeee net/ fo forum?ic
WQWV7YBQWa
Mohammad Reza Keshtkaran, Andrew Sedler, Raeed H. Chowdhury, Raghav Tandon, Diya Basrai, Sarah L. Nguyen, Hansem Sohn, Mehrdad Jazayeri, Lee E. Miller, and Chethan Pandarinath. A large-scale neural network training framework for generalized estimation of single-trial population dynamics. bioRxiv, 2021. doi: 01.11000/22221001133226666
URL hhtp:::////// http::/wwwwww.xxi.ooe biorxiv
13....3..... 426570
Michael Kyzar, Jan Kaminski, Aneta Brzezicka, Chrystal M. Reed, Jeffrey M. Chung, Adam M. Mamelak, and Ueli Rutishauser. Dataset of somaan-single neuron activity during sterbberg working memory task (version 0.240123. 1806), 2024a. StLhttps://ddee http://i....rrrg//// org. 48324/
dandi. 0000690000099900.00
1806
Michael Kyzar, Jan Kamiński, Aneta Brzezicka, Chrystal M Reed, Jeffrey M Chung, Adam N Mamelak, and Ueli Rutishauser. Dataset human-single neuron activity during a sternberg working memory task. Scientific Data, 11(1):89, 2024b.
Chengrui Li, Soon Ho Kim, Chris Rodgers, Hannah Choi, and Anqi Wu. One-hot generalized linear model for switching brain state discovery. In The Twelfth International Conference on Learning Representations, 2024. httphhtss////peett
net. forum 2id0rrm0m022
Scott Linderman, Matthew Johnson, Andrew Miller, Ryan Adams, David Blei, and Liam Paninski. Bayesian learning and inference in recurrent switching linear dynamical systems. In Artificial Intelligence and Statistics, pp. 914-922. PMLR, 2017.
Scott Linderman, Benjamin Antin, David Zoltowski, and Joshua Glaser. SSM: Bayesian Learning and Inference for State Space Models (version 0.0.1, 10 2020. URL htt.::ttts//iiii gitt::b//gttu.....// com/
linde rwa.mmaaaaa al// ssm
Scott W Linderman, Andrew C Miller, Ryan P Adams, David M Blei, Liam Paninski, and Matthew J Johnson. Recurrent switching linear dynamical systems. arXiv preprint arXiv: 1610.08466. 2016.
Noga Mudrik, Yenho Chen, Eva Yezerets, Christopher Rozell, and Adam S Charles. Decomposed linear dynamical systems (dlds) for learning the latent components of neural dynamics. Journal of Machine Learning Research, 25(5)): 144, 2024a.
Noga Mudrik, Gal Mishne, and Adam Shabti Charles. Sibblings: Similarity driven building- block inference using graphs across states. In Forty-first International Conference on Machine Learning, 2024b.
12

---

Published as a conference paper at ICLR 2025
Noga Mudrik, Eva Yezerets, Yenho Chen, Christopher Rozell, and Adam Charles. Linoes: Lookahead inference of networked operators for continuous stability, arXiv preprint arXiv:2404. 18267, 2024c.
K. P. Murphy. Switching Kalman filters. Technical report, Technical Report, UC Berkeley, 1998.
Mark E Nelson and James M Bower. Brain maps and parallel computers. Trends in neurosciences, 13(10):403-408 1990.
Marcel Nonnenmacher, Srinivas C Turaga, and Jakob H Macke. Extracting low-dimensional dynamics from multiple large-scale neural population recordings by learning to predict correlations. Advances in neural information processing systems, 30, 2017.
Chethan Pandarinath, Daniel J. O'Shea, Jasmine Collins, Rafal Jozefowicz, Sergey D. Stavisky, Jonathan C. Kao, Eric M. Trautmann, Matthew T. Kaufman, Stephen I. Ryu, Leigh R. Hochberg, and et al. Inferring single-trial neural population dynamics using sequential auto-encoders. Nature Methods, 15:805-815. 815, Sep 2018. doi: (0.1038/441592-018
Bijan Pesaran, Maureen Hagan, Shaoyu Qiao, and Ryan Shewcraft. Multiregional communication and the channel modulation hypothesis. Current opinion in neurobiology, 66:250-257 2021.
Kaare Brandt Petersen, Michael Syskind Pedersen, et al. The matrix cookbook. Technical University of Denmark, 7(15):510, 2008
Steven E Petersen and Olaf Sporns. Brain networks and cognitive architectures. Neuron, 88(1): 207-219, 2015.
Steffen Schneider, Jin Hwa Lee, and Mackenzie Weygandt Mathis. Learnable latent embeddings for joint behavioural and neural analysis. Nature, May 2023. ISSN 1476-4687. doi: 10. 1038/ $41586-023-06031-6
htt.httpss//do..ortt
1 1000/0888000 4 41866 023 - 00031-601161
João Semedo, Amin Zandvakili, Christian K Machens, M Yu Byron, and Adam Kohn. Cortical areas interact through a communication subspace. Neuron, 102(())249-255)) 2019.
Max Simchowitz, Horia Mania, Stephen Tu, Michael Jordan, and Benjamin Recht. Learning without mixing: Towards sharp analysis of linear system identification. In Conference On Learning Theory, pp. 439-44944477 473. PMLR, 2018.
Daniel Soudry, Suraj Keshri, Patrick Stinson, Min-hwan Oh, Garud Iyengar, and Liam Paninski. A shotgun sampling solution for the common input problem in neural connectivity inference. arXiv preprint 10fi::1309.3222 2013.
Christoph Studer and Richard G Baraniuk. Dictionary learning from sparsely corrupted or compressed signals. In 2012 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), pp. 3341-3344. IEEE, 2012.
Aneel Tanwani, Hyungbo Shim, and Daniel Liberzon. Observability for switched linear systems: characterization and observer design. IEEE Transactions Automatic Control, 58(4):891-904, 2012.
Srini Turaga, Lars Buesing, Adam M Packer, Henry Dalgleish, Noah Pettit, Michael Hausser, and Jakob H Macke. Inferring neural population dynamics from multiple partial recordings of the same neural circuit. Advances in Neural Information Processing Systems, 26, 2013.
Jacob L Yates, Il Memming Park, Leor N Katz, Jonathan W Pillow, and Alexander Huk. Functional dissection of signal and noise in mt and lip during drccision-making Nature neuroscience, 20(9):
1285-1299, 2017.
Semir Zeki. A massively asynchronous, parallel brain. Philosophical Transactions of the Royal
Society B: Biological Sciences, 370( (700166804000 (668):20140174, 2015.
Feng Zhu, Harrison A Grier, Raghav Tandon, Changjia Cai, Anjali Agarwal, Andrea Giovannucci. Matthew T Kaufman, and Chethan Pandarinath. A deep learning framework for inference of single-trial neural population dynamics from calcium imaging with subframe temporal resolution. Nature neuroscience, 25(12); 01177724417444 2022.
Hui Zou, Trevor Hastie, and Robert Tibshirani. Sparse principal component analysis. Journal of computational and graphical statistics, 15(2):265-28) 2006.
13

---

Published as a conference paper at ICLR 2025
Appendix
A
A
A2
A
As
6.7
-075
-030
0U002
0.50
0.1
3.13
sessions
6.5.
5.51
059
02
Area
Area
Ares
Area session 5
session
session 2
session
session
0:99
session
session
session 3
session
session 5
Neuron #
206
30
Time
Time
Time
Time
Time
Time
Figure 5: A: Distribution of neurons across areas (areas a, b, and c) over the 5 different synthetic sessions. B: Ground truth sub-circuits are 3 rotational matrices, each captures a different rotational direction. C: Time-coefficients of synthetic sub-cireuits D: Synthetic ensembles ( E: Latent (ensemble) dynamics. Different trajectories in the latent space emerge as distinct cross- session time-changing decomposition of the ensemble sub-cireuits. F: Synthetic Ground Truth observations ((Y)))
A
B
C
p(X" [1 [1]
p(x" [2].xd [2D
p(x" [3], [3]
session
session
session
0.86
0:85
session
0.67
0.39
0.75
0.93
0.59
session
0.26
0.84
0.85
0.99
0.96
session
8.75
0:33
0.99
session 2
0.61
0.85
0.87
0.38
session
0.25
0.92
0.58
session 2
session
session
secsion
session
session 3
session 4
esssoon1
session? 2
session
sexsion
session
zesaion
ression
6995000
Figure 6: Clustermaps of correlations of cross-session latent dynamics. Each subplot captures one dimension of the latent dynamics.
14

---

Published as a conference paper at ICLR 2025
A
Sculio
WEPREA har
bperee
Tegler Spared
hogor
GEMDO
Misou Dyneric Prpicon tiuaree
Single
Reytins
Session
C
Session 3
E
Seesion
1.0
05
1-1
L-1
B
Without
Beg
Singl:
Regarn
hbv
spard
-55
-.00
1000
1000
1000
-015
0.75
Session
0.75
-0.0
14
023
1025
-66
Session
No. Nouror
-0.75
-0.50
2
42
0,25
14
-06
-06
Seasion
-64
0A
C4
-02
522
666
4075
-0A4
-G4
1000
Coent
Figure 7: Comparison between CREIMBO to the baselines reveal that CREIMBO outperforms other approaches (see baselines details in Sec. DD. A: Learned latent dynamics (X) for each of the sessions (rows) for the different methods (columns). B: Learnt ensemble compositions (A) for each of the 5 sessions (rows) for the different methods (columns). C-E: SSog(MSE)) and correlation between the ground truth and the the identified ensemble compositions (in C), latent dynamics (in D), and sub-cireuits (in E). F: Heatmaps of the identified sub-circuits for the different methods.
15

---

Published as a conference paper at ICLR 2025
621
Figure 8: Further comparisons of the learned latent dynamics these learned by SLDS, rSLDS, and multi-region ISLDS, with their variations (see Sec. D for details). Different columns represent different methods. Different rows capture different sessions.
B
Ags
115-
150-
1254
100-
50
2.4 R.R
Aree
Armes
Area
Aree
Figure 9: Distribution of areas in second synthetic dataset. A: Histogram of neural counts per area and session. Ensemble masks (Amask) across sessions. C: Ground truth sub-circuits.
Algorithm CREIMBO training
Input: p, Kh B1, B2, 33, des des As, Aps Aotss Batch Size
2: Initialize: { (ca) {zald=
Pre-caleulate: {hy]
4: repeat
5: for each session in a random batch of 17 do
6:
Update ensembles
7:
Update hidden dynamics and coefficients
8: Select a random batch of sessions from D Update networks {ff}
9:
10:
if stuck in local minimum then Perturb {fr ca}
12: until convergence
16

---

Published as a conference paper at ICLR 2025
reb 33 boy 6600 resl
0.
-02
-0
204
272
Time (iiens)
Figure 10: Rate estimation example for specific session using a 30 ms Gaussian kemnl..
A
NOTATIONS
Symbol
Description
D
Number of sessions (encompassing both same and different subjects)
Ma
Number of trials for each session d
Y
Neural recordings firing rate estimation for trial and session d
Na
Number of neurons observed in trial m of session d
Td
Number of time points observed in trial m of session d
Yd
Overall set of observations
m Jm.dd
Ad
The neuronal composition of each sparse ensemble in trial m and session d
m
A block
Block of A representing neural dynamics specific brain area
Total neurons area
Maximum number of ensembles discernible within area
p
Total number of ensembles across all areas
J
Total number of distinct brain areas
n
Total neurons involved across all trials or sessions
Et
State vector at time
F
Transition matrix at time f, representing ensemble interactions
K
Number of global interacting units
Ckt
Influence of the k-th global interaction time
fk
Basic linear system representing the k-th global interaction
T
Duration of the trial or session
B SYNTHETIC DATA EXPERIMENT DETAILS
We simulated the synthetic data represent the non-stationary dynamics of brain activity including multi-region interactions. Current assumptions about brain activity suggest Zeki (2015): Nelson & Bower 1990) functions through distributed, parallel processing across multiple regions, potentially involving co-active circuits that can process a range of variables (c.g.. task demands, feedback signals, and sensory information), simultaneously. This biological reality thus motivated our data generation approach, to reflects the multi pollts-processieie nature of neural activity, where different processes can interact and contribute to brain dynamics. We simulated D sessions, each sharing a common set of basis rotational dynamics, denoted as /=3 with p; = 3 regions of interest. Different regions were associated with different ensembles, and the dynamics for each session were modeled using sparse decomposition of the basis dynamics. For each session, the dynamics' coefficients, cf were generated over T = 500 time points. These coefficients were generated by initially giving them binary values (active inactive states), with random switches occurring every æ time points, where
17

---

Published as a conference paper at ICLR 2025
session
session
session
session
session
session
session
session
Neuron #
R
2
R-
49
Time
Time
Time
Time
Time
Time
Time
Time
session
session session session session session
session
0.25
0
-0.25
Time
Time
Time
Time
Time
Time
Time
Figure 11: Data and results for the second synthetic experiment (overall D 15 sessions, J=
regions and pj 1 ensembles for each region j = 1....J A: Ground truth synthetic observations. C: Ground truth ensemble
x was uniformly drawn between 10 and 30. The sparsity constraint was maintained at each time point, ensuring that only one dynamic active any time (Icc 1vt).
Then, to introduce smooth transitions between switching states, a Gaussian convolution with a time window of points (standard deviation 0.7) was applied the coefficients, promote smoother more gradual changes between active states. Different sessions maintained the same sequence of active dynamics but allowed for varying durations of each state, introducing temporal variability across sessions.
Each session d augmented with unique projection matrix Ad that captures the mapping between the latent space and the observations (neural) space in the respective regions. These projection matrices A consisted of 3 columns, each corresponding to one of the three regions. The number of neurons in each region was randomly chosen between 4 and 9, with the count varying independently across regions and sessions. This variability ensured that different sessions were not identical in their neural configurations. Latent dynamics for each session were generated using the recursive equation 2; Σk k= Cki C starting from initial state Eg = [1, -1, 1] CREIMBO was trained based on the observations only, which were defined for each session as dd dd + where is addition of i.dd Gaussian noise with staviird- of Cnoise 0.1
18

---

Published as a conference paper at ICLR 2025
100
count
sessions
fs
3m
6
Area
Area
Area
Area
Area
Area
Area
Area
Area
Area
Figure 12: Multi-ensemble synthetic data. Top: Distribution of areas across sessions. Middle: Ground truth sub-cireuits. Bottom: Masks for ensembles identification.
The full set of parameters for the synthetic experiment are available in Table [
C PARAMETERS FOR HUMAN MULTIIREGIONAL EXPERIMENT
The full set of parameters for the human data experiment are available in Table
D INFORMATION ABOUT THE BASELINES:
We distinguish between two types of baselines: (1) ablation experiments and (2) comparison to the performance of entirely different methods. It is important note that the latter cannot fully capture the unique capabilities of CREIMBO, particularly in identifying basis ensemble interactions shared across sessions through sparse decomposition. Therefore, we only compare certain aspects that can be extracted from other methods (e.g.. multi-regional SLDS, which requires the same number of neurons per session). However, these comparisons do not imply that other models outperform or underperform CREIMBO; they simply offer different insights, with varying strengths and limitations depending on the application.
Importantly, since the sub-cireuid } is invariant to its ordering, used SciPy's implementation of the 'linear tinear_sum_ s_ss_ssmm_aaaaaeeeeem_aaeeeesee_ assignment problem [Crouse (2016) to re-order the identified sub-circuits before comparing them with the ground truth, matching the fs values by minimizing the error. Ablation Experiments:
"Without Dynamic Prior":
Running CREIMBO while removing the prior over the dy- namics from the inference. Instead, inferring the temporal traces of the ensembles using regularized least-squares with the addition of smoothness, de-correlation, and Frobenius norm regularization terms the dynamics.
2. "All Regions (sparse)": Running CREIMBO without applying the multi-regional diagonal mask over the ensemble compositions (A), thus supporting the finding of non-localized ensembles.
19

---

Published as a conference paper at ICLR 2025
A
plr. 91 0.81 (rew)
ply 91 0.82
ply: 72= 0.82
ply =0 0.13
Plr, - 0.15
(real)
Imea)
-0.
Irea
18-
-0.5
21-
-02
04
20
02
-000
Onferrec)
(nferred)
(nfered]
Inferred)
Onferred)
10
-0.15
-0.00
-04
-0.15
-0.10
-0.3
00.05
00.00
-0.0
0..00
000
23 2
23 2 22
Time
Time
Time
Time
Time
PIA 0.90
AA. 00.91 (real)
PIA 0.91 (sea))
014, Á) 0.93
AAA,A 0.93
(neal)
(real)
(real)
10
0.0
00.
82
0.1
-0.2
69-
54-
63
63c
(infrred)
(nfereed)
(nferred)
(nfereed
Onferedi
-0.6
-0.8
-0.8
10
00.0
8.4
02
-02
48-
0.2
0.0
0.0
63-
Insembles
Ensembles
Ensembles
Ensembles
Ensembles
C
.20. -1 0.78
-0.7
21 - (real)
plx. R 0.82
p0x, 21
(real)
(real)
(real)
(rea))
Excember lettvrty
-10
-0.5
-1.0
0050
-00
-000
-0.0
00.00
--0
(nfeered)
(ndemed]
(nferrec)
(Iferee!
Ensemblr activty
-0.55
00.00
80
0025
0.00
-0.00
Figure 13: Multi-ensemble synthetic data results. Top: Recovering the observations Y (five exemplary sessions). Middle: Identified Ensembles (five exemplary sessions). Bottom: Identified latent dynamics (five exemplary sessions).
3. "All Regions (non-sparse)": Similar the previous one but without applying sparsity on the ensemble matrix.
4. "PCA All Regions": Running CREIMBO while replacing the dimensionality reduction step with PCA.
5. "Single Session #" Running CREIMBO single session (view) of the data rather than leveraging cross-session information.
Other methods:
1. "SLDS": We used the SSM Python package described by Linderman et al. at Lin-
derman et al. (2020). For the (non-recurrent) SLDS option, we used the "gaus- sian cattooo emission parameter, the "bbvi" as the fitting method, "variational "eariationallposteiiorr" was set to "mf", and the number of iterations ("num_iters") was set to 500. The of the fitting parameters
left default. We were inspired by
20

---

Published as a conference paper at ICLR 2025
Stmbberg
Dorsal onterior crogu cottor
Hypocamaus right
wwpppemeettrryy.o.mmm.. moter right
Ventral medea prctiontl cortee
620
404404
Neurons Count
subject ea
subject
subject
subjectt
subject 5
subject 6
subject
subject
subject 9
subject 10
subject
subject 12
subject
subjest 14
subject 15
subject 16
subject 17
subject 18
subject 19
subject
subject 21
Figure 14: Areas distribution in real world human neural data.
anyydata
amygdnie, right-
daruel anterior cingulate cortus, Jult-
223
denaal anterior.rrrroooo cingulate, corter, fight-
Hppecampen Jeft
1
Aispocampus, fight
0000
pra supplementery motor arsa left
15 1060
2
pre supplementary. moter area right-
ventral medial ventral_mediat www.prrrrtaalooomm cortex Jeft
ventral medial pww.rootall cortex right
Figure 15: Number of neurons per-area and per-session in the real-world human data taken from Kyzar ct al. (024b)).
this notebook
by Linder
man al. (2020) for our comparison. This includes:
1) "SLDS (trials)" which captures the training of SLDS with "bbvi" fitting method from across all trials (by stacking trials information vertically).
2) "SLDS (per-trial)" which refers to training SLDS individually per-trial.
2. "rSLDS": Similarly to SLDS, for the recurrent version (rSLDS, (Linderman et al. 2016)), we used the same SSM package provided by Linderman al. (2020). Here, used the recur- rent option for the transitions (i.e., the "transitions" parameter was set to"rrccrrent__mm _otl""), using "diagonal gaussian" dynamics and "gaussian_orthog" emissions. This includes: 3) "rSLDS (trials)" which captures training rSLDS with the "bbvi" fitting method from across all trials (by stacking trials information vertically).
4) "rSLDS (per-trial)" which refers to training iSLDS individually per-trial.
3.
Multi-Regional rSLDS, as described Glaser et al.] 2020). To run comparison, we were inspired by the colab-notebook for this "mp_rslds" method, provided by the SSM Python Package Linderman et al. (2020) (at /https://colaate www.errre.rcc..oooo.r.e wwwpgithub/.//coottoccc..o linde rman ab/ ssm/
rtt.pppllla..ottrs.............. ipynb). This in-
cludes:
5) ")mm__re porrrSLD-Gauud (trials)" which captures training Lim__rSLDD from across all trials
(by stacking trials information vertically) under Gaussian statistics.
6) "mp rSL.D-Gauss (per-trial)" which refers to training imp_rSLDS" individually per- trial under Gaussian statistics.
7) "mp_ POm__SLDSS -SiDS-Poisoon (trials)" which captures training iip_rSLDS" with the "bbvi" fitting method from across all trials (by stacking trials information vertically) under Poisson statistics.
8) 'mp_rst -SLDS-PPisse (per trial)" which refers to training "mp_rS "mp_iSLDS" individually per-trial under Poisson statistics.
21

---

Published as a conference paper at ICLR 2025
Basuani
Susuao
Sacumi
Dascon
Area
Area
Area
Area
Area
Area
Area
Aree
Aree
Aren
Figure 16: Multi-regional Masks for worldworld human data.
MOVERAN LOCKAL
coana sne catto (catroccial
MopiCaa ONLAA
Whs squaece dary aue art MUNTAA
Vort par
DeteCiCole
DIMAIC PYSUMOAYR)
Daten SroRi
PYSuNoddd
Coln
A003
Hbo
PrS.Molri.
PYS.MolARi
Figure 17: Five Exemplary Ensemble Matrices of the Real World Data, as Identified by CREIMBO
ASSUMPTIONS
The CREIMBO model is based core set of assumptions the nature of the data. These assumptions draw on both known properties of neural processes and general well-know statistical models used widely in data science. We categorize these assumptions into two types:
1. Priors over the underlying neural processes.
2. Priors over Observational Constraints in order for CREIMBO to properly infer the underlying system by leveraging multi-session information.
D.1 PRIORS OVER THE UNDERLYING NEURAL PROCESSES
1. The neural dynamics in each session d are assumed to lie on a low-dimensional manifold, embedded in a P << N'd diw-dimeniionl space that is defined by functional groups of neurons (referred to as "ensembles").
2. We assume that these functional ensembles consist of neurons with co-activation patterns, where neurons can belong to more than one, but only a few (a sparse number) of ensembles. with varying degrees of membership.
3. We further propose that the interactions between these ensembles drive the evolution of the latent manifold over time and are key encoding changes in conditions and behavior.
4. For CREIMBO be effective, we assume these interactions arise from the joint synchronous activity of multiple co- co-occurring processes, captured by a limited-size set of "basis--nn maasis-ensemble- interactions". The time-van varying, sparse decomposition of these sessis-nnemmbb--neeeeee weighted by their lome-local contributions in every time point, can adequately describe the manifold's s evolution over time.
22

---

Published as a conference paper at ICLR 2025
Observations
Reconstructions
Setsion
Neurors
session
secsion
session
sossion
sess.on
session
session
secsion
Neuruns
session
Session
Session
Seasion
Neurons
session 13
session
session 15
session
Neurors
Neurors
session 17
Seision
Gession 10
Setsion 20
scation 17
Scssion 18
session
Scssion
Neurons
9e5sion
Sossion
secsion
Neurons
session 25
Session
Gecdion 77
Session 28
sension 25
seesion
secston 27
session 28
Nourons
Neurams
session 29
session
session 31
session 32
session
Newrons
33
35
Neuruns
sexdion 37
seaxion 39
seszion 40
scation 37
scssion 38
scasion 39
sossion
Neurons
Neurons
Neurons
Figure 18: Observations vs. CREIMBO's reconstructions for the real-world human data experiment.
5. We assume that different dynamics basis elements capture distinct processes behaviors. Some processes are required to be globally related to the cognitive task, while others capture session- subject-specifie processes (see Sec. D.2
23

---

Published as a conference paper at ICLR 2025
610
f22
B
6
p0..0
Figure 19: The sub-circuits underlying human data identified by CREMMBO. A: The identified 12 sub-circuits. B: Pairwise correlations between the identified sub-circuits.
Table I: Parameter values for synthetic experiments
Parameter
Value (range)
Explanation
same_c
False
Indicates if to shared sub-cireuits coefficients
across sessions.
step
0.1 0.5
Gradient descent step size for updating f
0.99 0.999
Decay over iterations of gradient dessent step size.
max_error
Ie-09
The error threshold to stop training the model.
max iter
500
Threshold on the maximum number of iterations. Indicates the inclusion of projection to latent space.
include D
True
step_D D
10-5 10-2
Range of Gradient Descent Step size for A.
seed
0-4000
Random seed to choose from
normalize_cre
True False
Indicates if to normalize the sub-circuits by
start sparse_
False True
Whether initialize to be sparse.
sparse_
True
Whether to apply sparsity the sub-circuits
num gradient_stee steps
-3
Number of gradient steps in an iteration.
add avg
True False
Whether add moving average an iteration.
sparsity _o___oo_________ f max
40-60
Percentile of sparsity applied on each Indicates if to take multiple steps of gradient descent.
take_ take_multiple_dd
False, True
D D__raphh_drin _drven
True
Indicates that D) is inferred with graph-driven way. Indicates the inference of æ and e occurred simultaneously.
infer x cifer_x_c_ttgeeel ___together
False, True
include_mace
True
Indicates the inclusion of a mask on A Indicates the normalization columns
norm_D_com
True
lambda_1
0.3
Regularization weight on updating A
step D d_ccc
0.99 - 0,9999
Decay rate for updating A parameters
num regions
Number of brain regions.
lambda_x
0.1.1
Sparsity term the latent dynamics.
latent dim
3
Latent dimension.
noise level
0.2-0.6
Standard deviation noise added to the data.
num_subdyns
3
Number of sub-cireuits.
lasso_so solvessolverrre
'spgll'
Pylops solver used for & regularization.
6. Each of these ensemble interactions may capture between- interactions, within-area interactions, or both.
D.2 PRIORS OVER OBSERVATIONAL CONSTRAINTS
To effectively learn a unified representation by leveraging information across sessions, we make the following assumptions on the statistics of the data.
24

---

Published as a conference paper at ICLR 2025
Screening Ceofficients (c)
Sternberg
Ceofficients (c)
subject
subject
3
sutjject
subject
subject
subject
subject
9
AANAAAAA A
subjact
MOLODA
subject
subject 12
subject 13
AM A MA
subject 14
in MAAAMAA
subject 15
subject 16
subject
subject
subject
subjact
subject
6.0
6.0
Time (s)
Time (s)
Figure 20: Sub-circuits Coefficients separated by task type ("Screening" vs "Sternberg") and subject.
1. Our ability to infer the underlying latent state (see Fig. 13) is contingent upon the overall observability the dynamical system. In traditional linear systems, the observability matrix guarantee that any state is visible in finite time horizon output. Non-stationary systems, as in CREIMBO, do not have succinct guarantees, the closest which come in the form of observability conditions switched linear systems (Tanwani et al. 2012
25

---

Published as a conference paper at ICLR 2025
5.bRI14 Subaz110
M
Apbr
11
MAA000
R
P
p
Rt
Marton
Time 6s)
Time Dcl
Time 150
Timne
Time 153
Timae (al
Time Eil
Time 160
Time Dal
Time (S)
Time
Tima Oil
Figure 21: Comparing the use of the same coefficient under 'Sternberg' vs 'Screening' within four exemplary subject, reveals usage of similar sub-circuits within subjects.
CrEIMBO
Session (SPARSE)
Session (SPARSE)
Dynamic Prior
plf 00
phh 5)
ptf. 9
6.6.6 6. AAMifo pth. f]
KENLA C.. LURRAA plf, f))
5664647 LAC.Ai pl6, f))
6.6
plh, f
Session 3 (SPARSE)
Session (SPARSE)
All Regions Single Session
All Recions
plr 0
pl.. f))
plf 60
pt..
plf.
plf. f)
plf,
Correlations between the brain sub-circuits identified by CREIMBO compared to the other methods. These results highlight the CREIMBO discovers more distinct sub-circuits than these idetified by the other methods.
Figure 22:
As the heart of the observability condition is that the dynamics "well mixes" the state such that any state element in the null space of the readout matrix will eventually be rotated into the its span, and thus visible. We thus assume that the spectral radius of each dynamical system (a measure of mixing in linear systems (Simchowitz et al. 2018)) in the basis is close one. This assumption is loose, due the lack of theoretical guarantees on general non-stationary systems, and further analysis should identify tighter, and less stringent, assumptions.
2. Another prior we want to consider is which fs capture ensembles unobserved in certain sessions, allowing to assert that CREIMBO can recover the activity of these unobserved ensembles during those sessions. To achieve this, we must ensure that in the sessions where these ensembles unobserved, at least some fs representing them also capture the activity of observed ensembles. For example, if certain missing ensembles are represented solely by fs that do not capture any observed ensembles, of these unobserved ensembles in the sessions where they are absent, because there are no observed ensembles linked through the dynamic prior.
will not be able to infer the activity
3. To ensure uniqueness of the sparse decomposition the dynamical systems model at each session and time point, assume that the effective spark of the dynamics basis is large (S* 2 2S). Essentially, there is no fx that can be linearly composed of 2s* other dynamical systems. This ensures that there cannot exist two (or more) equivalent
26

---

Published as a conference paper at ICLR 2025
Sultwes
Sudjjct (Shembbrr))
Subject (Scheningg
Sutjjctk Aacmemron
Subject (Stemberr
Sulowct
1oo0
(chembarol
20
00000
anggehee_ssn amppdaa right
ornel www.rreerrcoooommmomooo www.erttrrrggettt.ooo... corten, aft anterio of di dirulate cortex righ
pocampi Mtt poscampus right
Hr www.seeaataar...ooooo meto jra
ventodl Lrteeertttteeoomo
aeciemertary nater area fight
ventral medil prefirtal cortes right
Chatce
blodee brsn 20010
Enser neeentees
within atta
Subject 10
Subject 10
Subject 11
Subject 11
Subject
Sugject 12
Subject
Subject 13
Subject
ts
100
(oog
530
-0.0
H
pleeage
-0.
Figure 23: CREIMBO'S Robustness to Increasing Random Normal Noise Levels. A & B: Cor-
relations of identified ensemble compositions (Ad) under increasing noise, for individual sessions (A) and all sessions combined (B). Robustness decreases only when Rasis = 3. C & D: Correlations of ensemble compositions (from data concatenated across all conditions) displayed in a heatmap, with areas represented by rows of blocks (vertical) and noise levels by columns (horizontal). Block (i.j) captures the correlation of ensembles from area at the j j-th noise level. A sharp drop in within-ares correlations (D) occurs at noisc
Ganis
2.. & F: Correlations of ensembles'
trajectories (fxd) ) for each condition (E) or combined (mean/median) (F) show robustness until Gadia 0.27. G: Time-varying interaction coefficients for the 1st condition "(sub. 10, Screening)"
dan
under increasing noise reveal similar but more frequent changes, with a sharp frequency change at moise 5 Cume = 0.33). H & I: Correlations between corresponding (cr. K "(sub. 10, Screening)" over increasing noise (H), and compared to the coefficients identified from the original data (H).
da
decompositions ΣK_ fECK and EK facla for distinct coefficients that both compose the same dynamical matrix.
4. Lastly, to ensure accurate cross-session alignment, we propose that each dynamical system projects uniquely onto the obseverved subset of ensembles in each session. This assumption prevents two dynamical systems (e.g.. "ensemble - ensemble 2 > ensemble 3" and "ensemble 1 + ensemble 4 - ensemble 3"") from being indistinguishable from the data available that session (e.g., session with only ensembles and 3 recorded).
E COMPUTATIONAL COMPLEXITY
CREIMBO's learning process involves learning both the ensemble compositions per session ((A))))) the ensembles' temporal activity ((x4) and its underlying temporal evolution, which requires identifying the global (session-invariant) dynamic operators ((f) ( and their
per-session temporal coefficients ([)))
27

---

Published as a conference paper at ICLR 2025
Table 2: Parameter values for real-world human experiment
Parameter
Value False
Explanation
same_e
Indicates if sub-circuit coefficients are shared
across sessions.
step f
0.2
Step size for gradient descent in updating f.
GD deccy
0.992
Gradient descent step-size decay rate over iterations.
max error
1
10-9
Maximum allowable error.
max_iler
100
Maximum number of iterations.
seed
Random seed.
normalize cig
True
Normalize cigenvalues.
Spart_s
False
Start with sparse c.
include_D
True
Whether include A.
step_D D
0.0001
Step size for updating A.
num gradient steps
4
Number of gradient steps.
add avg
True
Add average.
sparsity_on Lacily__o____ _at____anaaaaa_ max
45
Maximum sparsity on f.
like__mltiple__01 gd
False
Take multiple gradient descents.
D_ D_graph_driin _ddiven
True
Graph-driven A.
infer x_ __toeehher
False
Infer and together.
include mask
True
Include mask.
norm D DD_cols
False
Normalize columns of A.
lambda D
0.3
Regularization parameter for A.
step_ Susp_Dddecayy
0.9999
Decay rate for step size of A. t regularization for A.
11_D
2.74
E COMPLEXITY OF ENSEMBLE-MATRI UPDATE
The loading matrix update relies on 4 main computational steps:
1) Channel Graph Construction:
This operation, performed once for all N channels of every state d 1..... D, generates a channel graph H N for each state d E [1, D] by concatenating within-state trials 1... Ma horizontally, resulting in N EMr, T matrix. For simplicity, let Ma Td. The computational complexity of calculating the pairwise similarities of this
T MM concatenated matrix for all D states is thus O (DÌ'N(N- - 1)).
2) The k-threshold step: involves keeping only the & largest values in each row while setting the other values to zero zero-the complexity will be O (Tlog k) per row for total computational complexity of O (DNTlog k) for N rows and D states.
3) State Graph Construction:
This is a one-time operation that involves calculating the pairwise similarities between each pair of states. For simplicity, if we assume the case of user-defined scalar labels, and as in this case there are D states (and accordingly D labels), the computation includes D(D 1) pairwise distances for (D2).
4) Ensemble Inference (Eq. equation 1: This iterative step involves per-channel re- ce-weighted E optimization. If the computational complexity of a weighted E is denoted C, then the computational complexity of the re- Weighted & Graph Filtering is N LNNk, where N is the number of channels, is the number of iterations for the RWLF procedure, and is the number of nearest neighbors in the graph. For the last term in Eq. equation 1 there are p² multiplicative operations involving the vector v/ and the difference in ensembles, arising from the lz norm. Additionally, there is an additional multiplication step involving Par. For each state d, this caleulation repeats itself D 1 times (for all d' f d). This process is carried out for every d = 1... D. In total, these multiplicative operations sum up to (p² 11) D(D 1), resulting in computational complexity of O pD2
28

---

Published as a conference paper at ICLR 2025
E.2 COMPLEXITY OF INFERRING THE ENSEMBLES ACTIVITY AND INTERACTIONS COEFFICIENTS (EQ. EQUATION 2
The inference of of for each 1 1 amd each condition d 1.. D, involves: 24
arg mille 2: 11 - az 1 + @ F & + Aelleillz where FX E Rex is the hor- izontal concatenation of (f
This simplifies Lasso problem of the form
1
min 2 ||5 M0l1 s.t. ||0|| S 7,
with dimensions & (N + 1, M: (N + (p + K), and 0 (p + K) 1, where
{
yt 0pxi
E <R(N+iix1
x Cr
E R(++K)11
and
M
A
ON*K FX
E E
I
This can be solved with SPGLI with up to num_iters iterations, resulting in a computational com- plexity of
O (N + (p + K)).
per session and time point and overall (across all sessions and time points):
D
( iters (N+ pp K)
T
E.3 COMPLEXITY OF INFERRING THE DYNAMICS' DICTIONARY
The inference of the basis-interactions dictionary {fa] (via Eq. equation involves solving
Fan arg min | Fall
IC
where + E Rpx Edchaard and ca E pKpxE Id For simplicity, assuming that Ap 0, The problem simplifies to a LASSO formulation, which is solved using SPGL11 for Fulll with an overall computational complexity per iteration of O iter Kp²r). where num iter represents the number of iterations required for convergence.
F APPLICATION TO MULTI-REGIONAL CROSS- MICE NEURAL ACTIVITY
We further tested CREIMBO's ability to infer meaningful task variables underlying whole-brain multi-session data from mice Chen et al. (2024a) while performing gvemory-guided movement task. The measurements were obtained using multicleectrdee extracellular clectrophysiology. with the overall dataset stored in DANDI Archive under the NWB format Chen et al. (2023).
Particularly, the data include electrophysiology recordings of neural activity from multiple brain areas per session (Figs. 26 27), across overall 175 sessions, with numerous trials in each. The task the mice performed during the recordings involved selecting one of two "lick ports" based on auditory cues, with the left right port leading reward depending whether they heard high low
29

---

Published as a conference paper at ICLR 2025
Table 3: Parameter table for the CREIMBO Mice Mesoscale.
Parameter
Value
max max_error
1e-09
seed
0
normalize Fis
True
sparse_f
True
sparsity_on__ maxx
40
increase in jn_sparsity_______ f
1.4318630048860998
norm_D norm_D_colls
True
D_ D_graph_ppaams
('with_kNN': :NN': True, 'with_norm': True, 'k': 15)
lambda_D
0.4392235894758376
pydate_type___
spgll
latent aitent_dyns___ sii__ddns__iiii initiiiiiiiaiinnionnnn
random
D with lasso
True
11
2.7440675196366233
params D num_rege _errreioon
("update_c_type' 'spgl1")
10
reg_ rrg_tpe_onn__
spgll
lambda_x
1.029122116748077
latent_dis dinntt_dmm__
30
num ensembles per region
30
number of sub-dynamics
10
seed_f
0
digma_mix_
0.1
sparse f params
('axis': '1' 'percentO": 20)
tone. Neural activity was recorded across several brain regions in each session, with some regions overlapping between sessions, though no two sessions had identical coverage (refer to Fig. 1., A-B. Chehn et al. (2024a)).
We first loaded data from all sessions and selected areas for CREIMBO that featured the highest number of sessions including them. We identified the parent area of each sub-area labeled in the dataset and removed neurons with unlocated or NaN areas. Inspired by our exploratory analysis of parent-area distributions across sessions and neurons (7B.C) we selected the 10 parent areas with the highest number of neurons across sessions We then selected 20 random sessions. One of these 20 sessions contained only (Fig. 2
neurons, leaving 19 sessions for this analysis
We pre-processed the data using 30 ms Gaussian kernel for firing rate estimation across all the sessions used (see FF.. We then selected only the neurons from the 10 parent areas and limited the analysis to the first 60 trials. Some sessions exhibited no activity or disconnection certain trials, meaning some sessions included fewer than 60 trials (but no fewer than 40).
Next, we CREIMBO with 10 sub-eircuits (K 10) and Pj 3 ensembles per area. The full list of parameters shown in 3 and an example exploration of a random session is shown in (Fig 28 A-C show the firing rate, trial start-end times, and distribution of trial durations, respectively).
We identified ensembles per region (exemplary in A9 with sub-circuits that capture both between- area (Fig 29C) and within-area interactions (Fig 29 showing diverse motifs with varying degrees
Ussed Parent Areas: Thalamus, Orbital Cortex, Secondary Motor Cortex, Basal Ganglia, Olfactory Cortex, Brain Stem, Somatosensory Cortex, Frontal Cortex, Hippocampus, Primary Motor Cortex
"Sessions DANDI names: 72255681--000--222444444--aaaaa-
"bb1039554-9991446611
c6f91a1e9313',
9033--772addfAccd: 44750
'22791d80-266 http:///2....4444000000000000000
0b4e-2233883bh14422
30

---

Published as a conference paper at ICLR 2025
of cross-cireuit correlations (Fig, 29B). Some cross-area interactions reveal effects originating from specific regions involved in the task. For example. fi demonstrates hippocampal effects, which align with the task's 's memory demands, while, e.g., f3 shows interactions involving the frontal cortex, including follow-up activations, corresponding to the task's need for planning and execution. Other sub-circuits reveal inputs converging on specific areas. For instance, fs shows motor cortex inputs, relevant for moving the tongue during the lick response, and fz reflects inputs into the secondary motor cortex, which is involved in coordinating motor actions.
When extracting only within-area interactions, we observe unique motifs as well. For example, fi shows further internal processing between basal ganglia ensembles, which may reflect the internal coordination and integration of motor planning and execution signals. fr and fz show strong effects of one ensemble on another within the secondary motor cortex, which could represent the modulation and refinement of motor commands for more precise movements. fs demonstrates brainstem effects, with two ensembles influencing two other ensembles, potentially supporting basic motor functions such muscle activation and coordination. f4 shows multiple self-activations of ensembles, which may capture self-regulatory processes like feedback inhibition or facilitation that modulate motor output and prevent activation..
We further examined CREIMBO's dynamics coefficients (ct) and identified task-related patterns with similar coefficient profiles across trials for the same task variable (Figl4) 4A for an example session, and B for hit vs. miss outcomes). Similar differences were observed across other task variables. To quantify the predictive power of the dynamic coefficients, we used them input for training simple one-vs-rest logistic regression classifier. Each trial was split into four equal-duration time windows, and the 10 dynamic coefficients were averaged within each window, resulting in 40 features (4 windows x 10 coefficients). These features were then used to predict task variables, including outcome (hit/miss/ignore), early lick (whether lick was early), lick side instructed, and lick side performed in practice. Additionally, we tested the model to predict 2-3 variables simultaneously, which was more challenging due to the need to distinguish between multiple options. The resulting accuracy levels (Fig 4C, black stars) were significantly higher than chance levels (gray bars). We further analyzed the confusion matrix for each variable's s prediction and applied the x² test to assess whether the distribution of predictions significantly deviated from the expected distribution under chance. We found p-values well below 1 x 10-10 (see subtitles of Fig 4 D. F, G, H, I, J), strongly supporting CREIMBO's predictive power.
Furthermore, we assessed the importance of subcircuits' activations across the four time windows 29 with each block of rows representing one subcircuit and the blocks corresponding to different time windows. CREIMBO revealed patterns where certain time windows and dynamic coefficients were pivotal for different task variables. For instance, C7 at 13 was critical, while other coefficients, like at t3. were specific to encoding the lick side. This specificity aligns with the fact that this time point likely reflects the final stages of processing during the task, including aspects of learning. When we examined cy (Fig 29C), we observed multiple inputs into the secondary motor cortex, related to the execution of the lick towards the trial's end. highlighting CREIMBO's ability to pinpoint region-specific interactions. Additionally,
showed increased importance in the initial window 40 contributing across multiple variables and capturing inputs into the hippocampus- an area integral to memory processing-une underscornnggggggggggggo its role in early trial stages. These findings further validate CREIMBO's s ability to capture and predict task-related neural dynamics, offering insights into circuit-level brain activity across time.
G INFERENCE ASSUMING POISSON STATISTICS
CREIMBO, as outlined in the main text, assumes Caussian-distributed i.i.d noise. This assumption is reflected in the inference procedure, which minimizes the Eg and Frobenius norms (e.g., equa-
tion [ equation D Particularly, these norms arise from minimizing the negative log- likelihiddo of
yd - -IA" 11.
Caussian-distribuut observations noise (p( (p(Y.d.)
V2nga
exp
2n
with some
noise standard-deviation Ø.
Electrophysiology, however, records trains of action potentials (spikes) that are often modeled as a Poisson processes. Binned spike counts are thus typically approximated as a Poisson random variable. However, when firing rates (FRs) are high enough, the Gaussian assumption is an adequate
31

---

Published as a conference paper at ICLR 2025
approximation, and thus the Gaussian assumption remains common in data-driven neuroscience models (c.g.. Linderman et al. (2016); Aoi & Pillow (2018)). However, some species (e.g., bats (Allen et al. 2021)) exhibit extremely low FRs, which may require Poisson likelihood functions. Thus we propose an extension to CREIMBO that with a Poisson likelihood function for ensemble and traces inference, that replace the Gaussian model in low spiking-rates scenarios.
Our Poisson extension first involves transforming the spike-timing data to binned spike counts for all sessions d 1... D (2))) Next, we consider the same goal of finding the ensembles A and their traces (XT) under the assumption that the data follows Poisson distribution.
To simplify, we focus on a single session d, with = (Y (T for transpose), and note that the extension below applies all sessions. In the Poisson model, for each neuron 1... N, the likelihood of certain spike count at time given a latent rate d is:
t.n A
pl P(ye,
yt,!
where V denotes the factorial of Ve,..
Here we model the unknown rate TxN as the low-dimensional representation captured via the
ensembles A E RN and their traces क == (X4)²* E RR'x (T for transposc), i.c., À X such that the likelihood becomes
ply A)
II
al De,n)
n=1 t=1
II 1= 1=
P P(te.. m [AA" [ ((,m)
N (GAT
Ve.n (I..)
ye.n!
#=1
This can be trained by iteratively updating ã and to minimize the negative log of the above likelihood, which can be represented mathematically as:
[ [ Ie. Me ((EAA le,a)
( A)
arg min
log
IIII
Vt
arg min A
[ I. Ven log( A [t,n]]
where the logarithm is taken using the natural exponential base. Notably, in the above arg min, chose to omit the constant term that emerge from taking the logarithm (i.e., (og(v...!) as it does not affect the argument minimization.
In this Poisson case, in contrast to the Gaussian, we can no longer solve for æ using e.g., least squares LASSO; instead, will update z via Gradient Descent. The first step is compute the gradient of the cost function in equation 4 with respect to ă.
To simplify the calculation, we first notate the two components obtained in equation 4 by the auxiliary terms
gi A) = Σ A(EAT1 1
f,m
and
( x. A) = Σ
Vt.n log( A le..).
1,m
The cost function from equation thus be written in terms of these functions
32

---

Published as a conference paper at ICLR 2025
( A}
arg min (ā, A) + 92 ã. A)]- a.A
(5)
Hence, for updating ø via Gradient Descent, we need to find the gradients of gi and g2.
We first establish two notations to make subsequent steps clearer:
let (m).)) denote a matrix of ones of size m x k.
let P),. denote matrix of shape M x P whose entries are all zeros except for the entry at index (m,p), p), which is set to 1.
Calculate
0g () Ohe
To begin, we can rewrite gi A) as 91 (z, A) Σ A In ( (
and follow the identity (taken from Petersen (2008) Eq. 20):
a Mb
ab"
where a, vectors, and M matrix.
(6)
Building on this, let (1),, correspond to the af vector from equation AT (T E RRx1 correspond to the vector of equation 6 and correspond the matrix M Then, the gradient of 91 (z, A) with respect to 2 is:
dgi (ã, A)
O ((Hla,,,āA"l (4)o)) 1 A = [1(.N
Oæ
dz
Calculate
Bir
First, will rewrite the term
[TA?] fnn
EA Jen
Σ
(i 6(1,7), EA 6(N,1)
(7)
(,,T),, is vector of length T with value of index t and O elsewhere, and (N,1)).. is column vector of length with a value of index n and 0 elsewhere.
We will use again the identity from Petersen et al. (2008) presented in equation 6 and now mark 5(1.T),,, with the equation 6 vector a i is again M, and A & When applying the identity in equation 6 to equation we thus obtain
7
T A I.n
5(T,1),, 61,N), 8(T,N) A
d
Given the above, we receive
a EN
ΣT 1=
Ven log (|ãA² (1,)))
8gp Bg2(y, 3. A)
(8)
N
T
Von
6 (T.,). A
Σ
Wee 6 (T.
Σ
A
1
v A
33

---

Published as a conference paper at ICLR 2025
where Pe,n
Ms. and V RTXN is matrix whose entry index (t, n) is dhe,n-
Going back to the expression in equation the full gradient of the Poisson loss with respect to æ is:
Ølg1 (z, A) + 92 क, A) [)]e,,)) VA (HHer,N) 4)
Oa
and the projected gradient step in x at iteration m thus follows:
m+ IT - - (Iler,N) v) A)
where 1I can project the columns of ã onto some constraint set, with the identity transform used if no projection is desired, and 17" the gradient descent step size at iteration m.
For inferring A, each iteration m, can then the Spiral-Tap package Harmany et al. (2012) to infer the sparse ensembles.
The above ensemble model (including A and æ) is the only part that directly connects the Poisson- distributed spike-counts (Y"), while the latent evolution of æ via F F is independent of the Poisson assumption the spike-counts y Particularly, as the transition from Dt to 21++1 is dictated by the ensembles' dynamic interactions (F, = Σx=1 Ckr which we continue to assume follows Gaussian distribution in the latent space, the dynamic inference part of CREMMBO remains unchanged from the main text.
34

---

Published as a conference paper at ICLR 2025
CrEIMBO
f
Amygdala left
Amygdala right
Dorsal anterior cingulate cortex left
Dorsal anterior cingulate cortex right Hippocampus left
fg
fg
f10
f11
f12
Hippocampus right
supplementary
Pre supplementary motor area right
Ventral medial prefrontal cortex left
Ventral medial prefrontal cortex right
Wiirett Dynumic Price
f
fa
fs
f6
(c))
[cg)
(cg)
T8
TO
110
111
4 (c4)
4 (cs)
4 [cc)
f (c))
(cc)
Scatn:: 1SMAO2)
[Gg)
fo.
ho (Gra
4 ((1)
f (Gz)
Fo
Sowarr (OWEET)
f
f2
fs
fo
fa
f9
110
f11
f12
f13
f14
f10
Sauson SANREO
fg
F10
Regjons Stngl Session
f
fs
ho
Figure 24: Comparison of sub-cireuits identified by CREIMBO with those identified by other ap- proaches using the real-world human neural recordings. CREIMBO reveals distinct motifs capturing multi-regional flows specific areas, including both cross-regional and multi-regional interactions. In contrast, other approaches identify sub-circuits with overlapping flows that fail to capture the full spectrum multi-regional interactions.
35

---

Published as a conference paper at ICLR 2025
11
Figure 25: Identified dynamics' coefficients (c for the human data experiment under increasing noise levels demonstrate robustness, with rapid increase in internal pattern frequency at higher noise levels.
36

---

Published as a conference paper at ICLR 2025
Figure 26: Number of neurons per area across all files
37

---

Published as a conference paper at ICLR 2025
A
SUSISEa5
Figure 27: number of neurons per area motivation for area-focusing. A: Overall number of neurons per sub-area (colored by parent area) across sessions. B: Distribution of areas with sub-area name across all sessions. C: Distribution of areas with parent-area name across all sessions. D: Number neurons per parent-area (across the 10 selected areas) in the 19 used sessions.
38

---

Published as a conference paper at ICLR 2025
Thiarna
Texel
Hasocarees
Prirary
Time (seconds)
Figure 28: Example Session Chen et al. 2024a). A: Data firing rate over time separated by trials in vertical lines. the color of the vertical line captures the relevant brain area (example random session). B: Trials start and end times (in seconds) for examples session. C: Trial duration distribution for example random session.
39

---

Published as a conference paper at ICLR 2025
A
Thalamy
0
A Awswood
A
Secondary
Orbital
Offactory
Secondary mo 68
Orbital-
Basal genglia
Olfac 193
151
Brain stem-
Brain stem-
Secondary metor-
Frontal
230
250
Offactory-
Frontal
Primary
Brain stem-
388
417
B
ensembbles
Frontal-
515
Primary metor
A-car
Thalamus
Secondary
Basal ganglia
Olfactory
Brain stem
ensembles
Somatosensory
Frontal
Hippocampus
Primary motor
C
f2.
f
5
f6
f8
f9
f 10.
D
f1
f2
f3
f4
fs
0
10
0
f6
f7
f8
f9
f10
Identified components of the Mesoscale data (Chen al. 2024a) using 19 random sessions with 40-60 trials within each. A: Identified ensemble matrices for three random sessions. B: Correlations between circuits. C: Between-arree arees interactions (extracted from | via block off-diagonal mask) D: Within-areat areas interactions (extracted from {fx 1K via block diagonal mask)
Figure 29:
40

---

Published as a conference paper at ICLR 2025
SLDS s all together
66est
6dcal
start
trial
early
miss left right
6teal
end
trial
Ignore
early
B
Example Example
SLDS single session
110
Figure 30: The components identified by the SLDS (Fig. 30 and rSDLS (Fig. 31 baselines Lin- derman et al. 2016 2020), when we applied them to the mouse mesoscale data from all sessions together (A) individual sessions (B). For each model, we showed both the sub-cireuits and their active coefficients for the first few trials. Compared to CREIMBO, the identified sub-circuits are more distributed and dense, and the dynamic coefficients capture only sharp transitions, i.e., unable to represent multiple co-occurring processes, in contrast to CREIMBO. Checking individual sessions with SLDS (B) rSLDS (C), resulted in per-session set of sub-cireuits (Fs) that cannot enable cross- session sub-circuits comparison. All models ran with K = 10 sub-circuits and pj 3 ensembles per-region, keep maximum consistency with our CREIMBO results from Figure 29
A
rSLDS tolgethrr
Bica)
Mel
6lcal
no early
miss left right
G(ca)
Gical
end
trial
Ignore
*
Mccl
early
B
Example
single session
rSLDS single
TRNGALLS CRAM
Example 2
Figure 31: Same as Figure 30 but for the rSLDS baseline Linderman et al. 2016
41

---

Published as a conference paper at ICLR 2025
A session
session 2
session 3
session 4 session session 6 session 7
1
11
1-
X1
1
11
1
X1
session session 9 session 10 session 11 session 12 session sessino 14 session 15
1-1
11
1-1
X1
1-1
42
1-1
11
1-1
1-1
X1
1-1
12
1-1
B
A1
A2
A3
A4
A5
A6
A7
A
8
104
ALE
S1B
06N
6
Area
Area
Area
Area
Area
Area
Area
Area
Ag
A
A
A
A
A 14
415
Area
Area
Area
Area
Area
Area
Area
Figure 32: Cont. of Figure 11
A: Ground truth latent dynamics ((X Já=l
B: Ground truth
ensemble compositions (fA 1215.
42

---

Published as a conference paper at ICLR 2025
A
session
session 2
session 3
session 4
0.5
0.0
session 5
session 6
session 7
session 8
1.0
0.5
session 9
session 10
session 11
session 12
1.0
0.5
session 13
session 14
session 15
1.0
0.5-
C2
C3
100
200
100
400
500
100
200
300
400
500
160
200
300
400
500
Time
Time
Time
p(x.3)
P(A.A)
PtP. FF)
mcan median
114 All
111
IF F11
Mly S11
Figure 33: Cont. of Figure A: Ground truth sub-cireuits coefficients (c) vary between sessions. B: Correlation between ground truth and identified components for A (left), x (middle left). (middle right), Y (right). C: l distance between the components learned by CREIMBO and the ground-truth components for A (left), (middle left), F (middle right), Y (right).
43


Reviews:
[{"review_id": "6F0lzFC8OJ", "reviewer": "Reviewer_8qa1", "summary": "Through the study, the authors have proposed a novel analytical approach (named as “CREIMBO”) for learning dynamics of latent representations from high-dimensional electrophysiology. The major advances of CREIMBO comes from that compression of high-dimension and extraction of dynamics were conducted simultaneously in CREIMBO, while keeping the interpretability. Through experiments with synthesized- and real data, the authors have properly demonstrated the validity of CREIMBO in the study. As CREIMBO contains conceptual novelty and its validated effectiveness, I believe this model can be one of useful candidate models for studying high-dimensional, partially overlapped, data, such as intracranial EEG or multi-array spike data. Thus, CREIMBO will be useful to neuroscientists.", "questions": "1. Please provide details of how synthetic data was generated. \n\n2. In Fig. 2M and O, I observed the high inconsistency in reconstruction accuracy in other models, e.g.,  mp_rSLDS_Gauss (per tial). This is synthetic data and, I assume, each session contains similar level of noise. It is very unclear how other baseline models nearly completely fails but performed near perfectly in some sessions. My question is, does this inconsistency supposed to support the robustness of CRIEMBO?  \n3. With real data, the authors demonstrated the validity of CRIEMBO. While I agree with the author’s claim, it does not necessarily lead to the strength of CRIEMBO. Are there any interesting differences or unique patterns that can be obtained from CRIEMBO vs. other baseline models? \n4. The authors checked the robustness of CRIEMBO over different noise levels, with real data. More common practice is evaluating the robustness of the models with synthetic data, one with grountruth. If there were any specific reasons why it was done with real data, please specify it.", "limitations": "", "rating": 8, "confidence": 4, "soundness": 4, "presentation": 3, "contribution": 3, "strengths": "The major strength of this study comes from its conceptual advances embedded in the proposed CRIEMBO. While there have been efforts to model the high-dimensional brain dynamics via deep learning, many of models failed to yield interpretable features, which is the most important aspect in the neuroscience field. Thus, the neuroscience field still tends to rely on relatively simpler yet interpretable models. In this regard, I believe CRIEMBO proposed in this study can be a good solution for this gap, clearly upholding the originality of this study. The authors thoroughly examined the validity of the proposed model using simulated- and experimental data, leading to the high quality of this work. Plus, the clarity of this paper is relatively high as the model was well described in the manuscript, whereas the reviewer believes there are some rooms requiring the attention of the authors. Altogether, the scientific significance of this paper is clear and can be interesting to the electrophysiology field.", "weaknesses": "While the overall strength of this study is obvious, there is a major weakness: Insufficient details in simulation study with synthetic data. It is unclear how the synthetic data was generated, e.g., what was the noise level, what kinds of parameters used for. Due to this uncertainty, it is nearly impossible to understand some of the intriguing findings in this study, especially with synthetic data. For example, there is very high inconsistency in performance across sessions in other methods, but not for CRIEMBO (Fig. 2M). While it can be interpreted as robustness of CRIEMBO, it is also possible the choice of benchmark models was not optimal for this type of synthetic data. Related to it, there was no comparison work with real data. Thus, the superiority of CRIEMBO needs further validation.", "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": "", "extra_scores": {}}, {"review_id": "QvG0EBoITH", "reviewer": "Reviewer_jy8o", "summary": "The recent developments in neural recording technologies allow recording from large populations of neurons from multiple brain regions. Latent space models are often used to analyze these datasets, but they are generally limited to the study of single or few populations of neurons recorded simultaneously. To overcome these limitations, the presented work introduces a new algorithm that can capture variability across recording sessions and across and within brain regions. The method assumes a shared latent representation across areas and structured priors given each session. The authors validated the method on simulated data and neural data, showing that the model can capture variability across and within brain regions.", "questions": "Can this model be applied to other organisms, like mice or rats? Extending the application to neuroscience animal models would greatly increase the impact of the presented work.\n \nHow difficult would it be to extend this model to use a Poisson observation model to better capture neural spiking activity?", "limitations": "", "rating": 6, "confidence": 4, "soundness": 3, "presentation": 4, "contribution": 2, "strengths": "The manuscript is clearly presented, referenced in the context of the relevant work and technically sound. The method was tested in simulations coming from the generative model, where the model was able to recover the parameter settings, and was able to better capture the variability compared to existing models. They also tested the model in neural recordings identifying functional connectivity between and within brain regions. Moreover, they tested the robustness to increasing levels of noise in the data.", "weaknesses": "Further work in needed to fully understand the impact of this work. The authors motivate their work from their ability to better understand behavioral tasks from asynchronous recordings of spiking neural activity. However, the authors only tested the model in a single dataset where they limit the analysis to functional connectivity. It would be relevant to assess behavioral variables, such as ability to decode the present image stimuli from the learnt representations. Since the emphasis of the model is brain functional connectivity, the method should be compared to alternative recording methods such as fMRI or the LFP present in the dataset. Moreover, the comparison to alternative models in the neural data is limited to one qualitative test. A comprehensive quantitative comparison, such as decoding or reconstruction performance, is needed to understand the capabilities of the proposed method. Along the same lines, it is not surprising that the model outperforms alternative methods when the simulated data is tailored to the given model. A more relevant comparison would be simulating neural data from a neural process with temporal, task and/or behavioral variability and fit the different models there. It would also be relevant to highlight the model strengths and weaknesses in simulated data. One of the motivations for this work is its application to spiking data, but the Gaussian assumption limits its application to this kind of observation model and must be handled in preprocessing. While the authors verbally list this and other limitations in the discussion, it would be informative to show the limitations and, more importantly, capabilities of the model as results.", "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": "", "extra_scores": {}}, {"review_id": "J7xOgMSuuL", "reviewer": "Reviewer_eTfH", "summary": "The authors propose CREIMBO that can extract interpretable neural subcircuits underlying multiregional neural signals by utilizing multi-session recordings that can have a variable number of recording units, trials, and durations. Through simulations, they show that their model successfully uncovers the ground truth dynamics when multi-session recordings are modeled together. In real data analysis, authors show that identified subcircuits can be sparse indicating specialized functionality, and reveal across-region interactions.", "questions": "- In line with my previous comment in weaknesses item 2, the authors show subject 10 in Fig. 3 which has a small number of available neurons compared to some other subjects (for Screening task). Are the example subcircuits and latent states extracted for the Screening task? Do the identified subcircuits and A matrices look denser for a higher number of neurons?\n- In the ablation study in Fig. 2 ‘All regions (sparse)’ A matrix, authors show that the CREIMBO cannot infer the true underlying connectivity matrix, which makes sense since the inductive bias on block-diagonal A matrix is removed. Do authors use the same A matrix initializations in this ablation study? If not, how would the results change if the same block-diagonal initialization is applied for ‘All regions (sparse)’ case? Can authors provide an intuition on why inferred latent factors deviate significantly from the true latents? Would the same hold in the K=1 case, in which, a similarity transform would exist between not block diagonal sparse, not block diagonal not sparse, and block-diagonal sparse A matrices? Also, are the authors showing trial and session averaged latent states in these figures? If so, how does single-session latents look like for CREIMBO?\n- The block structure imposed on A matrix seems like no regions have shared latent states and interregional interactions are captured by temporal dynamics. Did the authors try having some latent factors shared across all regions? How would it change the performance?", "limitations": "", "rating": 8, "confidence": 4, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "- Authors provide a new perspective over existing models of multi-regional neural models by allowing their model to leverage multi-session recordings that can help extract global and robust interregional interactions. While doing that, they keep the interpretability in-tact unlike deep-learning approaches. \n- Authors performed exhaustive simulation experiments to show the importance of their model formulation. \n- The paper is well-written and easy to follow. The proposed model architecture and training framework are intuitive and effective.", "weaknesses": "- Even though the paper is well written, the subfigures are very small and too crowded, which makes it hard to understand. Also, I think captions/labels for some appendix figures such as Fig. 11 should be improved. \n- It seems like the number of ensembles per region ($p_j$) is an important hyperparameter for CREIMBO. I wonder if the model would reveal some consistent subcircuits with small $p_j$, such that these subcircuits explain most of the variance and are consistent across sessions and subjects. Also, how would training times vary with large $p_j$ values such that $p \\approx N$? Overall, I think the scalability of such a model is an important aspect since modern neural recordings can include hundreds of neurons from one region, in such case, max($p_j$) = 7, can limit the interpretability of the identified subcircuits. \n- In simulations (Fig. 2F), the authors show that the single-session model underperforms CREIMBO by a large performance gap. This can be caused by a small number of trials in each simulated session and short trials, but I could not find these details in the text (if it is in Fig. 11, I think it requires explanations in the caption). If this is indeed the case, I wonder how single-session models' performance would increase with longer sessions.  \n- I think the biggest contribution of CREIMBO is its multisession modeling over other approaches like mDLAG. However, the modern recording session can have hundreds of trials of data that can be sufficient to train models to understand multiregional dynamics. Therefore, I think it would be nice to see how their model compares to mDLAG even if it operates on a single-session basis. Based on Fig. 22, for the real dataset considered in this study, using multiple sessions in modeling seems important, and a comparison to mDLAG would highlight the importance of multisession modeling. Even though mDLAG does not learn dynamic matrices for temporal evolution, its learned readout matrices and lag parameters would still indicate interregional interactions, and I wonder if interregional interactions learned by mDLAG would be as poor as 'Session # (SPARSE)' in Fig. 22. Also, did authors compare their model to SLDS variants in real data as done for simulations? Overall, I think this work would benefit from more baseline comparisons to existing approaches in real data.", "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": "", "extra_scores": {}}]

Author rebuttal:

Decision: Accept (Spotlight)