CREIMBO: Cross-Regional Ensemble Interactions in Multi-view Brain Observations

## Abstract

Abstract Modern recordings of neural activity provide diverse observations of neurons across brain areas, behavioral conditions, and subjects—thus presenting an exciting opportunity to reveal the fundamentals of brain-wide dynamics underlying cognitive function.
Current analysis methods, however, often fail to fully harness the richness of such data as they either provide an uninterpretable representation (e.g., via “black box” deep networks) or over-simplify the model (e.g., assume stationary dynamics or analyze each session independently).
Here, instead of regarding asynchronous recordings that lack alignment in neural identity or brain areas as a limitation, we exploit these diverse views of the same brain system to learn a unified model of brain dynamics. Specifically, we assume that brain observations stem from the joint activity of a set of functional neural ensembles (groups of co-active neurons) that are similar in functionality across recordings, and propose to discover the ensemble structures and their non-stationary dynamical interactions in a new model we term CrEIMBO (Cross-Ensemble Interactions in Multi-view Brain Observations).
CrEIMBO identifies the composition of the per-session neural ensembles through graph-driven dictionary learning and models the ensemble dynamics as a latent sparse time-varying decomposition of global sub-circuits—thereby capturing non-stationary dynamics.
CrEIMBO identifies multiple co-active sub-circuits while maintaining representation interpretability
due to sharing sub-circuits across sessions. Furthermore, CrEIMBO distinguishes session-specific from global (session-invariant) computations by exploring when distinct sub-circuits are active.
We demonstrate CrEIMBO’s ability to recover ground truth components in synthetic data and uncover meaningful brain dynamics, capturing cross-subject and inter- and intra-area variability, in high-density electrode recordings of humans performing a memory task.

## 1 Introduction

Identifying the interactions between and within brain areas is fundamental to advancing our understanding of how the brain gives rise to behavior.
Recent advances in neural recording present an exciting opportunity to study brain-wide interactions by enabling simultaneous recording of neural activity across many brain areas through multiple high-density electrodes.
Such experiments, repeated over many sessions with different implantation patterns, offer multiple asynchronous recordings of the brain system, with each session encompassing hundreds of distinct neurons across regions.

Individual recording sessions provide a 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 to 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 session-by-session, may hinder our ability to distinguish between computations that are session or subject specific from task related or global session invariant neural processes.

Compounding the difficulty of merging different sessions’ data are the non-linear and non-stationary nature of neural activity, as well as the presence of noise and variability across trials—both of which require additional model complexity to account for.
Moreover to extract scientifically meaningful insight, the models must be interpretable in that model parameters must be readily relatable to task variables or connections between recorded units.
This gap between the opportunities offered by modern neural data and the limited capabilities of current methods necessitates new approaches to leverage the richness in modern brain data and discover the fundamental neural sub-circuits governing brain activity.

One current approach used to combine multi-session data into a single model is the use of deep learning models. By training to predict all sessions’ data, deep learning models implicitly merge the datasets into the fit network weights.
While these powerful models can combine information across sessions, their uninterpretable ‘black box’ 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., [, ]) may be limited in their ability to merge sessions of varying durations or identify the interactions between the identified neural components.

Here, we re-frame cross-session variability in neurons recorded and experimental settings not as a drawback, but as a valuable advantage. Rather than viewing the different sessions as capturing distinct neural populations, we perceive such data as providing multiple, complementary perspectives of one brain system. This approach enables a more thorough and holistic discovery of the underlying system designed to extract joint information from the entire dataset.

In particular, we hypothesize that
the population dynamics at each session is driven by a common set of global sub-circuit interactions shared across sessions. To capture the non-linear and non-stationary aspects of the dynamics, we leverage recent advances in decomposed linear dynamical systems (dLDS), which enables us to model how the neural trajectories transverse a low-dimensional manifold via a time-varying combination of the sub-circuits.
While we assume that the representation of the latent dynamics is shared across sessions, the measured neural ensembles are not. We thus model each session’s ensembles as independent neural representations, i.e., the membership of neurons to ensembles, that captures how the latent trajectories project into each sessions’ view of the brain.

In this work we lay out this novel model, which we term Cross Ensemble Interactions in Multi-view Brain Observations (CrEIMBO), 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, same-area neural ensembles, deploying a specialized structured prior over per-session projections from circuit space to neural space.
- •
We accurately recover ground truth components in synthetic data and discover meaningful cross-regional brain interactions underlying human brain data utilizing a shared low-dimensional circuit space across recording sessions.

## 2 Related Work

Over the past decade, remarkable advances in neural data acquisition technologies have enabled the measurement of numerous neurons across brain regions, subjects, and experimental settings. This has spurred follow-up work to leverage this extensive data for novel neural discoveries, including merging non-simultaneous sessions and examining cross-regional brain interactions, all while accounting for the complex, non-linear, and non-stationary nature of neural activity.

Merging non-simultaneous neural recording sessions.
Existing work for combining non-simultaneous neural population recordings, here termed “sessions”, include
[] who introduced the term “stitching” for this problem and modeled neural activity with a shared latent dynamical systems model and [] who proposed to
stitch asynchronous recordings for inferring functional connectivity by assuming shared cross-session latent input and demonstrated the advantage in observing more neurons for inferring connectivity.
However, none of these methods capture multiple ($>2$) interacting neural sub-circuits active simultaneously spanning populations across multiple areas, nor do they distinguish between within- and between-area interactions, nor do they address temporal non-stationarity.
In [, ], Bishop et al. proposed integrating the neural activity structure into inference and expanding the proposed applications of stitching to study of communication between two neural populations of brain areas and proposed a theoretical framework of when stitching is possible.
In [], the authors proposed
extracting low-dimensional dynamics from
multiple sessions by learning temporal co-variations across neurons and time.
Both of these methods, through, assume at least a small overlap in the identity of cross-session observed neurons.
LFADS [, ] enables stitching
while addressing the need for neural overlap
by finding a shared non-linear dynamical model.
However, the latent dynamics may be uninterpretable in terms of interactions among neural populations, and due to its foundation on a variational auto-encoder, it may be challenging to tune.

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 capture non-stationarity through various methods such as ’black box’ deep learning models (e.g. [, , ], which are powerful but are often not interpretable with respect to neural interactions; linear models, which are interpretable but not expressive enough to capture the complex nature of neural data; or switching piece-wise linear models  [, , ], which fail to enable the identification of co-active dynamics.
dLDS (decomposed Linear Dynamical Systems) [] (with its improved learning procedure in []), addresses this gap by generalizing the switching assumption to sparse time-changing decomposition of linear dynamical systems, thus enabling
capturing co-active sub-circuits.
However, while the current dLDS framework can be modified for co-processing multi-subject recordings that present variability in neuron identity and distribution, this modification lacks the structured regional prior over the loading matrix that is necessary to partition multi-area recordings such that between and within area interactions can be isolated and studied.

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 refereed to as groups of neurons with similar functionalities—that determine the axes in the latent space. For instance, the mDLAG model used Automatic Relevance Determination (ARD) to promote population-wise
sparsity patterns of each latent variable.
In [], the authors presented an approach to address sparsity in loading matrices while accounting for multi-trial and task-conditions variability. They proposed a graph-driven regularization technique that facilitates the grouping of co-active neurons into the same ensembles and pushing apart neurons that present different functionality, while allowing for structural adjustments in ensemble configurations across different task conditions. They further enable per-trial variability in temporal activity of the ensembles. However, their model does not model the dynamic interactions between ensembles or give a closed form prior for the evolution of the ensembles’ temporal traces.
While other methods, e.g., clustering approaches like  [] or Sparse Principal Component Analysis (SPCA) [], 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 support controlled adaptations of structural membership across sessions.

Multi-regional brain interactions.
With the ability to record from multiple brain areas and the recognition of the importance of considering multi-area communication [], recent works have delved into understanding the communication between or within brain areas.
In [, ], the authors proposed models to identify time delays between cross-regional brain interactions, overcoming the challenge of distinguishing indirect, concurrent, and bi-directional 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 [] using dimensionality reduction, or by Generalized Linear Models (GLMs), with either Poisson [] or Gaussian [] 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, including [, , ] discussed the potential usage of their model to multi-regional interactions, have proposed studying multi-regional interactions using switching dynamical systems approaches, where the model switches between different states over time to capture distinct periods. However, this approach cannot disentangle interactions among multiple sub-circuits active simultaneously, nor can it identify interactions at multiple temporal resolutions that occur concurrently and distinguish between these or multiple co-processed processes whose encoding occurs gradually rather than through sharp transitions.

## 3 Problem Formulation and Approach

Problem Formulation

Let $\{\bm{Y}^{d}\}_{d=1}^{D}$ be a set of estimated firing rates over $D$ asynchronous sessions, indexed by ${d=1\dots D}$. The neural recordings of each $d$ session, $\bm{Y}^{d}\in\mathbb{R}^{N^{d}\times T^{d}}$ captures $T^{d}$ time observations of the activity of $N^{d}$ neurons from up to $J$ distinct brain areas, thus offering a partial view of the brain system. Note that sessions may refer to the same or different subjects, can vary in duration or number of neurons, i.e. $N^{d}$ and $T^{d}$ need not be equal to $N^{d^{\prime}}$ and $T^{d^{\prime}}$ for different sessions $d,d^{\prime}$, and may differ in the sub-set 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—hindering the direct application of existing analyses as tensor factorization.

We assume that these brain observations reflect the hidden activity $\bm{X}^{d}\in\mathbb{R}^{p\times T^{d}}$ of $p<<N^{d}$ functional neural groups that evolve and interact over time. Each of these groups encompasses a sparse set of same-area neurons, while each neuron can belong to more than one group with varying degrees of cross-group membership, encoded in the values of $\bm{X}^{d}$.
The low-dimensional latent space in which the groups interact is shared across sessions and projected to each session’s observed neurons via an unknown per-session projection $g^{d}:\mathbb{R}^{p}\rightarrow\mathbb{R}^{N^{d}}$
that captures the groups’ compositions in the visible neurons, such that $\bm{Y}^{d}=g^{d}(\bm{X}^{d})$.

The interactions between these functional groups follow a per-session non-linear and non-statioary dynamics $\bm{F}_{t}^{d}:\mathbb{R}^{p}\rightarrow\mathbb{R}^{p}$, where $\bm{X}^{d}_{t}=\bm{F}^{d}_{t}(\bm{x}^{d}_{t-1})$. As in general non-stationary systems cannot be fit without constraints, we follow the dLDS model and assume that these dynamics can be described by a set $\{\bm{f}_{k}\}_{k=1}^{K}$ of $K$ global time- and session-invariant “sub-circuits”. Here each sub-circuit ${\bm{f}_{k}\in\mathbb{R}^{p\times p}}$ represents a basis interaction between these groups that 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 with a single underlying mechanism.

The problem we address in this model is identifying the set of unknown multi-regional sub-circuits, including their neuronal composition and non-stationary activation levels, from asynchronously collected observations.

CrEIMBO (Cross Ensemble Interactions in Multi-view
Brain Observations)

A naive approach to identifying the unknown sub-circuits involves first projecting each session’s data into a low-dimensional space and subsequently fitting the dynamics model in the low-dimensional space. However, this two-step approach can be sub-optimal as the dimensionality reduction projections may not prioritize the same dimensions needed to preserve dynamic fidelity.
An additional advantage of simultaneously identifying the ensembles and dynamics is that in the cross-session model it automatically aligns the cross-session components based on the shared dynamics prior. Such alignment is unlikely given the sensitivity of the subspace ordering in dimensionality reduction methods to each session’s distinct statistics.

We therefore present CrEIMBO, an algorithm that simultaneously fits the global sub-circuits that underlie the inter-ensemble interactions across multiple sessions.
Our model is predicated on three key assumptions that
1) ensembles of neurons, rather than individual neurons, form the basic units that interact in brain dynamics,
2) it is feasible to identify functionally analogous neuronal ensembles across sessions, and 3) the full repertoire of interactions between these ensembles are governed by a global set of dynamical primitives in the form of a basis of linear dynamical systems.

Let $\bm{A}^{d}$ be a sparse matrix where each column, $\bm{A}^{d}[:,j]$, encodes the neuronal composition of the $j$-th ensemble in the $d$-th session, such that non-zero values within these columns represent the membership of neurons in the ensemble and their magnitudes reflecting the degree of membership.
Since we aim to capture multi-regional interactions, we design the ensemble matrix $\bm{a}:=\bm{A}^{d}$ as a block diagonal matrix where the $j^{th}$ block, $\bm{a}^{j}\in\mathbb{R}^{n_{j}\times p_{j}}$, contains the neural ensembles for the $j$-th area. We denote $p_{j}$ as the number of ensembles within that area (note that $p=\sum_{j=1}^{J}p_{j}$) and $n_{j}:=N_{j}^{d}$ as the number of neurons in that area (i.e., the total number of neurons observed in session $d$ is $N^{d}=\sum_{j=1}^{J}n_{j}$).

Furthermore, let $\tau:=\mathrm{T}_{m}^{d}$ indicate the time duration (i.e., the number of time-points) of the $m^{th}$ trial, for $m=1\dots M_{d}$, in session $d$. In each trial, the temporal activity for the $p$ ensembles are denoted as $\bm{x}:=\bm{X}_{m}^{d}\in\mathbb{R}^{p\times\tau}$, leading to a model where the observations $\bm{y}:=\bm{Y}_{m}^{d}$ in trial $m$ are assumed to arise from the joint activity of all these ensembles, up do a normally-distributed error $\epsilon\sim\mathcal{N}(0,\sigma)$, i.e., $\bm{y}=\bm{a}\bm{x}+\epsilon.$
The evolution of the ensembles’ activities ($\bm{x}$) reflect the latent interactions between ensembles, which we model with a non-stationaly linear dynamical system $\bm{x}_{t}=\bm{F}_{t}\bm{x}_{t-1}$ for $t=1\ldots\tau$. The time-changing transition matrix $\bm{F}_{t}\in\mathbb{R}^{p\times p}$ captures both non-linear and non-stationary brain activity.
A key desired property of our model is to capture the overlapping activity of multiple sub-circuits of interacting ensembles. We thus follow the dLDS concept and model the interactions ($\bm{F}_{t}$) through a time-varying sparse decomposition ${\bm{F}_{t}=\sum_{k=1}^{K}\bm{c}_{kt}\bm{f}_{k}}$, where $\{\bm{f}_{k}\}_{k=1}^{K}$ represents the $K$ global ensemble interactions, and $\bm{c}_{kt}$ are the time-varying coefficients that represent the modulation of each interaction pattern over time.

Fitting this model requires identifying the global operators $\{\bm{f}_{k}\}_{k=1}^{K}$ and the functional ensembles $\{\bm{A}^{d}\}_{d=1}^{D}$.
For this, we employ an alternating approach, where we iterate between updating the global operators, $\{\bm{f}_{k}\}_{k=1}^{K}$, and infer,
for each session $d$, 1) the ensemble compositions $\bm{a}:=\bm{A}^{d}$, 2) their dynamics $\bm{x}$, and 3) the sub-circuits’ temporal coefficients $\bm{c}_{t}$. These steps are iterated until the parameters converge.

Ensemble update: We update $\bm{a}$ per-row (neuron) $n$ ($\bm{a}_{n}$). Following the work of [], we infer the sparse ensemble structures with graph-driven re-weighted $\ell_{1}$ regularization that groups together neurons with shared activity patterns, while pushing apart neurons who do not behave similarly. The graph for each area $j$, $\bm{h}^{j}\in\mathbb{R}^{n_{j}\times n_{j}}$, is calculated based on a data-driven Gaussian kernel that measures the temporal similarity between neurons, such that $\bm{h}^{j}_{n1,n2}=exp(\frac{\|\bm{y}^{j}_{n1}-\bm{y}^{j}_{n2}\|_{2}^{2}}{
\sigma_{h}})$ where $\sigma_{h}$ is a hyperparameter.
Specifically, for each row $\bm{a}_{n,:}$ and considering only the ensembles in region $j_{n}$ that neuron $n$ is part of,

$$ $\displaystyle\widehat{\bm{a}}_{n,:}=\arg\min_{\bm{a}_{n:}}\|\bm{y}_{n:}-\bm{a} _{n:}\bm{x}_{j_{n}:}\|_{2}^{2}+\sum_{j=1}^{p}\bm{\lambda}^{d}_{n,j}|\bm{a}_{n, j}|\quad\textrm{ where }\quad\bm{\lambda}_{n,j}=\frac{\beta_{1}}{\beta_{2}+| \widehat{\bm{a}}_{n,j}|+\beta_{3}|\bm{h}_{n:}^{j}\widehat{\bm{a}}_{:j}|}.$ $$

Here, $\bm{x}_{j_{n},:}$ is the activity of the ensembles of the $j_{n}^{th}$ region and $\beta_{1}$, $\beta_{2}$, and $\beta_{3}$ are scalars that control the effect of the graph on the regularization. The graph adjacency $\bm{h}_{n:}\widehat{\bm{a}}_{:j}$ weights the sparse membership in the $j$-th ensemble based on temporal similarity. Specifically, if $\bm{h}_{n:}\widehat{\bm{a}}_{:j}$ is high (i.e., the neighbors neuron $n$ on the graph are in the current estimate of the $j$-th ensemble),
$\bm{\lambda}_{n,j}$ is small, resulting in lower sparsity regularization on neuron $n$ for that ensemble, promoting it to be included. Alternatively, if this correlation is small, the regularization weight $\bm{\lambda}_{n,j}$ is large, and neuron $n$ will be more likely not to be included in ensemble $j$.

Latent state and dynamics coefficient update: The ensemble activities $\bm{x}$ and coefficients $\bm{c}$ are updated iterative for each time-point $t$ via the LASSO optimization

$$ $\displaystyle{\widehat{\bm{x}}_{t},\widehat{\bm{c}}_{t}=\arg\min_{\bm{x}_{t}, \bm{c}_{t}}\|\bm{y}_{t}-\bm{a}\bm{x}_{t}\|_{2}^{2}+\left\|\bm{x}_{t}-\left( \sum_{k=1}^{K}\bm{f}_{k}\bm{c}_{k,(t-1)}\right)\bm{x}_{t-1}\right\|_{2}^{2}+ \lambda_{x}\|\bm{x}_{t}\|_{1}+\lambda_{c}\|\bm{c}_{t}\|_{1}}$ (1) $$

where $\lambda_{c}\in\mathbb{R}$ and $\lambda_{x}\in\mathbb{R}$ are regularization weights.

Dynamical system update: The global ensembles interactions, $\{\bm{f}_{k}\}_{k=1}^{K}$ are assumed to be sparse (i.e. ensembles interactions are not all to all) and are identified directly by:

$$ $\displaystyle\bm{F}^{\textrm{all}}_{n}=\arg\min_{{\bm{F}^{\textrm{all}}}_{j}} \|\bm{y}_{n}-\bm{F}_{n}^{\textrm{all}}\bm{cx}\|+\lambda_{F}\|\bm{F}^{\textrm{ all}}_{j}\|_{1}+\sum_{{k_{1}},{k_{2}},({k_{1}}\neq{k_{2}})}\rho({\bm{f}}_{k_{1 }},{\bm{f}}_{k_{2}})$ (2) $$

where $\bm{F}^{\textrm{all}}$ is an horizontal concatenation of all $\bm{f}$s, and
$\sum_{{k_{1}},{k_{2}}}\rho(\bm{f}_{k_{1}},\bm{f}_{k_{2}})$ is a de-correlation term used to ensure that distinct $\bm{f}$s are not too similar, Above, $\bm{cx}_{t}:=(\bm{c}^{batch}_{t}\otimes[1]_{1\times p})\circ([1]_{1\times K}
\otimes\bm{x}^{batch}_{t})$, and $\lambda_{F}$ is the weight on the sparsity-regularization on $\bm{F}$ and “batch” means that this operation is in practice taken on a random subset of sessions for computational complexity considerations.

Figure: Figure 2: A-E: Testing CrEIMBO over 312 iterations of random initializations and parameters (see Sec. D) reveals high correlations between ground truth and fit components for all unknowns, $\bm{Y}$, $\bm{A}$, $\bm{x}$, $\bm{F}$, $\bm{c}$ (subplots A-E respectively). F-G: A comparison between the ensemble trajectory ($\bm{X}$) and the ensemble compositions ($\bm{A}$) as identified by CrEIMBO and other approaches compared to the ground-truth. Details about the other approaches are found in Section F. Additional sessions in Figure 5. H-J: $-\log(MSE)$ between the ground truth latent dynamics, ensemble composition, and sub-circuit, to the corresponding components as identified by CrEIMBO vs. other approaches. The sub-circuits were re-ordered for each method to maximize correlation with the ground-truth circuits using Scipy’s linear_sum_assignment [].
Refer to caption: extracted/5623997/FIGURES/synth_simple/may_22_figure_synthetic_data_and_baselines6.png

Figure: Figure 3: CrEIMBO identifies cross-regional neural sub-circuits underlying multi-session human brain recordings. A: Two exemplary observations ($\bm{Y}$, right) compared to their reconstruction by crEIMBO ($\widehat{\bm{Y}}$, left). B: Reconstruction performance in terms of $R^{2}$ 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: Examplary identified sparse ensemble matrices ($\{\bm{A}^{d}\}_{d=1}^{D}$) 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 ($\bm{X}^{d}$) for the two exemplary sessions. H: The identified sub-circuits. Width indicates effect magnitude. Each circle corresponds to an ensemble with its color indicating the ensemble’s area. Black arrows near $\bm{f}_{8}$ and $\bm{f}_{11}$ highlight within-region interactions. I: Pairwise correlations between within-session sub-circuit coefficients (left, middle) and average within-session correlations across all sessions with std marked by the size of a “+” marker. J: Average and std of connection values in $f$s. K: Number of non-zero connections per sub-circuit $\bm{f}_{k}$.
Refer to caption: extracted/5623997/FIGURES/supp_human/figure_human_main8_with_A.png

## 4 Experiments

### 4.1 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$ sub-circuits represented by a set of 3 rotational matrices (Fig. 4B). These sub-circuits capture interactions among $p=3$ ensembles from $J=3$ regions ($p_{j}=1$ for all $j=1,\dots,J$). Each region consists of a random number of neurons drawn from a uniform distribution between $4$ to $9$. We finally generated the data by simulating 5 sessions with different neuron identities and varying numbers of neurons per region(Fig. 4A 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 random parameters from a range described in 1. Accurate recover included all parameters: the reconstructed observations $\bm{y}$ (Fig. 2A), ensemble composition $\bm{A}$ (Fig. 2B) and activity $\bm{x}$ (Fig. 2C),
sub-circuits $\bm{F}$ (Fig. 2D),
and coefficients $\bm{c}$ (Fig. 2E).
Comparing CrEIMBO’s results to components identified by other methods (baselines detailed in F), we found that approaches that either separate the ensembles and dynamics identifications steps, analyze each session separately, or do not consider the localized structure of the ensembles, fail to recover the ground-truth components, i.e., the resulting fit model parameters are less correlated with the ground-truth observations (Fig. 2F-J, 5).

Next, we tested CrEIMBO on a richer synthetic example with $D=15$ sessions. Each session had a varying number of neurons (a random integer between $14$ and $19$ per region from a uniform distribution) and distinct sub-circuit coefficients (see Fig. 8 for the generated ground truth components and 6 for regional distributions and masking). The sub-circuits were set as rotational matrices, similar to the previous experiment. Over $204$ random initializations with different parameters (see Sec. D) we again found that CrEIMBO robustly 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. 8E) with low relative error (Fig. 8F).

We further tested CrEIMBO on more advanced data that included $D=40$ sessions with a maximum of $4$ regions and $p_{j}=3$ ensembles per region (Fig. 9). We found that CrEIMBO manages to recover the components with high correlation with the ground-truth (Fig 10).

### 4.2 CrEIMBO Discovers Multi-Regional Dynamics in Human Electrophysiology Data

We then tested CrEIMBO on human neural recordings from a high-density electrode array provided by []. The data consists of neural activity from $p=10$ brain areas in the human medial temporal and medial frontal lobes, recorded while subjects performed a screening task (detailed in []). The data encompasses $21$ subjects across $D=41$ non-simultaneous sessions. For each subject (except Subject 19), the data offers two recording sessions: 1) a “Screening” session to identify natural images that elicit strong neural responses for that subject and 2) a “Sternberg test” session to assess working memory using the images identified in the screening session. For Subject 19, the data includes only the Sternberg session.
The data encompasses multiple regions with limited cross-session regional overlap, a few missing samples per session, and variability in the distribution and count of neurons over areas, with some areas having no neurons in some trials or subjects (Fig. 11, 12).

We first converted the spike-sorted units to estimated firing rate matrices by convolving the spike trains with a 30 ms Gaussian kernel (Fig. 7). Applying CrEIMBO to the firing rate data with a maximum of $K=12$ sub-circuits and $p_{j}=7$ ensembles (Fig. 13) per region, we found multi-regional motifs with variability across sessions.

We then tested CrEIMBO on the first $18$ seconds from whole sessions together with a maximum of $K=12$ sparse sub-circuits
and as most $p_{j}=7$ ensembles per region (full parameter list available in Tab. 2).
First we note that CrEIMBO can accurately reconstruct the data with high accuracy across all sessions (Fig. 3A, B). Interestingly, the ensemble compositions identified by CrEIMBO (Fig. 3E) presented sparse patterns with one dense ensemble for most regions that we hypothesize captures the ‘general’ mean field activity of that region (the right-most column of each diagonal block in Fig. 3E). The sparse ensembles we hypothesize capture the more specialized, nuanced functionality.
The ensemble dynamics ($\bm{X}^{d}$, Fig. 3G) reveal that some ensembles are constantly active, implying that they might play an important role in the observed neural activity.

The identified set of global sub-circuits ($\{\bm{f}_{k}\}_{k=1}^{K}$) exhibits distinct localized motifs (Fig. 3H), i.e., most sub-circuits present clear trends of either sourcing or targeting the same area. This is in contrast to other approaches (Sec. F), that yield overlapping sub-circuits (Fig. 19) that emphasize only limited circuitry (Fig. 20).
Interestingly, sub-circuits 10 and 12 ended up almost empty (mainly zero-ish values, Fig. 3H, 16). This highlights CrEIMBO’s ability to automatically nullify extra or redundant sub-circuits—i.e., a form of model selection—via the sparsity regularization. The number of active connections and average connection strength (Fig. 7J, K) present modest variability between sub-circuits. This underscores the model’s ability to identify distinct sub-circuits with different structural behaviors, capturing various interaction resolutions. Specifically, it can find sub-circuits of varying strengths and is not limited to only k-sparse sub-circuits. From a neuroscience perspective, this suggests that underlying 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-region interactions) with either multiple ensembles from the same area affecting together ensembles from a different area (e.g., $\bm{f}_{3}$), or multiple ensembles from diverse areas converging to the same target area (e.g., $\bm{f}_{5}$). However, these source or target groups of ensembles vary across sub-circuits, hence the full repertoire of activity cannot be captured by a single. In addition to these between-area interactions, we also observed, although to a lesser degree, within-area interactions (e.g., as in $\bm{f}_{8}$ and $\bm{f}_{11}$, as marked by black arrows). These dynamics emphasize that within-area amplification or regulation dynamics do exist in addition to global cross-region computations, however seem to be less salient than cross-area communication.

When exploring the activation patterns of these sub-circuits (Fig. 3D, Fig. 17), we reveal variability between subjects and sessions. However, the usage of overall similar sub-circuits, with some sub-circuits not being used at all in any of the exemplary sessions (Fig. 3C, F, e.g., $\bm{f}_{12}$ not being used), while others are constantly used with high variability over time (e.g., $\bm{f}_{8}$ light-green). Interestingly, when exploring the within-session pairwise correlations between sub-circuits coefficients (Fig. 7I), it seems that most coefficients present very low ($<0.1$) correlations, indicating that the activity of the corresponding circuits is not similar and they might capture different cognitive processes, while other pairs of sub-circuits may show similarities in activity. This might indicate that for a specific subject or session, they may be related. For instance, a specific image may evoke some emotional response in one subject but not in another (e.g., $\bm{c}_{9}$ & $\bm{c}_{11}$ in the second session of subject 10, however not in the first session of that subject nor in the average correlations across all sessions).

## 5 Discussion.

In this paper we introduced CrEIMBO (Cross Ensemble Interactions in Multi-view Brain Observations), 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 leveraging joint inference and sparsity regularization to capture the underlying neural sub-circuits governing brain activity.

CrEIMBO offers several advantages over other models of whole-brain multi-session data. Chiefly, it identifies a cross-session shared latent space where non-stationary ensemble interactions are governed by a set of global sub-circuits.
By structuring the time-varying dynamics in terms of a basis set of sub-circuits, CrEIMBO allows for the exploration of active sub-circuits across sessions and aiding in the identification of subject-to-subject and cross-session variability in neural activity. Additionally, constraints on the loading matrices $\bm{A}$ enable the identification of area-specific ensembles, facilitating the understanding of within-area compared to between-area interactions.
Moreover, the simultaneous fitting of the ensembles and dynamics automatically aligns cross-session ensembles in terms of functionality, which is critical for studying cross-session variability.
We demonstrated the efficacy of CrEIMBO through both synthetic experiments and neural recordings. On synthetic data, CrEIMBO consistently recovered ground truth components with high accuracy. In our neural recordings analysis, CrEIMBO found cross-regional motifs that span interactions across sessions and subjects.

Limitations and future work:
While CrEIMBO leverages information from multiple brain sessions, the most modern Electrophysiology tools are still far from observing every neuron or region simultaneously. Consequently, and as with all dynamical systems models, obscured interactions with missing (i.e., unrecorded) neurons within the brain may limit the precision of the identified interactions.
Additionally, CrEIMBO’s current implementation does not include an offset term. Adding this would increase complexity and require thorough exploration; hence, we will address it in future work.
We limit CrEIMBO here to linear projections from the ensemble-space to per-session observation space. Extending this aspect of the model to more flexible, non-linear projections while maintaining interpretability is challenging and a future goal.

CrEIMBO’s effectiveness in aligning ensembles across sessions relies on the assumption that these sessions share commonalities in their active sub-circuits. If sessions operate with entirely different sub-circuits, cross-session ensemble alignment is not guaranteed.
Another limitation is that we assume Gaussian observation statistics, i.e., spiking data must be converted to rates to use CrEIMBO. Future research could incorporate non-Gaussian likelihoods, such as Poisson for spike counts.
Finally, CrEIMBO rests on the mathematical foundation of dictionary learning, which can be computationally intensive at higher dimensions. To enhance computational efficiency, we plan to implement parallel processing, data batching, and GPU processing in future iterations.

We also look to leverage CrEIMBO to identify pathological or cognitive states during which specific sub-circuits manifest exclusively (e.g., neural circuitry changes during epileptic attacks). Another extension for the current method includes extending the global cross-session circuits set to enable small, controlled structural adaptations of the sub-circuits between sessions, e.g., through an $\ell_{2}$ regularization on the distance between cross-session corresponding sub-circuits $\|\bm{f}_{k}^{d}-\bm{f}_{k}^{d^{\prime}}\|_{F}^{2}$.

## Acknowledgments

N.M. was supported by the Kavli Discovery Award of the Kavli Foundation.

## Appendix A Notations

| Symbol | Description |
| --- | --- |
| $D$ | Number of sessions (encompassing both same and different subjects) |
| $M_{d}$ | Number of trials for each session $d$ |
| $\bm{Y}_{m}^{d}$ | Neural recordings firing rate estimation for trial $m$ and session $d$ |
| $N_{m}^{d}$ | Number of neurons observed in trial $m$ of session $d$ |
| $\mathrm{T}_{m}^{d}$ | Number of time points observed in trial $m$ of session $d$ |
| $\{\bm{Y}_{m}^{d}\}_{m,d}$ | Overall set of observations |
| $\bm{A}_{m}^{d}$ | The neuronal composition of each sparse ensemble in trial $m$ and session $d$ |
| $\bm{A}_{\text{block}}$ | Block of $\bm{A}$ representing neural dynamics in a specific brain area |
| $n_{j}$ | Total neurons in area $j$ |
| $p_{j}$ | Maximum number of ensembles discernible within area $j$ |
| $p$ | Total number of ensembles across all areas |
| $J$ | Total number of distinct brain areas |
| $n$ | Total neurons involved across all trials or sessions |
| $\bm{x}_{t}$ | State vector at time $t$ |
| $\bm{F}_{t}$ | Transition matrix at time $t$, representing ensemble interactions |
| $K$ | Number of global interacting units |
| $\bm{c}_{kt}$ | Influence of the $k$-th global interaction at time $t$ |
| $\bm{f}_{k}$ | Basic linear system representing the $k$-th global interaction |
| $\tau=\mathrm{T}_{m}^{d}$ | Duration of the trial or session |

## Appendix B Differences from dLDS

While our proposed model draws inspiration from dLDS [], it significantly diverges in its approach, particularly in ensemble identification and the details thereof. Below we list the key differences:

- 1.
Incorporation of sparsity constraints on the operators $f$.
- 2.
Extension to multi-sessions through shared $f$s and distinct projections ($D$) instead of coefficients ($c$).
- 3.
Imposition of regional structure on $D$.
- 4.
Integration of graph-driven sparsity on $D$.
- 5.
Modifications to the inference procedure inspired by linocs.
- 6.
Ability to capture inter-subject variability.
- 7.
Adaptation to handle variations in trial duration, fluctuations in neuron count, and changes in distributions between sessions.

## Appendix C Differences from SiBBlInGS

Although our proposed model draws inspiration from SiBBlInGS [] for ensemble identification, it introduces several fundamental differences:

- 1.
SiBBlInGS does not consider interactions between ensembles, whereas our model incorporates ensemble interactions.
- 2.
SiBBlInGS does not address differences in the number of neurons between subjects, whereas our model accounts for this variability.
- 3.
SiBBlInGS lacks consideration for region-specific structure, focusing primarily on a single block of $D$, whereas our model imposes region-specific structure on $D$.
- 4.
SiBBlInGS does not account for the evolution of the ensemble activity space, whereas our model incorporates this aspect.

## Appendix D Parameters Range for Synthetic Data Experiment

The full set of parameters for the synthetic experiments are available in Table 1.

**Table 1: Parameter values for synthetic experiments**
| Parameter | Value (range) | Explanation |
| --- | --- | --- |
| same_c | False | Indicates if to use shared sub-circuits coefficients across sessions. |
| step_f | 0.1 - 0.5 | Gradient descent step size for updating $\bm{f}$ |
| GD_decay | 0.99 - 0.999 | Decay over iterations of gradient-descent step size. |
| max_error | 1e-09 | The error threshold to stop training the model. |
| max_iter | 500 | Threshold on the maximum number of iterations. |
| include_D | True | Indicates the inclusion of a projection to a latent space. |
| step_D | $10^{-5}-10^{-2}$ | Range of Gradient Descent Step size for $\bm{A}$. |
| seed | 0 - 4000 | Random seed to choose from |
| normalize_eig | True , False | Indicates if to normalize the sub-circuits by |
| start_sparse_c | False , True | Whether to initialize $c$ to be sparse. |
| sparse_f | True | Whether to apply sparsity on the sub-circuits |
| num_gradient_steps | 1 -3 | Number of gradient steps in an iteration. |
| add_avg | True , False | Whether to add a moving average to an iteration. |
| sparsity_on_f_max | 40-60 | Percentile of sparsity applied on each $\bm{f}$ |
| take_multiple_gd | False, True | Indicates if to take multiple steps of gradient descent. |
| D_graph_driven | True | Indicates that D is inferred with graph-driven way. |
| infer_x_c_together | False , True | Indicates if the inference of $\bm{x}$ and $\bm{c}$ occurred simultaneously. |
| include_mask | True | Indicates the inclusion of a mask on $\bm{A}$ |
| norm_D_cols | True | Indicates the normalization of columns of $\bm{A}$ |
| lambda_D | 0.1 - 0.3 | Regularization weight on updating $\bm{A}$ |
| step_D_decay | 0.99 - 0.9999 | Decay rate for updating $\bm{A}$ parameters |
| num_regions | 3 | Number of brain regions. |
| lambda_x | 0 - 0.1 | Sparsity term on the latent dynamics. |
| latent_dim | 3 | Latent dimension. |
| noise_level | 0.2-0.6 | Standard deviation of noise added to the data. |
| num_subdyns | 3 | Number of sub-circuits. |
| lasso_solver | ’spgl1’ | Pylops solver used for $\ell_{1}$ regularization. |

## Appendix E Parameters for human multi-regional experiment

The full set of parameters for the human data experiment are available in Table 2.

**Table 2: Parameter values for real-world human experiment**
| Parameter | Value | Explanation |
| --- | --- | --- |
| same_c | False | Indicates if sub-circuit coefficients are shared across sessions. |
| step_f | 0.2 | Step size for gradient descent in updating $\bm{f}$. |
| GD_decay | 0.992 | Gradient descent step-size decay rate over iterations. |
| max_error | $1\times 10^{-9}$ | Maximum allowable error. |
| max_iter | 100 | Maximum number of iterations. |
| seed | 0 | Random seed. |
| normalize_eig | True | Normalize eigenvalues. |
| start_sparse_c | False | Start with sparse $\bm{c}$. |
| include_D | True | Whether to include $\bm{A}$. |
| step_D | 0.0001 | Step size for updating $\bm{A}$. |
| num_gradient_steps | 4 | Number of gradient steps. |
| add_avg | True | Add average. |
| sparsity_on_f_max | 45 | Maximum sparsity on $\bm{f}$. |
| take_multiple_gd | False | Take multiple gradient descents. |
| D_graph_driven | True | Graph-driven $\bm{A}$. |
| infer_x_c_together | False | Infer $\bm{x}$ and $\bm{c}$ together. |
| include_mask | True | Include mask. |
| norm_D_cols | False | Normalize columns of $\bm{A}$. |
| lambda_D | 0.3 | Regularization parameter for $\bm{A}$. |
| step_D_decay | 0.9999 | Decay rate for step size of $\bm{A}$. |
| l1_D | 2.74 | $\ell_{1}$ regularization for $\bm{A}$. |

## Appendix F Information about the baselines:

- 1.
“Without Dynamic Prior”: Running CrEIMBO while removing the prior over the dynamics 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 on the dynamics.
- 2.
“All Regions (sparse)”: Running CrEIMBO without applying the multi-regional diagonal mask over the ensemble compositions ($\bm{A}$), thus supporting the finding of non-localized ensembles.
- 3.
“All Regions (sparse)”: Similar to 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 on a single session (view) of the data rather than leveraging cross-session information.
Since $\{\bm{f}_{k}\}$ is invariant to its ordering, we used SciPy’s implementation of the ‘linear_sum_assignment‘ problem [] to re-order the identified sub-circuits before comparing them with the ground truth, matching the $f$ values by minimizing the $\ell_{2}$ error.

## Appendix G Code and Data Sharing

The neural data used in this study [] is publicly available with License: “spdx:CC-BY-4.0”. We used version 0.240123.1806 of the data which can be accessed through the DANDI Archive at [https://dandiarchive.org/dandiset/000469/0.240123.1806](https://dandiarchive.org/dandiset/000469/0.240123.1806) (Access Information: dandi:OpenAccess).
The code will be shared via GitHub upon publication.

## Appendix H Technical Details

The whole code is written in Python. We utilized the following Python packages along with their respective versions: scipy (1.8.0), matplotlib (3.8.2), numpy (1.23.5), pandas (1.5.0), seaborn (0.11.2), scikit-learn (1.0.2), pylops (1.18.2), and statsmodels (0.13.2), all within a Python 3.10.4 environment. For the $\ell_{1}$ regularization, we used the “SPGL1” solver of pylops []. During the writing process, we used ChatGPT very sparingly to improve wording and fix grammar mistakes at the word or sentence level.
The synthetic experiment can run locally on a desktop (simple ASUS machine). However, the robustness tests with hundreds of random repeats and the human data experiments used the entry node with CPUs only of NERSC’s “Perlmutter” supercomputer. Running time for each synthetic instance ranges from 0.8 to 2 minutes (with the capability to run in parallel). For the human data, processing all 41 sessions with 1000 time points per session (and varying neuron counts as described in the text and figures) takes between 17 to 25 minutes, depending on parameters and result quality.

## Appendix I Societal Impact and Potential Harmful Consequences

With respect to Societal Impact, the development, application, and future work of CrEIMBO can advance our understanding of brain function. Future work can potentially lead to improved diagnostics of neurological disorders. We believe our method does not present any significant harmful consequences. We explicitly discussed the method’s limitations, ensured proper use of the data, and adhered to the NeurIPS code of ethics.

## Peer Reviews

confidence: 4
contribution: 3
initial_rating: 8
presentation: 3
questions: 1. Please provide details of how synthetic data was generated. 

2. 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?  
3. 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? 
4. 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.
soundness: 4
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.
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.
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.

---

confidence: 4
contribution: 2
initial_rating: 3
presentation: 4
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.
 
How difficult would it be to extend this model to use a Poisson observation model to better capture neural spiking activity?
soundness: 3
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.
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.
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.

---

confidence: 4
contribution: 3
initial_rating: 6
presentation: 3
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?
- 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?
- 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?
soundness: 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. 
- Authors performed exhaustive simulation experiments to show the importance of their model formulation. 
- The paper is well-written and easy to follow. The proposed model architecture and training framework are intuitive and effective.
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.
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. 
- 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. 
- 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.  
- 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.

Meta-review: The authors propose CREIMBO, a novel model that extracts interpretable neural subcircuits from multi-region, multi-session electrophysiological data, maintaining both compression and interpretability. CREIMBO demonstrates its effectiveness in both simulated and real data, revealing specialized subcircuits and interregional interactions, making it a promising tool for neuroscience applications, particularly for analyzing complex, high-dimensional neural recordings.

The major weakness lies in insufficient details regarding the simulation study, particularly the generation of synthetic data, which raised concerns about the validity of some findings. There were also queries about the robustness of CREIMBO and its comparison to real data and other baseline models. However, these concerns were largely addressed in the rebuttal, with the authors clarifying the synthetic data generation process and providing more context around model comparisons. Another suggestion is that the font size in the figures is excessively small, making them difficult to read. Additionally, the paragraphs are overly compact. The introduction to the approach feels too vague and lacks mathematical detail. The models should be presented with clear mathematical formulations rather than relying solely on textual descriptions. Please address these points as thoroughly as possible in the final revision.

The model's conceptual novelty and its ability to simultaneously compress data while preserving interpretability are major strengths, especially compared to deep learning approaches often used in neuroscience. CREIMBO was thoroughly validated through simulations and real data analysis, showcasing its potential in understanding multiregional neural dynamics. The paper is clear, well-written, and makes significant contributions to the field, demonstrating both technical rigor and scientific relevance.

Despite some weaknesses in the simulation study and comparisons with other models, the novelty, clarity, and demonstrated effectiveness of CREIMBO make it a valuable contribution to the neuroscience field, and I recommend it for acceptance.

Decision: Accept (Spotlight)