Accurate Identification of Communication Between Multiple Interacting Neural Populations
Authors: Belle Liu
Jacob Sacks
Matthew D. Golub
Venue: ICML 2025 poster (2025)
Keywords: Computational Neuroscience
Dynamical Systems
Interpretability
Brain-Wide Communication
Primary area: applications->neuroscience_cognitive_science

Abstract: Neural recording technologies now enable simultaneous recording of population activity across multiple brain regions, motivating the development of data-driven models of communication between recorded brain regions. Existing models can struggle to disentangle communication from the effects of unrecorded regions and local neural population dynamics. Here, we introduce Multi-Region Latent Factor Analysis via Dynamical Systems (MR-LFADS), a sequential variational autoencoder composed of  region-specific recurrent networks. MR-LFADS features structured information bottlenecks, data-constrained communication, and unsupervised inference of unobserved inputs--features that specifically support disentangling of inter-regional communication, inputs from unobserved regions, and local population dynamics. MR-LFADS outperforms existing approaches at identifying communication across dozens of simulations of task-trained multi-region networks. Applied to large-scale electrophysiology, MR-LFADS predicts brain-wide effects of circuit perturbations that were not seen during model fitting. These validations on synthetic and real neural data suggest that MR-LFADS could serve as a powerful tool for uncovering the principles of brain-wide information processing.

Accurate Identification of Communication Between Multiple Interacting Neural Populations
Belle Liu
Jacob Sacks
Matthew D. Golub
Abstract
et al., 2018; Stringer et al., 2019; Musall et al., 2019; Allen et al., 2019; Jia et al., 2022). Consequently, there has been growing interest in the design of data-driven communica- lion modely that seek infer the pathways and content of communication between the recorded regions.
Neural recording technologies now enable simul- tancous recording of population activity across many brain regions, motivating the development of data-driven models of communication between brain regions. However, existing models can struggle disentangle the sources that influence recorded neural populations, leading to inaccu- rate portraits of inter-regional communication. Here, we introduce Multi-Region Latent Factor Analysis via Dynamical Systems (MR-LFDDS). a sequential variational autoencoder designed to disentangle inter-regional communication, inputs from unobserved regions, and local neural pop- ulation dynamics. We show that MR-LFADS outperforms existing approaches at identifying communication across dozens of simulations of task-trained multi-region networks. When applied to large-scale electrophysiology, MR-LFADS pre- dicts brain- wine effects of circuit perturbations that were held out during model fitting. These val- idations synthetic and real neural data position MR-LFADS as promising tool for discovering principles of brain-wide information processing.
Accurately identifying inter-regional communication is chal- lenging for at least four reasons (Biswas et al., 2020; Kang & Druckmann, 2020; Keeley et al., 2020; Perich & Rajan, 2020; Semedo et al., 2020; Kass et al., 2023). First, com- munication signals are not directly observed in multi-region recordings. Although some recorded neurons might project to other recorded regions, the identity and targets of these projection neurons are typically unknown. Second, models may need to account for inputs to the recorded brain regions from other regions that were not recorded during the experi- ment. Third, models should faithfully reconstruct activity within each recorded region by accounting for communi- cation between recorded regions, inputs from unrecorded regions, and local neural population dynamics- capturing complex features such as structured trial-to-trial variabil- ity (Goris et al., 2014) and nonlinear, nonstationary, and state-dependent population dynamics (Shenoy et al., 2013; Vyas et al., 2020; Duncker & Sahani, 2021; Durstowitz et al., 2023). Fourth, accurate reconstruction of the recorded data does not guarantee accurate inference of the underlying communication. Many different models may sufficiently ex- plain the recorded data, leading to ambiguities about which, if any, should be trusted for scientific interpretation.
1. Introduction
Large-scale neural recording technologies, such as high- density electrophysiology (Steinmetz et al., 2019; Siegle et al., 2021; IBL et al., 2023; Chen et al., 2024; Bennett et al., 2024) and calcium imaging (Sofroniew et al., 2016; Song et al., 2017; Allen et al., 2017), now enable simulta- neous recording of neural population activity across many brain regions. These advances have revealed that many sensory, cognitive, and motor processes engage spatially distributed networks in the brain (Makino et al., 2017; Gilad
In this work, we introduce Multi-Region I Latent Factor Anal- ysis via Dynamical Systems (MR-LFADS), a multi-region communication model that directly addresses all of the chal- lenges outined above. MR-LFADS is a probabilistic model that represents each recorded region with a distinct set of stacked recurrent neural networks (RNNs) that capture the region's 's potentially nonlinear and nonstationary population dynamics. MR-LFADS represents communication between observed regions and inputs from unobserved regions as disentangled sets of latent variables. Structured informa- tion bottlenecks encourage the model to infer inputs from unobserved regions only when their effects cannot be plained by communication among the recorded regions. MR-LFADS infers single-trial initial conditions and time-
'Graduate Program in Neuroscience. University of Washing- ton; "Paul G. Allen School of Computer Science & Engineering, University of Washington. Correspondence to: Matthew Golub <mgolub@co weteeeecsssssstsn...omm washington.edu>.
Proceedings of the 4dnd International Conference Machine Learning, Vancouver, Canada. PMLR 267, 2025. Copyright 2025 by the author(s).

---

Accurate Identification of Communication Between Multiple Interacting Neural Populations
varying inputs that together account for trial-to-trial vari- ability in the recorded activity. This automatic inference of inputs eliminates the need to manually specify input signals, thereby avoiding strong, tfffiullt-o--vaidaat assumptions about how external signals influence each region. Finally, MR-LFADS constrains communication to originate from model-reconstruuted neural activity, rather than from more flexible latent representations design choice that, as we will show, enables more accurate inference of communica- tion without sacrificing the quality of data reconstruction.
two such techniques that, like MR-LFADS, are coupled nonlinear dynamical systems, with each dynamical system representing one recorded region: multi-population sticky recurrent SLDS (mp-srSLDS) (Glaser et al., 2020), and Multi-Region Switching Dynamical Systems (MR-SDS) KKaniol-Tambbuu et al., 2024). In mp-srSLDS each region is modeled by an SLDS, while MR-SDS each is modeled as a switching nonlinear dynamical system.
While these existing approaches to communication model- ing address some of the challenges outlined in Section 1, none, to our knowledge, address all four challenges. In particular, none of these existing methods support inferring inputs from unobserved brain regions. This functionality is crucial because inputs from unobserved regions might in- fluence the target region's 's population dynamics and modes of communication. Some approaches attempt to account for unobserved inputs by explicitly providing task-related signals as inputs to each region or by removing condition av- erages and modeling the residual single-trial neural activity. However, these manual strategies impose strong assump- tions about the content and targets of input signals. As we will show, misspecifying such inputs risks confounding in- ferred population dynamics and communication, leading to models that accurately reconstruct neural activity through incorrect mechanisms.
To evaluate MR-LFADS, we developed 37 synthetic multi- region datasets that capture real-world challenges in commu- nication modeling across a range of neuroscience-relevar scenarios. On these datasets, MR-LFADD consistently out- performs existing models in recovering the pathways and content of communication. Through targeted ablations of key design features, we demonstrate that these features indeed improve the identification of communication. We then applied MR-LFADS to multi-region electrophysiologi- cal recordings in mice performing a decision- meriiiiinmmninggggg task (Chen et al., 2024). In a subset of trials that were held out during model fitting, photoinhibition was applied to the an- terior lateral motor cortex. MR-LFADS predicted the brain- wide effects of these circuit perturbations, suggesting that MR-LFADS inferred an accurate account of inter-regional communication. Moreover, MR-LFADS infers consistent communication across multiple training runs from differ- ent random initializations, demonstrating its robustness and reliability in real-data settings.
To address such unobserved inputs, Pandarinath et al. (2018) introduced Latent Factor Analysis via Dynamical Systems (LFADS), a sequential variational autoencoder (sVAE) for modeling single-trial neural population dynamics within a single recorded brain region. LFADS jointly identifies non- linear dynamical system, implemented as an RNN, along with the single-trial initial conditions and time-varying un- observed inputs needed to drive the system to reconstruct single-trial neural population recordings. We henceforth use SR-LFADS to refer to a single-region LFADS model. MR- LFADS builds on this foundation to support multi-region modeling with explicit disentangling of communication, un- observed inputs, and local population dynamics.
2. Related Work
Existing communication models can be broadly categorized as either static or dynamic. Static methods predict each timestep of neural activity in a target region from a corre- sponding timestep of activity in one or more source regions and then interpret predictive source activity as inter-regional communication (Kaufman et al., 2014; Perich et al., 2018: 8; Ruff& & Cohen, 2019; Veuthey et al., 2020). Reduced-rank regression (RRR) is a prominent static technique that mod- els target-region activity as a low-rank linear function of the nource-region activity (Semedo et al., 2019; MacDowell et al., 2025). While static methods are straightforward to fit and interpret, they are typically limited capturing in- stantaneous, linear dependencies. Consequently, they do not readily account for nonlinear nonstationary relationships, or temporal structure that may arise due to neural population dynamics (Vyas et al., 2020).
3. Multi-Region LFADS (MR-LFADS)
MR-LFADS is composed of a set of SR-LFADS modules (Fig. Ia) that interact through constrained communication channels (Fig. lb). At a high level, MR-LFADS is cou- pled set of driven nonlinear dynamical systems that jointly trained to reconstruct all single trials of a multi- region dataset. Each recorded brain region i is modeled as a dynamical system that attempts to reconstruct the region-i recorded neural activity æf at each time = 1,..., T. Each region-i dynamical system evolves from a single-trial ini- tial state śó and is driven by (1) single-trial time-varying communication messages m| from other recorded brain regions j and (2) single-trial time-varying inferred inputs
Dynamic methods explicitly model temporal dependencies using, for example, switching linear dynamical systems (SLDS) (Linderman et al., 2016), RNNs (Perich et al., 2020), or Gaussian processes (Yu et al., 2008; Gokcen et al., 2022; 2024). We will pay particular attention to
2

---

Accurate Identification of Communication Between Multiple Interacting Neural Populations
(a)
SR-LFADS
(b)
MR-LFADS
Area
Encoders
my N-1
RNN
80
r
X
m1 1-2
m 2-11
u4
Area 2
1.1+100
m I-N
gn
Ponalty
KL Penalty
M
Prioe: (0,o'))
80 Me
Prior: (0(0,011))
: m 2-N
m N--2
KL. Penalty
RNN gen
Area N
Prior 11
8d
Reconstruction Loss
81
Figure 1. MR-LFADS architecture. (a) Single-region LFADS, as adapted for this work. (b) MR-LFADS with N regions. KL penalties from SR-LFADS in panel (a) are included in MR-LFADS, but are omitted in the diagram for elarity.
w/. representing input from unobserved brain regions.
rf WF( where r/ and Σ are the region-i predicted mean and co- variance of the time-l recorded neural activity, respectively, and are each computed via separate affine transformations, Wi and W For spike count observations, as modeled in the electrophysiology experiments of Section apply Poisson output distribution:
E diag (f)))
(5)
Notation. All time-indexed variables and parameters are also indexed by trial, though we omit trial indices for notational simplicity. We use 4/ to denote the inclu- sive sequence of integers f + 1,..., 12)- We use ('(r)==W'z W'x + b' to denote an affine transformation with weights W and offsets b².
Generative Model. MR-LFADS treats all initial states gó- communications m and inferred inputs of as latent vari- ables. The prior distributions over these latent variables are modeled as:
P(a .) = Poisson(r))
rf T exp (WF
(6)
where the exponential nonlinearity ensures non-negative predicted firing rates r/-
Po. ui, u4.mm ~ NN(00021
(1)
Inference Model. Following VAE conventions (Kingma & Welling, 2013), MR-LFADS approximates the intractable true posterior distributions over the latent variables using variational posteriors, denoted g(-l-).
Given these quantities, the neural population dynamics in region are modeled by a "generator" gated recurrent unit (GRU) network, GRU gen with internal states of that evolve
MR-LFADS defines the approximate posteriors over com- munication messages from observed regions j to as Gaus- sian distributions:
according to:
g:= = GRU
gen (gf- "
q(m) +i = = g(m) |
N (p m E!!!)
(7)
set of region-i factors fi are defined as an affine readout from the corresponding generator states:
with parameters derived from the region-j predicted firing rates r/:
f = - (ot)
1 = pm (rZ)
(8)
These factors are then transformed into parameters of time- varying output distributions a manner dependent on the nature of the neural recordings. For continuous-valued observations, as in calcium imaging, a Gaussian zero- inflated Gamma distribution may be appropriate (Zhu et al., 2022). In the synthetic-data experiments of Section 4, we apply a Gaussian output distribution:
M = diag exp (Wi
om
Constraining communication to be derived from rf anchors it to the neural recordings and in doing reduces ambiguity in system identification. We refer to this rate-based commu- nication model as MR-LFADS(R) In Section 4, we also explore generator- based MR-LFADS(G) and factor-based MR-LFADS(F)) communication models, which replace all instances of ri in Eq. 8 with gf and f(. respectively.
P(# gb, uhe {m 1iz1) = N(rr, E/
(4)
3

---

Accurate Identification of Communication Between Multiple Interacting Neural Populations
Approximate posteriors over the region-i initial generator states gó and inferred inputs (from unobserved brain regions) nn are defined as the following Gaussian distributions:
4. Results I: Synthetic Multi-Region Datasets
Here, we evaluate MR-LFADS' ability to recover the ground truth inter-regional communication across broad range of synthetic multi-region datasets. Each dataset was generated by a unique data-generan geaaeettngninggggg network (DGN): an ensemble of noisy RNN modules jointly trained to perform specified cognitive neuroscience task, with each module representing a distinct brain region. Prior to training, we explicitly specified the presence or absence of directed, low- rank communication channels between each region pair.
d | ===+0)=== N )ne q(u, | x = M ((uu, I,,))
(9)
(10)
where the mean and covariance parameters are computed from the corresponding conditioning recorded neural ac- tivity of via a set of region-i-specific "encoder" and "con- troller" GRU networks (see Appendix A.I). Our approach here slightly modifies the original SR-LFADS specification, which allows acausal inference via a bidirectional encoder network that processes each entire T-timestep neural record- ing F1:T to infer 9o and each element of "1:T. In contrast, we infer & (Eq. 9) using bidirectional encoder applied only to past neural activity x'-,.0- preserving causality, and we infer un (Eq. 10) using a unidirectional encoder RNN that processes sist in strictly forward, causal manner. This formulation ensures that all predicted firing rates rt. and thus all derived communication signals m² are inferred causally from neural activity recorded up to time L.
In Experiments and 2, we manually designed the DGNs to impose specific challenges outlined in Section 1. In Experiment 3, we randomly generated dozens of DGNs, each trained to perform a randomly selected cognitive neu- roscience task. Across all experiments. we treated each module's 's hidden-unit activity as recorded neural activity and retained all external inputs, inter-module connectivity, and communication signals as ground truth. These ground truth quantities are crucial for evaluating communication models but are typically not directly observable in real-data settings.
Model Fitting, Following VAE conventions, MR-LFADS is trained by maximizing the evidence lower bound (ELBO), a variational lower bound on the data bogllikelhhood. The ELBO a sum of two terms: (1) the expected log-likelihood
To benchmark MR-LFADS, we compare it to three estab- lished communication modeling techniques- RRRR (see Ap- pendix A.2), mp-srSLDS, and MR-SDS (see Section 2). To demonstrate the importance of specific MR-LFADS design features, we also compare against ablated MR-LFADS vari- ants that selectively exclude those features. We focus eval- uation on each method's 's ability to recover a causal model of each DGN CGGNiincludng both the pathways and content of inter-regional communication-givenn a multi-region dataset generated by the DGN.
T Eg log P(1x)) I (ot), (m) [D]
(II)
and (2) the negative Kullback-Leibler (KL) divergence DKL between the approximate posteriors (Eqs. 7-10) and the priors (Eq. 1) over the latent variables.
To assess communication pathways, consider recovery of "effectome" (Pospisil et al., 2024) deseribing the causal flow of effects along the inter-regional connectome. We rep- resent this effectome as a matrix with each element (i, (izz) indicating the volume of directed communication flow from region j to i. The effectome reflects both the inter-regional connectivity and the magnitude of communication flow over each directed connection. For each dataset, we compare model-inferred effectomes to the ground truth effectome by computing cosine similarity Sc.. between the vectorized effectome matrices.
The expected log- leeellkei measures reconstruction ac- curacy and is estimated by running samples from the ap- proximate posteriors (Eqs. 7-10) through the generative model to evaluate the experiment-dependent output distri- butions from Eqs. 4-6. The DKL, term acts as a regular- izing information bottleneck on the latent variables. To control this regularization, we allow rescaling of the DKL. term (Higgins et al., 2017; Keshtkaran et al., 2022). Noting that the Dxi. term decomposes into contributions from the three sets of MR-LFADS latent variables {gb}, { and {mm we weight each contribution differently and treat the weights (5p, Bu., Bm) as hyperparameters. To encour- age MR-LFADS to infer inputs from unobserved regions only when that information cannot be obtained as communi- cation from observed region, we propose structured KL bottleneck with B, 10,,.. Other choices of KL regular- ization structure might be appropriate if priori knowledge is available about information flow or anatomical connectiv- ity between recorded regions. See Appendix A.I for further detail on MR-LFFADS.
To assess communication content, we quantify how well the model-inferred messages capture the information in the ground truth messages. Specifically, we apply lin- ear regression predict the ground truth messages m! from the inferred messages μ (Eq. 7), and we re- port cross- valosd-v coefficient of determination, denoted R2 (u/ht, , m For MR- MR-LFADS, we use an analogous procedure to compare model-inferred inputs yhu, (Eq. 10) to ground truth external inputs. This comparison is not applica-
4

---

Accurate Identification of Communication Between Multiple Interacting Neural Populations
ble RRR, mp-srSLDS, or MR-SDS, as these methods do not infer inputs from unobserved regions. Additional details
sages. In the ground truth DGN at time 4, region-i receives only of and m An accurate communication model should therefore infer inputs and messages that encode only these time-t quantities. However, due to structure of the memory task, the DGN's region-i time-l activity contains information about s/-4:t and m = 5 6tt-2 2. To re-
on experiments, model hyperparameters, and evaluation metrics are provided in Appendix B-D.
4.1. Experiment I: Inferring Unobserved Inputs
This experiment demonstrates that inferring inputs, rather than manually specifying inputs, enables more accurate identification of inter-regional communication. To eval- uate design implications related to imput specification, we designed a DGN that implements a dynamical memory func- tion. Each region of this "memory network" receives unique stimulus information from one observed region and one observed region, and is tasked with remembering recent history of those signals (Fig. 2a, left). With each region of the DGN receiving information from both observed and unobserved sources, this setup poses the challenge of disen- tangling whether each signal arises due to communication or due to external imput. A common modeling choice is to provide all known external inputs to all model regions and to let model fitting determine which inputs are needed by each region. However, in this case, such manual input specification can result in a model completely forgoing com- munication (Fig. 2a, right) because the manually specified inputs contain the information that was actually transmitted as communication in the DGN.
(a)
Memory Network
Ground Truth
Incorrect
f- 1-2
S
comumnncction
2
A
A
A
A
readout
Si 1-2
S-ate s
Sia 1-2
Model Fit
(c)
Inferred Effectome
8
(x.r)
Srow
MR LEADS (R)
MR-L-FAS (5)
M9 505
mo
MR-SDS
RRR
We evaluate MR-LFADS(R), which does not use the stim- ulus signals during training evaluation but rather auto- matically infers external inputs in an unsupervised manner. We also evaluate an ablated MR-LFADS variant that does not automatically infer inputs. Termed MR-LFADS(S), this model receives all stimulus signals as manually specified external inputs to each region's generator GRU gen (replac- ing uf with s! 1:3 in Eq. 2). By comparison, MR-SDS and mp-srSLDS also receive all stimulus signals as manually specified inputs, and RRR neither receives nor infers exter- nal inputs.
(d)
Inferred Inputs
(e) Inferred Communication
No lag
2-step lag
2-step lag
4-step
Inferred Input Mio
Mis
Target Repion
A2
Rin 9)
R2
A
A2
A1
A
Ground Truth Signal
Source Region
Figure 2. Experiment 1. (a) Left: DGN setup. Each region, area A', receives private stimulus of (red) and communicates a two-step delayed version st (blue) to downstream region. Each region is trained to recall (green) the last five time steps of its private stimulus and its received communication. Right: Potential incor- rect communication model capable of accurately reconstructing these synthetic neural data. (b) H² scores for data reconstruction. Box plots describe distributions values across 10 models fit from distinct random initializations (seeds). Boxes represent the interquartile range (IQR), and whiskers extend to the most extreme points within 1.5 IQR from the quartiles. (c) Left: Cosine similar- ity between inferred effectomes and ground truth. Right: Example fitted models. Color intensity indicates the relative message norm, computed by concatenating all multidimensional messages across trials and time, taking the 2-norm, then normalizing across com- munication channels. 2 indicates the ground truth input. (d) n² of linear prediction of ground truth stimulus inputs (with time lag E {0, 2}) from inferred inputs. (c) R2 of lincar prediction of ground truth messages (with time lag E {2, 4}) from inferred messages.
The MR-LFADS variants (R and S) reconstructed the simu- lated multi-region activity more accurately than MR-SDS, mp-rSLDS,, and RRR (Fig. 2b, Fig. Sla). Critically, MR-LFADS(R) accurately infers the effectome (Fig. 2c, left, Fig. Slb). By contrast, all models that are manually provided stimulus signals (MR-LFADS(S). MR-SDS, mp- iSLDS) infer less accurate effectomes and demonstrate the failure mode mentioned above, forgoing communication and instead relying on the specified inputs to provide the corresponding signals. RRR also infers a less accurate effectome. These results suggest that manually specify- ing inputs can discourage models from utilizing-and thus By automatically inferring inputs, MR-LFADS(R) avoids this failure mode.
Next, we assess the accuracy of inferred inputs and mes-
5

---

Accurate Identification of Communication Between Multiple Interacting Neural Populations
construct the data, a communication model must account for how this time-lagged information becomes represented in region-i at time-t. MR-LFADS(R) correctly infers the current timestep ground truth stimuli s as inputs to each region (Fig. 2d, left), while correctly avoiding inferring time- lagged versions of those stimuli (Fig. 2d, right), despite their utility for data reconstruction. Similarly, MR-LFADS(R) accurately recovers the current timestep messages m. - (Fig. 2e, left) and avoids incorrectly inferring time-lagged versions (Fig. 2e, right). By contrast, all other models infer no communication or communication with incorrect tempo- ral lags (Fig. Slc). See Appendix B.2 for further details.
ular, accurately identifying and localizing the ground truth pass-through and integration computations. Though the inte- gration dynamics are localized to A in the DGN, an overly flexible model might instead learn integration in A,, yielding accurate data reconstruction but mislocalizing the computa- tion, e.g., if A" integrates 5 into de and communicates d: as mn P-D to AP (Fig. 3a, right).
On this pass-decision dataset, all MR-LFADS models achieved comparably high-fidelity data reconstruction, slightly outperforming MR-SDS and mp-srSLDS (Fig. 3b, Fig. S2a). RRR reconstructed the data poorly, likely due to its inability to capture the timescale of the integration computation, i.e., de cannot be predicted from any single- timestep value Spr-
Taken together, these results indicate that MR-LFADS(R) learned region-specifie population dynamics consistent with those that implement the memory functions in the ground truth DGN. Moreover, these results demonstrate the unique ability of MR-LFADS(R) to disentangle region-specific ex- ternal inputs, communication between recorded regions, and local population dynamics, all in an unsupervised man- ner that mitigates biases associated with manual specifi- cation of external inputs. See Appendix 1 for further analyses linking this disentangling (Miller et al., 2024) to MR-LFADS(R)'s structured information bottlenecks.
MR-LFADS(R) inferred the most accurate effectome (Fig. 3c, Fig. S2b), identifying A" communication but not AD-AP. AP. By contrast, all other models, includ- ing MR-LFADS(F) and MR-LFADS(G), identified spurious
Pass-Decision Network
Ground Truth
Incorrect
A"
D sign(d)
D
d
4.2. Experiment 2: Cata-Constrained Communication
5
Se
This experiment demonstrates the implications of message- inference design choices. We designed MR-LFADS(R) to infer messages affine functions of the source-region pre- dicted firing rates rt (Eqs. 7-8), which are tied to the ob- served source-region neural activity of through the data reconstruction term in the ELBO (Eq. 11). This data- constrained communication architecture contrasts with that of MR-SDS and mp-srSLDS, which infer communication as a function of less-constrained, source-region latent dynami- cal states (see Table S3). To directly evaluate the implica- tions of this design choice, all within the MR-LFADS frame- work, we designed model variants with costtrnntrr factor-based and generator-based communication, termed MR- (R--FFAD((F) and MR-LFADS(G), respectively.
(b)
Model (c)
Inferred Effectome
MR. LFADS
R r)
S
d f
(R)
A
A
0
MR-- LFADS
A
HR L
Mi- (AOREE)
90'
mop MoorSSLDS DS
888
(G)
We 1000
A^
AP
SDS
(d) Inferred Inputs (e)
Inferred Communication
Ri .)
(me
To highlight the significance of this design choice, we evalu- ate models data from two-region "pass-decision" DGN that computes perceptual decisions based on time-varying sensory evidence (Fig. 3a, left). An upstream region, area A", receives a white noise stimulus s, and is trained to pass that stimulus through to a readout, effectively learning an identity function routing input to output. A downstream decision region, area A", receives this routed stimulus as communication m. P-D integrates that stimulus over time into decision variable de, and reports # choices indicating the sign of that decision variable (Mante et al., 2013).
MR sos000
Mn-DD
mp Sanss
PRLLFAOSIN
MR 2e022000 muss Seeett sua Re6
MI-- FRDSII
MR F
MIK Mn 305
mip- SLDS
RPR
WiLL
mp
Figure 3. Experiment 2. (a) Left: DGN setup, with stimulus (red), communication (blue), trained readouts (green). and computations (I: identity; f:: integration). Right: Potential failure mode for a learned model. (b) R² scores for data reconstruction. (c) Left: Cosine similarity between inferred effectomes and ground truth. Right: Example fitted models. Color intensity indicates relative message norm. (d) R² of linear prediction of ground truth input to region P, from inferred inputs. (e) Left: R2 when predicting ground truth messages m P-D =s from inferred messages, um P.D Right: R² when predicting the decision variable d from pm indicating mislocalization of integration.
This setup again challenges models to disentangle external inputs, communication, and local dynamics- and in partic-
P-DD
6

---

Accurate Identification of Communication Between Multiple Interacting Neural Populations
AP communication.
MR-LFADS(R)'s data-constrained communication, which improves identification of communication by reducing am- biguity inherent to overly flexible models. See Appendix E.2 for further analyses into this excessive flexibility.
MR-LFADS variants correctly inferred inputs to A that encoded s, (Fig. 3d). Because MR-SDS and mp-srSLDS do not infer unsupervised inputs, we estimated their effective inputs by passing the manually specified inputs through their corresponding trained input mappings. In MR-SDS, these effective inputs to A 4 carried markedly less information about 8, relative to MR-LFADS and mp-srSLDS.
4.3. Experiment 3: Generalization Across Random Regti- Networks
The previous experiments utilized datasets synthesized by DGNs designed to highlight specifie failure modes of com- munication models. However, real brain-wide networks exhibit a broad range of architectures and computations. To assess generalization across such a broad range of settings, here we evaluate communication models on datasets gener- ated by a wide variety of randomly configured multi-region DGNs (Fig, 4a, top), each trained to perform a randomly selected cognitive neuroscience task (Fig. 4a, bottom). We generated 35 multi-region datasets, each from a unique task- trained DGN consisting of three or four regions and random- ized inter-regional connectivity. Tasks were drawn from the set described by Yang et al. (2019), spanning multiple vari- ants of decision-making, working memory, categorization, and inhibitory control (see Appendix D).
Accurate identification of communication requires inferred messages Df to encode the stimulus 8t- By contrast, mf instead encoding the decision variable de would imply mislocalization of the integration dynamic (Fig. 3a, right). Only MR-LFADS(R)), mp-srSLDS, and RRR cor- rectly encoded s, in m m (Fig. 3e, left) without incor- rectly encoding de (Fig. 3e, right).
Taken together, these Experiment 2 results again demon- strate MR-LFADS(R) outperforming existing methods at disentangling external inputs, communication, and lo- cal population dynamics, without sacrificing data re- construction. The comparisons to MR-LFADS(F)) and MR-LFADS(G) specifically highlight the importance of
(a) Random Multi-Region
(b)
Model Fit
Inferred Effectome rel. MR-LFADS (R)
(d)
Inferred Communication rel. MR-LFADS (R)
Task- Trainer Networks
Best
Network Configuration
1
0.25
18
Noisy Stimules Random Number of
R-
0.00
L m)
0.0
AR Pa)
0.0
Noisy Fizarion
Ragions
-0.5
-0.55
-05
Worse than
Randor
MRR Soss0
MR 5D5
mp-orsips
ARR
MR. LFADSS MACLLEEA FADSUE NR- MIRS miL SO00
MR-LFADS (R)
MR LEADSIF!
MiK 1000
MADsi
MM 120 (I8SS
min ALD
RAR
Raidoo Connectivity
Cog-Neuro
Task
mip-iS
NIR
Exumple Task:
Multisensory Decision-No Making (Raposo et al., 2014)
(e)
Inferred Effectome
Ground Truth
MR-LFADS (R)
MR-LFADS (F)
MR-LFADS (G)
MR-SDS
mp-ort.SDS SDS
Fixation
Saccade
Exsule
A
Regast
Stimulus
Response
(1, oily
Exsrlee
Choose larger
Figure 4. Experiment 3. (a) Top: Each DGN is configured with random number of areas (3 or 4), random inter-regional connectivity, and is trained on a randomly selected task. Bottom: Example task setup. On each trial, the DGNs received a noisy fixation stimulus $ns,, and two noisy task stimuli, represented in polar coordinates as : (a("), o(*)) for & E 21. DGNs were trained to process these inputs into task-dependent outputs: response angle (reP and task-dependent saccadic eye movement. Here, multisensory decision- maiisinnnnmaiiggiii is depicted as an example task. When the fixation cue disappears, the output area A must saccade and report the g( value corresponding to the stimulus with larger (9) (b) R² scores for data reconstruction. Box plots describe distributions of values across 35 DGN fits. (c) Cosine similarity (S.w) of inferred effectomes relative to ground truth, compared to that of MR-LLADS(RR) with LFDSSO) One-tailed t-test p-values: p 0.00016, 0.12, 0.006, 0.0, 0.01. NS; not significant. (d) Left: R² scores, A.S... Sos cos
relative to MR-LFADS(R), for linear prediction of ground truth messages from inferred messages. One-tailed t-test p-values: = 0.0, 0.0. 0.001, 0.0, 0.1 Right: n² scores, relative to MR-LFADS(R), for linear prediction of inferred messages from ground truth messages. One-tailed t-test p-values: p = 0.0004, 0.0, 0.08, 0.0, 0.0. (c) Inferred effectomes from models fit to datascts from a three-region DGN (top) and four-region DGN (bottom). We chose the datasets on which MR-LFADS(R) achieved its median Scss scores (across the 35 generated datasets). Color intensity indicates relative message norm
7

---

Accurate Identification of Communication Between Multiple Interacting Neural Populations
We fit MR-LFADS(R), (F), and (G), along with MR-SDS, mp-srSI.DS, and RRR, to each of these datasets. Aggre- gating results across all datasets, the MR-LFADS mod- els achieved the best data reconstruction, which was in- distinguishable across model variants (Fig. 4b, Fig. S3). MR-LFADS(R) and MR-LFADS(G) inferred the most ac- curate effectomes, with statistically indistinguishable Sans distributions (Fig. 4c). To evaluate the aceuracy of in- ferred message content, we attempted to linearly decode the ground truth messages from the inferred messages and quantified accuracy using the R² scores of these predictions (Fig. 4d, left). We also performed the reverse, predicting the inferred messages from the ground truth and interpreted lower R2 scores as an indication that inferred messages contained additional information beyond that present in the ground truth (Fig. 4d, right). MR-LFADS(R) was the only model that performed best across both of these metrics. In- ferred effectomes from two example datasets are shown in Fig. 4e. Taken together, these results demonstrate that MR- (F-DSR))(R) outperforms existing communication mod- els across a broad range of neuroscience--ellvan synthetic multi-region datasets.
were held out from training. We also assess the reliability of inferred inter-regional communication across random ini- tializations, comparing MR-LFADS(R) to MR-LFADS(G)..
We trained a 5-region MR-LFADS(R) model on simulta- neously recorded population activity from anterior lateral motor cortex (ALM), midbrain reticular nucleus (MRN), superior colliculus (SC), thalamic regions with strong, recip- rocal connections to ALM (Thal(A)), and other thalamic re- gions (Thal In a subset of trials excluded from training, ALM was transiently photoinhibited (Fig. 5a). MR-LFADS was trained only on unperturbed ("control"") trials, with pho- toinhibition trials held out for validation. See Appendix F for further detail.
We first verified that MR-LFADS accurately reconstructed multi-regional neural activity in held-out control trials (Fig. 5b). We then adapted MR-LFADS(R) to predict the effects of ALM photoinhibition on the remaining recorded regions (Fig. 5c). To mimic ALM photoinhibition in silico, we ablated MR-LFADS communication from ALM to the other regions by zeroing all outgoing messages from ALM. We summarized the temporal influence of ALM photoin- hibition by computing the differences between condition- averaged population activity in photoinhibition (F)) ) and control (ri) trials. MR-LFADS predicted these photoinhibi- tion effects (Fig. 5d), despite never seeing photoinhibition trials during training. Namely, MR-LFADS predicted that Thal(A) would be most affected by ALM photoinhibition, MRN least affected, and SC and Thal(O) intermediately af- fected. These results demonstrate that MR-LFADS learned
5. Results II: Multi-Region Electrophysiology
Here, we apply MR- (R-LFFAD(R) to large-scale electrophysi- ological data from multiple simultaneously recorded Neu- ropixel probes in mice performing a decision-making task (Fig, 5a) (Chen et al., 2024). We test whether MR-LFADS can predict the effects of causal circuit perturbations that
(a)
Experimental Task (Chen et al., Cell. 2024)
(b)
Example Single-Tree Mriill Neural Activity
(c)
Brdin-wide Effects of ALM Photoinhibition
ALM Cakes
Thalamus (A) Da
Thalamos (O) Cui
ATM
Thalamos (A)
Thalamos (O)
f
Photoinhibition applied ALM
h High
Law
Fintng Rone 200 Itt
MRN
11 spikes
MRN Evars
SC 1
Inter- Trial
Sample
Delay
Go
Response
Joaphnacc Mle
spike count
smmotthed
Ruse
P-inh: Control: 7
Interval
0.5 s
LRADS
Time
(d)
Predicted Photoinhibition Effects
(c)
Effectome
(1)
MR-LFADS(R) Inf. Comm. Across Seeds
(2)
Inf. Comm. Correlation
Left- Righe Iit
1 65 (R))
Data
MR-LFADS
11 (a.u)
Seed 155
Seed 2
3/5
95
X" P' (
Thation MINN
0.595
08-
0.992
02-
4.1
0.0
0.1
0.7
0.1
0.0
0.1
0.2
Time (3)
ThakAA
MEBUN
MER
ALM
Thal/A)
ALM
Time p-inh
MINN
ThaiAs
ThalAA)
Figure 5. MR-LLADDS(R) applied to multi-region, high densiyy electrophysiology. (a) Mice receive high- or low-tone auditory stimulus ("sample") and respond by licking left or right ("response"". (b) MR- LFADS single-trial predicted firing rates in held-out control trials (blue), recorded spike times (black vertical ticks). and smoothed, binned spike counts (black; causal, exponential filter) for example neurons in each modeled brain region. (c) Condition-averaged smoothed spike counts of example neurons in control and photoinhibition trials (left-hit condition). Shaded regions indicate the standard deviation across trials. (d) Photoinhibition-rclaa -elaaatttiio--laeeedd changes in population recordings, as observed experimentally (left) and as predicted by MR-LFADS (right). (c) Cosine similarity of inferred effectomes across models with different random initializations (seeds). (I) Message norms inferred by MR-LFADS(R) for all connections across example seeds k. (g) Correlation of inferred message norms across pairs of seeds (k, l). Box plots describe distributions of values across 5 models fit from distinct random initializations.
8

---

Accurate Identification of Communication Between Multiple Interacting Neural Populations
model of inter-regional communication that is accurate enough predict multi-regional effects of causal circuit perturbations.
MR-LFADS models were also more reproducible across random model initializations compared to a model variant that removed data-constraints on inferred communication.
Finally, we evaluated the consistency of MR-LFADS across random initializations of the model parameters. We com- pared MR-LFADS(R) and MR-LFADS(G),, each trained us- ing five different random seeds, on simultaneous population recordings from ALM, thalamus, and MRN (see Appendix F2). Both models inferred consistent effectomes (Fig. 5e), but message content was more consistent across seeds in MR-LFADS(R) (Fig. 5f, g), further highlighting the ben- efits of data-constrained communication- particulary for improving the reproducibility of scientific conclusions de- rived from the model.
Despite its advantages, MR-LFADS has potential limita- tions. MR-LFADS can be sensitive to hyperparameter (HP) settings, mirroring a known limitation of SR-LFADS (Keshtkaran et al., 2022). Thus, significant computational resources might be required to adequately optimize HPs. As in SR- SR-LFADS, HPs specifying the prior distributions (Eq. 1) can shape the representations of inferred latent variables, which can in turn shape the inferred population dynamics. Although inferred inputs have been shown to reflect the pres- ence, identity, and timing of external inputs (Pandarinath et al., 2018), future work is needed to interpret the represen- tations of those signals (Sedler et al., 2023; Versteeg et al., 2024). The data-constrained communication propose in this work might mitigate this representational sensitivity in the context of inferred communication.
6. Discussion
Understanding how brain regions interact to support dis- tributed computation requires communication models that can disentangle inter-regional communication from local population dynamics, accounting for region- region-specific inputs from unrecorded regions. In this work, we identified critical failure modes that limit existing communication models- including misidentification due to manually specified ex- ternal inputs and mislocalization of neural dynamics due to overly flexible communication architectures. We intro- duced MR-LFADS, a communication model specifically designed to mitigate such failures through three key de- sign features: (1) automatic inference of eggion-specffic inputs from unobserved sources, (2) data-constrained com- munication inferred from reconstructed firing rates, and (3) structured regularization that promotes disentangling and prevents inferring inputs from unobserved sources when the same information can be obtained via communication from observed regions. These features discourage the model from learning spurious solutions that explain the data but misrepresent inter-regional interactions. While MR-LFADS is one conerete implementation, the three design principles are more general, and alternative architectures- for exam- ple, different formulations of the local-region dynamical mystems-may also succeed if they incorporate these same principles.
Another set of MR-LFADS HPs (the B's from Section 3) control the model's preference to infer communication rather than inputs from unobserved regions whenever pos- sible. While the settings we chose were effective in this study, different choices may be needed in other scenarios especially when an unrecorded region directly drives multi- ple recorded regions. This classical correlation ccassation confound cannot be resolved from passively observed neu- ral activity alone. While MR-LFADS recovers more causal structure than existing methods, its greatest value may be in generating hypotheses about communication and motivating targeted causal circuit perturbations to test them.
A Acknowledgments
This work was supported by NIH award T32-MH132518 (IS), the Paul G. Allen Foundation (MDG), and NIH award R00-MH121533 (MDG). We thank O. for valuable discussion and support with MR-SDS (Karniol- Tambour et al., 2024). We are grateful to S. Chen, N. Stein- metz, F. Shea-Brown, and J.N. Kutz for insightful input on the project, and to H. Gurnani, T. Kim, and J. Pemberton for thoughtful feedback on the manuseript. We are also grateful to D. Sussillo, K. Shenoy, W. Newsome, and C. Pandarinath for many years of guidance and inspiration relevant to this project.
Using synthetic datasets designed to rigorously test com- munication models, we demonstrated that MR-LFADS out- performs existing approaches in accurately recovering both the structure and content of inter-regional communication. Crucially, ablated model variants indicated that these perfor- mance gains stem directly from MR-LFADS design features. Applying MR-LFADS to real multi-region electrophysio- logical recordings further validated its utility. MR-LFADS inferred inter-regional interactions that accurately predicted brain-wide effects of causal perturbations, despite these perturbations being absent during training, In this setting,
Impact Statement
This work aims to advance the fields of machine learning and neuroscience, with potential societal impacts includ- ing the development of neuroengineering technologies and treatments for neurological injuries, diseases, and neuropsy- chiatric conditions.
9

---

Accurate Identification of Communication Between Multiple Interacting Neural Populations
References
Glaser, J., Whiteway, M., Cunningham, J. P., Paninski, L., and Linderman, S. Recurrent switching dynamical sys- tems models for multiple interacting neural populations. Advances in Neural Information Processing Systems, 33: 1 14867-14878, 2020.
Akiba, T., Sano, S., Yanase, T., Ohta, T., and Koyama, M. Optuna: A next- next-generation hyperparameter optimization framework. In Proceedings of the 25th ACM SIGKDD International Conference on Knowledge Discovery and Data Mining, 2019.
Gokcen, E., Jasper, A. I., Semedo, J. D., Zandvakili, A.,
Kohn, A., Machens, C. K., and Yu, B. M. Disentan- gling the flow of signals between populations of neurons. Nature Computational Science, 2(8):512-523 525, 2022.
Allen, W. E., Kauvar, I. V., Chen, M. Z., Richman, E. B., Yang, S. J., Chan, K., Gradinaru, V., Deverman, B. F.,, Luo, L., and Deisseroth, K. Global representations of goal-directed behavior in distinct cell types of mouse neocortex. Neuron, 94(4)88919907, 2017.
Gokcen, E., Jasper, A., Xu, A., Kohn, A., Machens, C. K., and Yu, B. M. Uncovering motifs of concurrett signaling across multiple neuronal populations. Advances in Neural Information Processing Systems, 36, 2024.
Allen, W. E., Chen, M. Z., Pichamoorthy, N., Tien, R. H., Pa- chitariu, M., Luo, L., and Deisseroth, K. Thirst regulates motivated behavior through modulation of brainwide neu- ral population dynamics. Science, 364(6443))ccaav33322 2019.
Goris, R. L.., Movshon, J. A., and Simoncelli, E. P. Parti- tioning neuronal variability. Nature Neuroscience, 17(6): 358-865, 2014.
Higgins, I., Matthey, L., Pal, A., Burgess, C., Glorot, X., Botvinick, M., Mohamed, S., and Lerchner, A. beta- vae: Learning basic visual concepts with a constrained variational framework. In International Conference on Learning Representations, 2017.
Bennett, C., Ouellette, B., Ramirez, T. K., Cahoon, A., Cabasco, H., Browning, Y., Lakunina, A., Lynch, G. F,, McBride, E. G., Belski, H., et al. Shield: Skull-shaped hemispheric implants enabling large-scale electrophys- iology datasets in the mouse brain. Neuron, 112(17): 2869-2885, 2024.
IBL, Benson, B., Benson, J., Birman, D., Bonacchi, N., Bougrova, K., Bruijns, S. A., Carandini, M., Catarino, J. A., Chapuis, G. A., et al. brain-wide map of neu- ral activity during complex behaviour. bioRxiv preprint. 2023.
Biswas, T., Bishop, W. F., and Fitzgerald, J. E. Theoreti- cal principles for illuminating sensorimotor processing with wide neuronal recordings. Current Opinion in Neurobiology, 65:138-145, 2020.
Jia, X., Siegle, J. H., Durand, S., Heller, G., Ramirez, T. K., Koch, C., and Olsen, S. Multi-regional module-based signal transmission in mouse visual cortex. Neuron, 110 19855-1115 2022.
Chen, S., Liu, Y., Wang. Z. A., Colonell, J., Liu, L. D., Hou, H., Tien, N.-W., Wang, T., Harris, T., Druckmann, S., et al. Brain-wide neural activity underlying memory- guided movement. Cell, 187(3):676-99.. 2024.
Kang. B. and Druckmann, S. Approaches to inferring multi-regional interactions from simultaneous popula- tion recordings. Current Opinion in Neurobiology, 65: 108- 2020.
Dai, B., Wang, Y., Aston, J., Hua, G., and Wipf, D. Con- nections with robust pea and the role of emergent sparsity in variational autoencoder models. Journal of Machine Learning Research, 19(1)):1-42 2018
Karniol- Tamibuur, O., Zoltowski. D. M.. Diamanti, E. M.. Pinto, L., Brody, C. D., Tank, D. W., and Pillow, J. W. Modeling state-dependent communication between brain regions with switching nonlinear dynamical systems. In The Twelfth International Conference on Learning Repre- sentations, 2024.
Duncker, L. and Sahani, M. Dynamics the manifold: Identifying computational dynamical activity from neural population recordings. Current Opinion in Neurobiology, 70:163-170, 2021.
Durstewitz, D., Koppe, G., and Thurm, M. I. Reconstructing computational system dynamics from neural data with recurrent neural networks. Nature Reviews Neuroscience, 24(11):693-710, 2023.
Kass, R. E., Bong, H., Olarinre, M., Xin, Q., and Urban, K. N. Identification of interacting neural populations: methods and statistical considerations. Journal of Neuro- physiology. 130(3):755-96,, 2023.
Gilad, A., Gallero-Sale Salas, Y., Groos, D., and Helmchen, F. Behavioral strategy determines frontal or posterior location of short-term memory in neocortex. Neuron, 99 (4):814 8288, 2018.
Kaufman, M. T., Churchland, M. M., Ryu, S. I., and Shenoy, K. V. Cortical activity in the null space: permitting prepa- ration without movement. Nature Neuroscience, 17(3): 440 440. 2014.
10

---

Accurate Identification of Communication Between Multiple Interacting Neural Populations
Keeley, S. L., Zoltowski, D. M., Aoi, M. C., and Pillow, J. W. Modeling statistical dependencies in multi-region spike train data. Current Opinion in Neurobiology, 65: 202, 2020.
Perich, M. G., Gallego, J. A., and Miller, L. E. A neural population mechanism for rapid learning. Neuron, 100 (4):964-976, 2018.
Perich, M. G., Arlt, C., Soares, S., Young, M. E., Mosher, C. P., Minxha, J., Carter, E., Rutishauser, U., Rudebeck, P. H., Harvey, C. D., al. Inferring brain-wide inter- actions using data Cataaconstraiied recurrent neural network models. bioRxiv preprint, 2020.
Keshtkaran, M. R., Sedler, A. R.., Chowdhury, R. H., Tandon, R., Basrai, D., Nguyen, S. I.., Sohn, H., Jazayeri, M., Miller, L. E., and Pandarinath, C. A large-scale neural network training framework for generalized estimation of single-trial population dynamics. Nature Methods, 19 (12):1572-1577, 2022.
Pospisil, D. A., Aragon, M. J., Dorkenwald, S., Matsliah, A., Sterling, A. R., Schlegel, P., Yu, S.-c., McKellar, C. E., Costa, M., Eichler, K., al. The fly connectome reveals path to the effectome. Nature, 634((0032:201-209, 2024.
Kingma, D. P. and Welling, M. Auto-encoding variational bayes. arXiv preprint, 1312.61 14, 2013.
Liaw, R., Liang, E., Nishihara, R., Moritz, P., Gonzalez, E., and Stoica, I. Tune: A research platform for dis- tributed model selection and training. arXiv preprint arXiv: 1807.05 / 18, 2018.
Raposo, D., Kaufman, M. T., and Churchland, A. category-free neural population supports evolving de- mands during decision-making. Nature neuroscience, 17(12):184411792 2014.
Linderman, S. W., Miller, A. C., Adams, R. P., Blei, M., Paninski, L., and Johnson, M. J. Recurrent switching linear dynamical systems. arXiv preprint, 1610.08466, 2016.
Ruff, D. A. and Cohen, M. R. Simultaneous multi-area recordings suggest that attention improves performance by reshaping stimulus representations. Nature Neuro- science, 22(10):6669-1670 2019.
MacDowell, C. J., Libby, A., Jahn, C. I., Tafazoli, S., Ardalan, A., and Buschman, T. J. Multiplexed subspaces route neural activity across brain-wide networks. Nature Communications, 16(1):3359, 2025.
Sedler, A. R., Versteeg, C., and Pandarinath, C. Expressive architectures enhance interpretability of dynamics-basse neural population models. Neurons, Behavior, Data Anal- ysis, and Theory, 2023, 2023.
Makino, H., Ren, C., Liu, H., Kim, A. N., Kondapaneni, N., Liu, X., Kuzum, D... and Komiyama, T. Transformation of cortex-wide emergent properties during motor learning. Neuron, 94(4):880-88) 890, 2017.
Semedo, J. D., Zandvakili, A., Machens, C. K., Byron, M. Y., and Kohn, A. Cortical areas interact through a communication subspace. Neuron, 102(():249-22)) 2019.
Mante, V., Sussillo, D., Shenoy, K. V., and Newsome, W. T. Context-dependent computation by recurrent dynamies in prefrontal cortex. Nature, 503(7774):78-844 2013.
Semedo, J. D., Gokcen, E., Machens, C. K., Kohn, A., and Byron, M. Y. Statistical methods for dissecting interac- tions between brain areas. Current Opinion in Neurobiol- ogy, 65:59 69,999999990 2020.
Miller, K., Eckstein, M., Botvinick, M., and Kurth-Nelson, Z. Cognitive model discovery via disentangled rnns. Ad- vances in Neural Information Processing Systems, 36, 2024.
Shenoy, K. V., Sahani, M., and Churchland, M. M. Cortical control of arm movements: a dynamical systems perspec- tive. Annual Review of Neuroscience, 36(1) 36(1):337-359 2013.
Musall, S., Kaufman, M. T., Juavinett, A. L., Gluf, S., and Churchland, A. K. Single-trial neural dynamics are domi- nated by richly varied movements. Nature Neuroscience, 22(10):1677-1666, 2019.
Siegle, J. H., Jia, X., Durand, S., Gale, S., Bennett, C., Graddis, N., Heller, G., Ramirez, T. K., Choi, H., Luviano, J. A., et al. Survey of spiking in the mouse visual system reveals functional hierarchy, Nature, 592(7852):86-92. 2021.
Pandarinath, C., Oheaa J., Collins, J., Jozefowicz, R., Stavisky, S. D., Kao, J. C., Trautmann, E. M., Kaufman, T., Ryu, S. I., Hochberg, R., al. Inferring single- trial neural population dynamics using sequential auto- encoders. Nature Methods, 15110)805-815, 2018.
Sofroniew, N. J., Flickinger, D., King, J., and Svoboda, K. A large field of view two-photon mesoscope with sub- cellular resolution for in vivo imaging. eLife, 5:e14472, 2016.
Perich, M. G. and Rajan, K. Rethinking brain-wide interac- tions through multi-region 'network of networks' models. Current Opinion in Neurobiology, 65:146-151, 2020.

---

Accurate Identification of Communication Between Multiple Interacting Neural Populations
Song, A., Charles, A. S., Koay, S. A., Gauthier, J. L., Thiberge, S. Y., Pillow, J. W., and Tank, D. W. Volu- metric two-photon imaging of neurons using stereoscopy (vTwINS). Nature Methods, 14(4):420 426, 2017.
Steinmetz, N. A., Zatka-Haas, P., Carandini, M., and Harris, K. D. Distributed coding of choice, action and engage- ment across the mouse brain. Nature, 57677786)666))7 2019.
Stringer, C., Pachitariu, M., Steinmetz, N., Reddy, C. B., Carandini, M., and Harris, K. D. Spontaneous behaviors drive multidimensional, brainwide activity. Science, 364 :eav7)3,, 2019.
Versteeg, C., Sedler, A. R., McCart, J. D., and Pandarinath, C. Expressive dynamics models with nonlinear injective readouts enable reliable recovery of latent features from neural activity. In Proceedings of the 2nd NeurlPS Work- shop on Symmetry and Geometry in Neural Representa- tions, volume 228 Proceedings of Machine Learning Research, pp. 258. PMLR, 2024.
Veuthey, T., Derosier, K., Kondapavulur, S., and Ganguly, K. Single-trial cross-area neural population dynamics during long-term skill learning. Nature Communications, 11(1):4057, 2020.
Vyas, S., Golub, M. D., Sussillo, D., and Shenoy, K. V. Computation through neural population dynamics. An- nual Review of Neuroscience, 43(1) (3(11:249-222) 2:249-275, 2020.
Watanabe, S. Tree-structured parzen estimator: Understand- ing its algorithm components and their roles for better empirical performance. arXiv preprint arXiv:2304.11127 2023.
Yang, G. R., Joglekar, M. R., Song, H. F., Newsome, W. T., and Wang, X.-J. Task representations in neural networks trained to perform many cognitive tasks. Nature Neuro- science, 22(2):297-306, 2019.
Yu, B. M., Cunningham, P., Santhanam, G., Ryu, S., Shenoy, K. V. and Sahani, M. Gaussian-proceer factor analysis for low-dimensional single-trial analysis of neu- ral population activity. Advances in Neural Information Processing Systems, 21, 2008
Zhu, F., Grier, H. A., Tandon, R.. Cai, C., Agarwal, A., Giovannucci, A., Kaufman, M. T., and Pandarinath, C. A deep learning framework for inference of single-trial neural population dynamics from calcium imaging with subframe temporal resolution. Nature Neuroscience, 25 (12):1724-1734, 2022.
12

---

Accurate Identification of Communication Between Multiple Interacting Neural Populations
A. Models
A.1. Multi-Region LFADS (MR-LFADS)
Reconstruction validation. During training, we held out 10% of the neurons in each real or synthetic brain region for validation. After training MR-LFADS, we fit a separate linear decoder post hoc on the training trials, regressing the inferred factors onto the held-out neurons' actívity. Since the decoder is trained end-to-end with MR-LFADS, held-out neurons do not influence model fitting. Goodness of fit is reported as the R between predicted and observed held-out-neuron activity on validation trials.
KL penalty. A critical hyperparameter during training is the scale of the KL penalty. The KL penalty coefficient for the inferred inputs, Su. is always set higher than that for communication, 3m. an implicit assumption that encourages the network to prioritize learning information via communication channels whenever possible. The timing of when KL penalties are introduced also impacts results, though to lesser extent. The schedule we found to work well begins with an initial stage where no KL penalty is applied, allowing the model to overfit to the data. Next, the penalty for inferred inputs is introduced, discouraging the model's reliance on these inputs. Finally, the penalty for communication is added, limiting the model from learning excessive information through communication channels.
Weight regularization. We apply light L2 regularization to all GRU network recurrent weights.
Hyperparameters. Table S1 summarizes key hyperparameters used in the MR-LFADS models. Overall, we find that KI. coefficients have the most impact on held-out neuron loss, San, and R² scores for messages compared to other hyperparameters. SR-LFADS was originally described with factor layer that is potentially lower dimensional than the number of generator units or modeled neurons. Here, remove the rank constraint from the generator hidden states rates by setting Nice Nucu- Additionally, the inferred input and message channel dimensions only need to exceed the estimated true dimensionality of these quantities, as KL penalties naturally suppress redundant channels by driving their activity to zero, as discussed in Appendix E.I.
Table S1. Key hyperparameters for MR-LFADE models. Experiment refers to applications to multi-region electrophysiology data.
Hyperparameter
value
Description
learning rate
E [10-5, 0.004]
Scheduled by PyTorch's ReduceLROnPlatee
initial value: 0.004 Total time steps used for inferring inferred inputs
T
190
10
Total time steps used for inferring the initial condition
total epoch
350
Total number of epochs
3u
KL penalty coefficient for u: performs search for this hyperparameter
3,, start epoch
50
Epoch which 8,, starts increasing from
3n increase epoch
200
Number of epochs for 3u, to reach the maximum value
3m
KL penalty coefficient for m; performs search for this hyperparameter
8m start epoch
150
Epoch at which 3. starts increasing from 0
3m increase epoch
100
Number of epochs for 8m to reach the maximum value
8gn
8,
KL penalty coefficient for In
104
L2 penalty coefficient
a start epoch
0
Epoch at which a starts increasing from 0 Number of epochs for a to reach the maximum value Number of neurons, (64, 16, 64) for exp 1, 2 and respectively
a increase epoch
80
Ni
neu
Nr
2N
Generator size
yen
Nou
Ni
Neu
Factor size
Ni
Inferred input dimension, (4,8, 6, 8) for exp 1, 2, and 4 respectively
N
Inferred message dimension, (4, 16, 10, 8) for exp 1, 2, 3 and 4 respectively
msg
13

---

Accurate Identification of Communication Between Multiple Interacting Neural Populations
Unidirectional Encoder and Controller. To ensure that inferred inputs reflect only causal information-so that messages, in turn, causal and thereby allow mechanistic interpretation in MR-LFADS- modify the original LFADS model so that both the encoder and controller used for input inference are entirely unidirectional:
? GRU ene, 15
(12)
C GRU con (c--1 (2
The inferred inputs are then given by:
a( xin) - q(ua) I (t) = N (Khe,e Ei,1) M Wi (c) Ey diag( exp (Wi (c))
(13)
By contrast, the encoder for the initial condition can remain bidirectional, since it operates on data preceding = 1 and thus does not violate causality:
e = GRU
cac, (e 1(
(14)
e, enc.g (e 22,)
GRU
( | = 4(96 |6= = =N (1ig D
13 = W Moo (ee">:
(15)
Σx = diag exp (Wa (6e-- (b+)))
A.2. Reduced-Raane Regression
The RRR model used in this study is based on the inter-regional communication subspace model of MacDowell et al. (2025), which integrates reduced-rank regression with ridge regression. The key difference our implementation that we apply the rank constraint separately each source brain region, allowing us to disentangle the contributions of individual areas. Specifically, rather than concatenating activity from all regions into single input matrix governed by shared rank-constrained weight matrix which conflates signals across regions- we assign dedicated weight submatrix each source region A3, each with its own rank constraint r'.
Therefore, for the communication subspace model, we have:
W arg min | -W [[]] mnk(W)-r
which is equivalent to:
Wriidge arg min ||Y - XWI + a ||
W arg min Y Xwridgo W
(16)
XW XW|
which is then equivalent to:
W Wridgge V, V.V.T. where UEVT X Widdge
where X E T represents the activity of all source regions concatenated together, and Y E Tx represents the activity of the target region. Wridgg the weight matrix obtained after applying ridge regression, and a is the ridge regularization parameter. WTT is the reduced-rant regression matrix obtained after ridge regression is applied. In the final step, singular value decomposition (SVD) is applied W where UEVT is the decomposition, and V,, corresponds to the top r components of V.
In this version, the Wridgg matrix is divided into chunks corresponding to different regions A
W
1 1rdige. W Ni
(17)
14

---

Accurate Identification of Communication Between Multiple Interacting Neural Populations
where each Wridge.J corresponds to the contribution of source region A', and [W ;..;; ...] represents a vertical stack of the matrices. SVD is then applied to each individual Wridge matrix:
U wridge.s.
(18)
The reduced-rank version of Wridge. is computed as:
W vydigo.
pd (V
(19)
Finally, all W'n. matrices are concatenated form the complete nn matrix. This ensures that the rank reduction applied to each Wridge. only compresses the information within that specific region's contribution, preserving the interpretability of the communication pathway from A to the target region.
The hyperparameters for this model include a and a matrix R E RNXN, where each element r² represents the rank associated with the communication from source region A to target region A'. Additionally, since time delays may exist between regions- and and such delays are explicitly configured the synthetic datasets- a delay parameter d' is introduced for each nouree-target communication channel.
For fitting the memory and pass-decision networks, tune the model, we perform iterative search based cross- validation performance, optimizing the hyperparameters in the following order: D {d' : = ,,..., N, * j}, a and R = {1 L, = 1...., N, i + j}. This process is repeated until the hyperparameter values converge. The final values used are provided in Table S2. While the rank values R for both networks did not converge exactly to the true ranks of the messages, they were close. Notably, providing the true number of latents did not necessarily lead to better results. The delay values D were accurately learned for both models.
For fitting the networks in Experiment 3, the true delay is directly provided, and other hyperparameters are iterated in the same order for 10 epochs.
To increase the robustness of the RRR model fit, we implemented bagging approach. For each model, 10 trials were bootstrapped from the training set, with each trial containing 200 time steps. A total of 87 fitted models were averaged to obtain the matrix W This specific number of models was chosen to ensure that the total number trials used during training remained consistent with other models.
Table S2. Hyperparameter search results for RRR models.
Memory
-essiseciiion
O
0.055
0.01
12
24
R.
12
32
24
18
d3 i Fj
2
O
A.3. Multi-Region Switching Dynamical Systems
We consider two variants of multi-region switching dynamical system models. The first is mp-srSLDS (Glaser et al., 2020), which consists of linear transitions, dynamics, and emissions. The relevant hyperparameters are the number of latent states per region, the number of discrete switching states, amount of J,, and L2 regularization on the weights, and the learning rate. Additionally, we consider MR-SDS (Karniol- Tambour et al., 2024), which is an extension that uses nonlinear transitions, dynamics, and emissions. It consists of two components: an inference network and latent state-space model. The inference network is a transformer that performs the amortized inference of latent variables given observed neural activity. The latent state-space model is composed of number of networks and functions as structured prior on the latent variables. Specifically, we consider additive communication and input terms to the latent dynamics of the state-space model. That is, messages from other regions and external inputs affected the latent dynamics via additive terms. Relevant hyperparameters include the number of latent states per region, the number of discrete switching states, and the sizes of each sub-network.
For both models, we did extensive hyperparameter tuning to find the best model for each of the synthetic datasets and then computed all metries on held-out test set. We used the Tree-structured Parzen Estimator (TPH) algorithm (Watanabe, 2023)
15

---

Accurate Identification of Communication Between Multiple Interacting Neural Populations
with the Optuna backend (Akiba al., 2019) in Ray Tune (Liaw al., 2018). The algorithm fits two Gaussian mixture models (GMMs), one to the set of parameter values associated with the best objective and another the remaining ones. It chooses new parameters to explore by maximizing the ratio of the likelihood between these two GMMs. As such, it is search strategy which uses results from prior tested hyperparameters to inform the next choice of hyperparameters to test. We used the TPE algorithm to search over all relevant hyperparameters above. Additionally, for MR-SDS, we used dropout for regularization with the default settings in the provided implementation. We also manually picked a good learning rate and number epochs for training. Finally, we also made use of co-smoothing for evaluation. This holds out set of neurons from the inference network, and computes the fit on the reconstruction of these held-out neurons.
A.4. Model Comparisons
We outline the key design features of MR-LFADS variants and existing communication models in Table S3. Models are compared across four criteria: (1) eegion-specffic dynamics, (2) unsupervised inferred inputs, (3) data-constrained communication, and (4) structured information bottlenecks. Only MR-LFADS(R) incorporates all four features.
MR-LFADS(S)) ablates the controller, removing inferred inputs and instead using manually specified external inputs for each region. MR-LFADS(F) and MR-LFADS(G) both communicate via latent variables not directly grounded in observed data-using factors and generator states, respectively.
MR-SDS and mp-srSLDS include region-specific dynamics but rely on external inputs and latent-variable-basedd messaging, without any regularization on inputs or communication. RRR infers communication from observable quantities but lacks dynamics and inputs altogether. While it enforces rank constraint the communication subspace, this is not equivalent to explicit regularization on messages.
Table S3. Design features of MR-LFADS variants and existing communication models.
Unsupervised Inferred Inputs
Cata-Constrained Structured Information Communication
Dynamics
Bottlenecks
MR-LFADS(R)
MR-LFADS(S)
x
MR- MR-LFADS(F)
X
MR-LFADS(G)
V
V
X
MR-SDS
V
X
X
X
mp-srSLDS
V
X
X
X
RRR
X
X
X
B. Evaluation Metrics
B.1. Quantifying Effectome Similarity
In a trained MR-LEADS model, we define the inferred effectome to be a matrix of pairwise message norms, M, with element M., as the average value of || 2|2 (Bq. 7) across all trials and timesteps. In Experiments 1-2, we compared model-inferred effectomes to the corresponding ground truth connectivity matrix Mine. consisting of ones and zeros to indicate the presence absence of a communication channel in the DGN, respectively. In contrast, Experiment features networks in which not all connections are actively used; in this case, define Miue analogously M, but computed using ground truth messages m instead of inferred messages /. To assess similarity between M and Mine, we flatten these matrices into vectors-nñ and mugue and compute their cosine similarity:
(m, Maue)
Scos
E [0, 1],
(20)
la
Manoll l|2
where perfect alignment is indicated when Scos = 1.
To visualize an inferred effectome (c.g., Fig. 2c, right), we plot arrows whose color intensity reflects the relative message norm, computed by concatenating all multi-dimensional messages across trial and time, taking the 2-norm, and then
16

---

Accurate Identification of Communication Between Multiple Interacting Neural Populations
normalizing across all channels, with inferred communication and inputs normalized separately. The arrows in Fig. 4e, bottom further scaled by a sigmoid function for visual contrast. All reported values elsewhere are computed without applying any thresholds scaling, and heatmap visualizations of the effectome are also provided in Fig. Slb, Fig, S2b.
B.2. Evaluating Information Encoded in Learned Messages
In Experiment we tested whether the inferred inputs and communications encoded information about the past ground truth values, the network's s hidden units activity contains information about past inputs. A correct model should only learn the ground truth inputs or communications, For the results shown in Fig. 2d-e, right, the r-squared values were calculated with a time lag d as 122 3- m m))d)
C. Synthetic Datasets for System Identification Issues
Synthetic datasets for networks from Experiments 1-3 have 1024, 1024 and 820 total trials respectively, of which 85% is used for training, and 15% for validation. The length of each trial is 200 time steps.
C.1. Memory Network
In this synthetic network, each region A is modeled GRU network with 64 units that receives a private stimulus sf (0, L) with dimensions r'. To simulate communication channels carrying different amounts of independent information, the dimensions are set as (r' ,2, r²) (2, 3, 4). Each region has linear readout W and the outputs are required to encode information about the history of all inputs (i.e., stimuli and communication) for up to 5 time steps, enforced using a mean squared error loss. Similarly, the messages transmitted between regions are trained to match s and are also optimized via a mean squared error loss. Additionally, each region is subjected to dynamic noise & ~ ((0 0.01 introduced as perturbations to the RNN activity each time step.
After training the synthetic network, for MR-LFADS, we performed hyperparameter sweep over the KL. penalty coefficients for inferred inputs (B, E {0.01, 0.1, 1, 103) and communication (3m. E {0.001, 0.01)). The coefficient pair that resulted in the lowest held-out neuron loss, (Bw, fm) = (0.1, 0.01), was selected. Using these optimized coefficients, we ran the MR-LFADS fit across 10 different seeds, which randomizes model initialization and subsequent sampling of inferred quantities during training, but does not change the allocation of training versus validation data.
Comparing Experiments and 2, it is shown that mp-srSIDDS and MR-SDS activity reconstruction performance underper- forms compared to the MR-LFADS variants in Experiment (Fig. 2b), but not in Experiment 2 (Fig. 3b). One possible explanation for this discrepancy is the amount of information that the latent variables must encode each time step. For example, in area A the private stimulus is s E IR2 and the incoming message from area A is m = Sy_2 E R.A. As result, the latent representation at time t must capture information spanning 5 time steps and 6 variables in total. Under a standard hyperparameter tuning scheme (Section A.3), latent variable models like mp-srSLDS and MR-SDS may struggle to represent all this information accurately.
C.2.
Network
In this synthetic network, each area is an GRU that contains 16 units. The stimulus s, is two-dimensional and independently sampled from an exponential distribution with rate parameter of 3 time steps. Using a non- non-Gaussian distribution tests whether the model can learn the correct solution despite structural mismatches between the data and the model, as both the prior and posterior of inferred inputs and messages are Gaussian.
Each region has linear readout Wouts whose output is required to match its corresponding latent variables (s, for area AP and d, for P). The message sent from A to A is trained to represent 8, Additionally, A must encode whether d, is greater or less than 0 at all times, mimicking binary decision-making process.
For the dessioni network, a low KL penalty for inferred inputs (B,, = 0.0075) was necessary to achieve good held-out neuron loss, while the KL penalty for communication (8m = 0.001) was set to be approximately one order of magnitude smaller. Using these coefficients, we ran the MR-LFADS fit across 10 different seeds.
17

---

Accurate Identification of Communication Between Multiple Interacting Neural Populations
Activity Reconstruction
MR-LFADS (R)
MR-LFADS (S)
MR-SDS
mp-srSLDS
RRR stipeee
Model
(b)
Inferred Effectome
Norm. 11 11 Norm. |ellls
A²
A
0
Region (
1.00
Tagget Rogion!
A-
A3
A3
Sourcee Region )
Inferred Communication
2-step lag
4-step lag
2-stop lag.
4-stop lag
2-step lag
4-step lag
2-step lag
4-stop lag
2-step lag.
4-step lag
Target Ragiont
R
4
A3 Region )
A
Figure SI. Memory network. (a) Example traces of data reconstruction. Row: example neuron from area, column: different trials. (b) Inferred effectomes visualized as heatmaps. Top: Inferred inputs, bottom: communication. Color intensity represents message norms normalized by the largest input (max, [| message norm (max,, 1| (2) within cach model. These largest inferred inputs are, from left to right: 76, 175, 60, 188, O. The largest messages are: 195, 22, 95, 1431, 5, (c), R of lincar prediction of ground truth messages (with time step lag the left, 4 time step lag the right) via inferred messages. MR-LFADS(R) results (left) replicated from Fig. 2e.
18

---

Accurate Identification of Communication Between Multiple Interacting Neural Populations
(a)
Activity Reconstruction
MR- LRADS (R)
MR-LFADS (F)
MR-LFADS (G)
A"
P MMA
Myp mmA
Awme
nwm
mmA
mwm
12
MR-SDS
mp-srSLDS
RRR
AP
A
Data
Model
(b)
Inferred Effectome
MR-LFADS
MR-LFADS (F)
MR-LFADS (G)
MR-SDS
mp-srSLDS
RRR
A
AP
4
4
4
Nomn || lb
A
A"
A
Aa
AD
A
A
Region
1.00
Target Region
0.75
Norm. 1%
0.50
0.25
0.00
46
A6
A0
46
A"
40
Source Region
j
Figure S2. Pass-decision network. (a) Example traces of data reconstruction. Row: example neuron from an area, column: different
trials. (b) Inferred effectomes visualized circuit diagrams and heatmaps. Top: circuit diagram, middle: inferred inputs, bottom: communication. Color intensity represents message norms normalized by the largest input or message norm within each model.
D. Randomly Generated Multi-Region Networks
We generated a distribution of networks designed to perform computational tasks inspired by Yang et al. (2019). Each network consists of either 4 regions. connection probability (0.5, 0.6, 0.7) specified, and the connectome is randomly drawn. To be considered valid, the connectome must meet two criteria: (1) each region must have at least one input connection and one output connection to ensure region is redundant, and (2) all regions must be within maximum distance of 2 steps from the output region. Once valid connectome is generated, the network is trained on one of the computational tasks. During training, dynamic noise & ~ (0, 0.01 T) is applied to all regions, and the only loss based on whether the output region produces the correct response. Networks that meet the performance thresholds (train accuracy 0.8 and validation accuracy > 0.6) are selected as synthetic datasets.
For this case, since aim collect distribution of results, we do not perform hyperparameter sweep the KI. penalties. Instead, we fit one instance of each communication model to each synthetic dataset using single random seed. The computational tasks inspired by Yang al. (2019) are described below.
All trials are 200 time steps in length. Each task receives a fixation input 8ix.r, stimuli from two channels s = (a 000) and o (a (()), and requires saccade response r-acc and an additional response Tresp.t (1, Oees)) which differ based on the specified task. For different tasks, stimuli may come from one both of the channels. Both the stimuli and responses are expressed in polar coordinates, with a resolution of 10 degrees per angle.
The tasks are described in terms different families: the go task family, context-dependen decision decision-making family, and matching family. For all tasks, faa E [30, 50) denotes the onset of the first (and sometimes only) stimulus (Table S4).
19

---

Accurate Identification of Communication Between Multiple Interacting Neural Populations
Table $4. Task parameters for all familics in multi-region data generating networks.
Task
Doffset
Aeeay
sace
grep
Go
N/A
0
Esta + Adur
((1)
Anti-Go
N/A
Estatt + Adur
(())
Delay-Go
N/A
[30, 50)
Istan + Adur + Adelay
0(x)
Delay-Anti-Go
N/A
[30, 50)
Istan + Aur
Dielay
(())
DM1
0
[30, 50)
Isaa + Adur
(
DM2
0
[30, 50)
Estant + Adur
(
MultSen DM
0
[30, 50)
Estast + Adur
0
arg max, r(()
Delay-DNl DM1
[10. 20)
[30, 50)
Estart + Adur + Doftset
(())
Delay-DMe DD2
[10. 20)
[30, 50)
Estart + Adur + Doffsct
((2)
Delay MultSen DM
[10, 20)
[30, 50)
Isturt + Adur + Doffset
0 i arg max; r(()
Angle
[10, 20)
[30, 50)
Estart + Adur + Doffer if ((t) (12)
(
Anti-Angle
[10, 20)
[30, 50)
Latan + Adur + Dotker if (tL) (112)
π + ()22)
Category
[10, 20)
[30, 50)
Ismn + Aur Doffset if sign (0t)) sign(0(2))
(f2)
Anti-Category
[10, 20)
[30, 50)
start + Adtr Dofisen if sign( (0()) sign(0(2)
π + ((2)
The duration of all stimulus pulses in a trial is represented by Aeur E [30, 50), and Aotiset specifies the offset between the two stimuli, applicable. The time between the last stimulus offset and the fixation cue offset is given by Adelay- where deeayy indicates that the fixation cue never disappears. The onset of the saccade is denoted as f Each parameter-taart Adur Doffiet and dis drawn independently and uniformly from its specified half-open interval (21.
D.1. Go Task Family
The common characteristic of tasks in this family is that only one of the stimulus channels contains the signal, which varies between trials. Depending on the specific task, the network must saccade dependently or independently of the fixation cue. The response is required to be either in the direction of the signal pulse, (($), or in the opposite direction, æ + 0(0) The individual tasks are summarized in Table $4.
D.2. Context-Dependent Decision-Making Family
For this family of tasks, stimulus pulses occur in both channels, and the network must report either ((1) or ((22, depending on the specific task type. In some tasks, the pulses occur at different times, requiring the network to maintain memory of the stimuli. The task parameters for this family are summarized Table $4.
D.3. Matching Family
In the matching tasks, the network determines whether to saccade based on whether the two stimulus angles "match." the "Angle" tasks, the network saccades only if ((1) (12) under the given resolution (10 degrees per angle). the "Category" tasks, the network saccades if sign(0")) = ( meaning both angles are either positive or negative, Task details are provided in Table S4. Regardless of whether the angles match, the response is always set report the angle ((2) (or the opposite of it).
20

---

Accurate Identification of Communication Between Multiple Interacting Neural Populations
(a)
Activity Reconstruction (3 Region)
(b)
Activity Reconstruction (4 Region)
MR-LFADS
MR-SDS
MI-LFADS (R)
MR-SDS
O)-LIADD
SSLDS
MR-IFADS
MR-LFADS (33)
ARR
FR-IFADS
RRR.
Data
Model
Model
Figure Randomly generated multi-region network: example traces of data reconstruction networks with median MR-LFADS(R) performance for 3 regions (a) and 4 regions (b).
21

---

Accurate Identification of Communication Between Multiple Interacting Neural Populations
E. Implications of Constrained Architectural Choices in Restricting Message Content
E.1. Input and Message Inference with KL Penalties
KL penalties with standard Gaussian priors in variational autoencoders are known to reduce latent space dimensionality by pruning unnecessary dimensions (Dai et al., 2018; Miller al., 2024). Consequently, with sufficiently high KL penalty coefficients (B., 9,), MR-LFADS is incentivized
only the communication channels essential for data reconstruction. This effect is evident when comparing the most active channel (i.e., the one with the highest input message across trials and time) to the least active ones (Fig. S4a, b). Examining input and message norms across all channels further confirms that some channels are effectively silenced (Fig. S4c).
(a) Inferred Inputs yu
Communication pm
(c)
Channel norm for p Mm
40
Least - Most active channels
[8]
1
20:
1
2
3
100
75
50
25
0
k
Figure $4. The effect of KI. penalties restricting information across inferred inputs and communication channels in Experiment 1. (a) Inferred imput norm over time, ||p for the most (dark red) and least (light red) active channels, where & denotes channel number. (b) Inferred message norm over time, ||100) la, for the most (dark blue) and least (light bluc) active channels. (c) Inferred input (top) and message (hottom) norm across channels for area A in descending order of its scalar norm across trials and time.
E.2. Message Inference via Rates versus Factors
Since MR-LFADS factors and reconstructed rates share the same dimensionality and are related by linear transformation, it may not be immediately clear how message inference from these variables leads to substantial differences. To investigate this, we examine MR-LFADS(F)) trained on the pass-decision synthetic data. We perform SVD on projection matrices (Fig. S5a) from factors to rates, W, (Eq. 5), and from factors to communication, W... (Eq. 8 with r replaced by fi). Our analysis shows that W, has small singular values, indicating that certain factor dimensions contribute minimally the reconstructed rates (Fig. S5b). In contrast, Wy... exhibits relatively uniform singular values (within the same order of magnitude), suggesting that all factor dimensions are utilized in communication (Fig. S5c). This implies that some factor dimensions play a role in message inference while being largely detached from rate reconstruction (Fig. S5a).
To further investigate, we re-express the factor space of area using the left singular vectors of Wr, denoted U. In Experiment 2, MR-LFADS(F)) shown to encode both the stimulus s and decision variable d in its m P-D messages (Fig. 3e). To examine how these variables are distributed across the subspaces of U and whether this aligns with the under-constrained dimensions identified earlier, we project the factors fP onto the subspace spanned by the top singular vectors of U, denoted [[I:)),, varying k from projected values I' (Fig. S5d, cyan lines). We repeat this process in the reverse direction, projecting onto the bottom k singular vectors instead (Fig. S5d, e. orange lines). The results reveal a clear separation: decoding accuracy for improves more when projecting onto the top singular vectors, whereas decoding accuracy for d increases more rapidly when projecting onto the bottom singular vectors. This suggests that information not used for rate reconstruction such as d-is preferentially encoded in the under-constrained dimensions of the latent space, reinforeing the idea that message inference
Ning We then compute the R² values for predicting and from the
22

---

Accurate Identification of Communication Between Multiple Interacting Neural Populations
utilizes latent dimensions beyond those needed for rate prediction.
W
W Ma
Rates
Factors
Comm.
r
f :
[T]
Pm
Only top k dims are used in reconstructing rates
All dims are used in messages
fik]
Under- constrained information
pllN nenu
flN
IN H
(b)
Singular vals. for W,
(c)
Singular vals. for W
Kx
log(o,)
-2
log(o,
-0.5
-1.
88
-2.0
(d)
5 decodability from subspace projection
(e)
d decodability from subspace projection
0.6
R² [1x]] 5)
18 d))
0.8
0.4
0.6
R
R²
0.4
bottom
top
0.2
0.0
k
k
Figure SVD Analysis on MR-LFADS(F), (a) Illustraioon factors containing unconstrained information. Rates, factors, and communication are all from area i.c. r r", f= f", Mm p-- PDDD & is the effective rank of W,, (b) Ranked singular values for W.. (c) Ranked singular values for Wmm. Shaded regions indicate standard deviation across different random initializations (seeds). (d) π?² value for decoding from ff projected onto subspaces spanned by top bottom left singular vectors of W, of example seed, (c) Same as (d), but for decoding
F. Application to Multi-Region Electrophysiological Data
Data analyzed are from mouse ID:440959. For the photoinhibition experiment, we analyzed session with the following recorded regions:
Anterior Lateral Motor cortex (ALM): MOs2/3, MOs5, and M0s6 (layer
Thalamus (ALM, A): VM and VAL
Thalamus (Other, O): Anterior Ventral (AV) and Lateral Dorsal (LD)
23

---

Accurate Identification of Communication Between Multiple Interacting Neural Populations
Midbrain Reticular Nucleus (MRN): MRN
Superior Colliculus (SC): intermediate gray (SCig), intermediate white (SCiw), optic (SCop), superficial gray (SCsg), and zonal layer (SCzo)
For the model consistency experiment, we analyzed session with the these recorded regions:
ALM: secondary motor cortex, layers 2/3 and (MOs2/33, MOs5)
Thalamus: Ventral Medial (VM) and Ventral Anterior-Lateral (VAL)
MRN: MRN
Trials were filtered to include only those with durations between 4.5 and 5.5 seconds. Each trial was binned into 500 time steps, with each bin corresponding to a 10 ms interval.
F.1. Comparison of Photoinhibition Effects
We selected one session the data involving five brain regions previously implicated in decision-making task, as identified in Chen et (2024). For both control and photoinhibition trials, we only used trials from the same condition per comparison- either left-hit or right-hit-where "hit" indicates making correct choice, and "left" or "right" refers to the correct choice. For photoinhibition trials, we focused those with perturbations within the delay period, aligning all such trials to the photoinhibition onset. Firing rate f. is smoothed from raw spike counts using causal exponential filter with rate parameter of 7.1. For each region, we computed the absolute difference in trial-averaged activity between photoinhibited and control trials, averaged over neurons, to estimate the influence of photoinhibition.
F2. Consistency of Model Inference Across Random Seeds
To ensure a fair comparison between models, we evaluated each using the same hyperparameters and same number of random seeds. We then computed all pairwise similarities between inferred effectomes and messages across seeds to assess the consistency of model inference.
24


Reviews:
[{"review_id": "w81cuegGgR", "reviewer": "Reviewer_DByg", "summary": "The paper presents a modeling approach for inferring dynamics and communication between brain areas in multi-region neuroscience recordings. The approach extends the LFADS model to multiple regions via the introduction of messages between brain regions that are functions of the inferred firing rates. This has an effective low-rank communication structure, and differs from previous dynamics models that have modeled communication between brain areas in the latent space. The authors nicely point out potential limitations of that approach. The model is demonstrated in three synthetic settings, where the proposed approach performs well compared to others in terms of recovery of the underlying communication structure.", "questions": "", "limitations": "", "rating": 2, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for their thoughtful and constructive feedback. We also appreciate the reviewer’s positive assessment of our effort to address limitations in common architectural choices for communication models. Below, we respond to each of the points raised. Rebuttal figures referenced throughout our response can be accessed at the following link: https://drive.google.com/drive/folders/1R94up1vl04bkkE12tpT80YeHo39vEbqw?usp=share_link\n\n1. On the use of Jensen-Shannon divergence (JSD): Indeed, JSD is classically defined over probability distributions. We used it heuristically by treating the normalized message-norm vectors as categorical distributions and used JSD to assess how closely they match in shape. We agree that this is not the most principled application of JSD. We have now re-evaluated our approach using cosine similarity, $S_{\\cos}$, as an alternative metric for comparing inferred message-norm vectors to ground truth. We find that all model comparisons and trends remain consistent when using cosine similarity. We will update our paper, exchanging JSD for cosine similarity throughout. We have provided the corresponding comparisons between JSD and cosine similarity in Figure 1 of our rebuttal figures PDF. \n\n2. On inferred latent inputs in MR-SDS: We appreciate the reviewer's insight into inferred latent noise in MR-SDS (Karniol-Tambour el al., ICLR, 2024). Eq 4 of the MR-SDS paper states: $x^k_t = f_{kk}^{z_t}(x^k_{t-1}) + \\ldots + f_{ku}^{z_t}(u_t) + \\epsilon^k_t,\\quad \\epsilon^k_t \\sim \\mathcal{N}(0, Q_k^{z_t}).$ In the absence of provided external inputs, $u_t$, the components that could account for unobserved inputs are the dynamics, $f_{kk}^{z_t}$, or the covariance of the noise term, $Q_k^{z_t}$. Since the dynamics function is agnostic to external inputs, as the reviewer points out, $Q_k^{z_t}$ is the only parameter that can reflect the influence of unobserved inputs, potentially through a low-rank parameterization to reflect low-dimensional inputs. However, this term captures only the variance of the inferred inputs, ignoring the time-varying mean component. That is, the inputs in LFADS are parameterized as $\\hat{u}_t \\sim \\mathcal{N}(u_t | \\mu_t^u, \\sigma_t^u)$ (Eq 23 in the LFADS paper), with both  $\\mu_t^u$ and $\\sigma_t^u$ being  inferred. Therefore, it is not entirely clear to us whether the noise in MR-SDS serves an identical role to the inferred inputs in LFADS. We would greatly appreciate  further clarification from the reviewer in case we are misunderstanding important details here. \n\n3. On generative model notation (Eq. 5–7):  \nWe appreciate the reviewer’s attention to this detail and we appologize for the confusion. We do have a proper generative model, and we will update our paper to clarify as much. In our submission, Equation (6) describes the prior over the messages (and other latent variables), and Equation (7) provides the approximate posterior over the messages (which is used for inference and not if generating data from the model). The latter was not clear in our initial submission, and we will clarify this in the revised paper. In addition, we have updated Equation (5) in our paper: $\\hat{x}_t^i \\sim \\mathcal{N}(\\mu_t^{x_i}, \\mathrm{diag}(\\sigma_t^{x_i}))$, to indicate the diagonal covariance, as suggested.\n\n4. On effectome visualizations:\nWe thank the reviewer for this helpful suggestion. We agree that the thresholded diagrams serve primarily a qualitative purpose. In response, we have added unthresholded heatmap visualizations (as in Fig. 2e of our paper) to provide a more rigorous presentation of the inferred effectome (Figure 3 of our rebuttal figures PDF).", "claims_and_evidence": "Many of the claims are supported. Overall, this is important work. However, there are some issues that should be addressed before I can recommend acceptance of this paper. I have detailed these in other comments, but specifically I am referring to: issues with the generative model and questions or limitations of the evaluation techniques. \n\nSeparately, in the motivation the submission states that existing models have not incorporated unobserved inputs. However, I think this is incorrect and overinterprets the inferred inputs $\\hat{u}_t^i$. Statistically, these are latent noise inputs and other models such as MR-SDS also have inferred latent noise inputs. The primary difference here is that in this case, the latent noise inputs are assumed to be much lower dimensional than the size of the RNN state (as in the LFADS setup), which induces something akin to a low-rank noise covariance.", "methods_and_evaluation_criteria": "Overall, the methods and evaluation criteria generally make sense. The simulations are important to evaluate the model's ability to recover ground truth quantities across multiple complex settings. \n\nHowever, there is an inconsistency in the generative model that must be addressed. As written the proposed model does not appear to be consistent with a single generative model. Eq. (6) places a prior over the messages but Eq. (7) then defines a different distribution over the messages. This does not correspond to a proper generative model, as there should only be one generative distribution over the messages. In other responses, I discuss potential issues with the JS metric and the limited ability to draw conclusions from the inferred effectome plots.", "theoretical_claims": "N/A", "experimental_designs_or_analyses": "Yes, I checked the procedures for the simulated experiments and have noted potential issues in other sections.", "relation_to_broader_scientific_literature": "Neuroscience technologies are enabling neural recordings across brain areas, and it is important to develop computational methods for the analysis of such data. This paper presents an approach that is related to other work in this area. The primary innovations are extending the LFADS framework with low-dimensional inferred inputs to the multi-region setting and constraining the communication between brain areas to occur through the inferred firing rates.", "essential_references_not_discussed": "N/A", "other_strengths_and_weaknesses": "It is a nice idea to constrain the communication to be a function of the rates.", "other_comments_or_suggestions": "In equation (5), $\\Sigma_t^i$ appears to be a vector that represents the diagonal of the covariance matrix. I suggest modifying the notation here to ensure the MVN distribution has a matrix-valued covariance.  \n\nThe ground truth and inferred effectome plots with arrows between areas are helpful for qualitative presentation. However, as these are thresholded and clipped, they are not a good approach to visualize for rigorous evaluation. I'd recommend using the imshow plots as in 2e for the inferred communication.", "questions_for_authors": "Eg. (8) is not very clear. The text refers to vectors of message norms, however the JS divergence is defined in terms of distributions. If distributions are being used to compute the KL divergences, which ones are they? From my understanding, P, Q, and M are vectorized norms and it is not clear whether this is a reasonable equation.", "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": {"ethical_review_concerns": "N/A"}}, {"review_id": "pqgXZpDgR8", "reviewer": "Reviewer_KKUM", "summary": "This paper addresses the challenge of identifying connectivity patterns between brain regions from neural recordings. It introduces Multi-Region LFADS (MR-LFADS), an extension of the LFADS framework (Latent Factor Analysis via Dynamical Systems) for modeling multiple interacting neural populations. The key innovation of MR-LFADS is the introduction of constrained communication channels, ensuring that each region is driven not only by its own inferred inputs but also by message signals from other recorded regions. The model is rigorously evaluated using synthetic data across three experiments, assessing the impact of unobserved inputs, the effectiveness of communication learning, and the ability to infer random connectivity patterns. MR-LFADS demonstrates superior performance compared to a static linear model and two recent dynamic state-space models.", "questions": "", "limitations": "", "rating": 4, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for the thoughtful and insightful comments. We are encouraged that our model was recognized as a valuable contribution and that the key claims were seen as well supported. Below, we respond to each point raised. Rebuttal figures referenced throughout our response can be accessed at the following link: https://drive.google.com/drive/folders/1R94up1vl04bkkE12tpT80YeHo39vEbqw?usp=share_link\n\n1. On number of regions: We agree that evaluation on larger networks is important for assessing relevance to neuroscience applications. We will invest in improving MR-LFADS infrastructure to support larger networks, and we will provide these ongoing improvements to the community by releasing and continually updating an MR-LFADS code package. To provide immediate neuroscience relevance, see point 2 below.\n\n2. On real-world applicability: We agree that evaluation on experimental data is essential. We will update our paper to include demonstrations of MR-LFADS on multi-region electrophysiology from 5 simultaneously recorded Neuropixel probes in mice (Chen et al., Cell, 2024) performing a decision making task (Fig. 2a in our rebuttal figures PDF). We show that MR-LFADS predicts experimentally perturbed firing rates across multiple brain regions following photoinhibition of the anterior lateral motor (ALM) cortex (Fig. c-e). Critically, MR-LFADS was not trained on photoinhibition trials. Thus, MR-LFADS' ability to predict photoinhibition effects suggests identification of an accurate model of inter-region communication. Additionally, we show that MR-LFADS(R) yields more consistent solutions across random initializations (seeds) on these experimental data, as compared to the under-constrained MR-LFADS(G) (Fig. 2f,h).\n\n* Photoinhibition study:\nWe applied MR-LFADS-R to a $5$-region network spanning ALM, thalamus(A), thalamus(O), midbrain reticular nucleus (MRN), and superior colliculus (SC). Thalamus(A) includes the ventral medial (VM) and ventral anterior lateral (VAL) nuclei, which are strongly reciprocally connected with ALM (Guo et al., Nature, 2017), while Thalamus(O) includes other thalamic sub-areas, such as anterior ventral (AV) and lateral dorsal (LD) nuclei. In a subset of trials, ALM is briefly photoinhibited (Fig. 2a,c). We only fit MR-LFADS to unperturbed \"control\" trials and reserve the photoinhibition trials for post-training validation.\n\n    We first confirm that MR-LFADS fits held-out control trials well (Fig. 2b). We then simulate photoinhibition in MR-LFADS by setting the MR-LFADS ALM-region rates to be the trial-average time-varying firing rates from the photoinhibition trials from the experiment (Fig. 2d), denoted $r^\\text{inh}_*$.\n\n    These perturbed firing rates then propagate through MR-LFADS via inferred messages, ultimately predicting photoinhibited firing rates $\\hat{r}_*(\\hat{m}^\\text{inh})$. We find that the predicted ordering of ALM photoinhibition influence aligns with the experimental data (Fig. 2e).\n\n* Consistency across seeds: We fit MR-LFADS-R and MR-LFADS-G to a $3$-region network (ALM, Thalamus, MRN). Both models produced consistent effectomes across different random seeds (Fig. 2f), but MR-LFADS-R exhibited significantly greater consistency in message content (Fig. 2g,h).\n\n3. On inclusion of Granger causality (GC) and Dynamic Causal Modeling (DCM): We agree that both GC and DCM offer valuable approaches for modeling inter-region interactions. While we were very interested in including these methods, a careful and fair implementation -- particularly of the many nonlinear variants of GC and the more involved generative modeling in DCM -- was beyond the time frame of this rebuttal. Rather than risk an incomplete comparison, we will update our Related Work section to highlight these approaches, and we will pursue them in our ongoing work in this space. \n\n4. On spurious connections: We agree that it would be helpful to further analyze MR-LFADS’s behavior -- specifically, whether the model tends to miss true connections or infer spurious ones, as the reviewer suggests. We have added the following to our paper: \n\"It would also be interesting, from an interpretability perspective, to assess whether MR-LFADS tends to under- or over-estimate connectivity -- that is, whether it is more likely to miss true connections or to infer spurious ones -- and if these mismatches are driven by particular statistical properties of the signals or the dynamics of the downstream areas.”\n\n    Regarding the reviewer’s question about Exp.1: yes, the model does infer weak connections even when no ground truth connectivity is present. However, these inferred message norms are two orders of magnitude smaller than those associated with true connections. To illustrate this, we have added heatmap visualizations showing the full distribution of inferred message norms (Fig. 3).", "claims_and_evidence": "In general, the results are well supported by the experiments. The key claims made in the paper are as follows:\n1. **Inferring inputs, rather than assuming all stimuli are known to all regions, is critical for accurate communication mapping.**\nThis is primarily demonstrated by Experiment 1, which shows that when stimuli are provided to MR-LFADS (S), the model fails to recover the true communication between nodes.\n2. **Constraining messages to be derived from observable rates (MR-LFADS-R) improves identification accuracy compared to unconstrained latent communication schemes.**\nExperiment 2 supports this claim by demonstrating that the rate-based (R) version outperforms other versions. It correctly encodes the stimulus while minimizing correlation with the decision variable, ensuring more reliable communication inference.\n3. **MR-LFADS enhances the identification of true communication signals compared to prior models.**\nThis claim is supported by Experiment 3, which shows that MR-LFADS outperforms benchmark models across more than 30 randomly generated connectivity patterns. However, it is important to note that its applicability to real neuroscience studies may be limited, as the number of nodes tested (<5) is relatively small for a neuroscience application.", "methods_and_evaluation_criteria": "The benchmark and metric selections are reasonable. The innovation of inferring unknown inputs and constraining message pathways has proven effective in improving accuracy under various conditions. \n\nHowever, incorporating Granger causality as a benchmark could help broaden the paper’s appeal to a wider neuroscience audience. Additionally, it would be interesting to analyze whether the model tends to miss edges or detect spurious ones.", "theoretical_claims": "There are no formal mathematical theorems in this paper; all claims are supported empirically.", "experimental_designs_or_analyses": "The experimental design is carefully controlled, with each experiment focusing on a specific aspect of the problem. The study considers different time lags, varying connectivity weights, and different numbers of regions. However, my main concern is the real-world applicability of the method, as all testing is performed on synthetic datasets.", "relation_to_broader_scientific_literature": "The paper builds on prior models, particularly LFADS, and demonstrates superior performance compared to existing methods, making it a valuable contribution to the literature.", "essential_references_not_discussed": "The references cover the necessary background, but it would be beneficial to discuss classic connectivity methods such as Granger causality analysis and Dynamic Causal Modeling (DCM) to enhance the paper’s relevance to a broader neuroscience community.", "other_strengths_and_weaknesses": "Although all claims are well supported, as noted by the authors in the conclusion, these problems are inherently ill-posed, and the experiments rely on strong prior assumptions. As a result, some claims may not fully generalize to real-world applications.", "other_comments_or_suggestions": "NA", "questions_for_authors": "I am curious in Experiment 1, if the ground truth follows the pattern in Figure (a.ii), where there is no actual communication between nodes, can the model correctly identify this scenario, or it will still predict some spurious connections?", "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": "SIYbFN0ZON", "reviewer": "Reviewer_K91y", "summary": "In the paper \"Accurate Identification of Communication Across Multiple Interacting Neural Populations\" the authors propose an extension of LFADS to multiple regions to model region-specific inputs and cross-region interactions. The authors test their model on synthetic data with ground truth network connectivity, external inputs and cross-region interactions. They show that their model outperforms existing approaches in identifying region-specific inputs and inter-region communication.", "questions": "", "limitations": "", "rating": 5, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for their thorough and constructive feedback. We appreciate the positive assessment of the manuscript’s contribution, clarity, and experimental design. Below, we respond to each of the points raised. Rebuttal figures referenced throughout our response can be accessed at the following link: https://drive.google.com/drive/folders/1R94up1vl04bkkE12tpT80YeHo39vEbqw?usp=share_link\n\n1. On Poisson outputs:\nWe thank the reviewer for this helpful suggestion. We agree that it is important to demonstrate that the model can operate on spike train data, not just continuous activity. To this end, we applied MR-LFADS to large-scale electrophysiology recordings from Chen et al. (Cell, 2024) (Fig. 2a in our rebuttal figures PDF), which consist of spiking data and were modeled using a Poisson output distribution (Fig. 2b). We elaborate on these results below.\n\n2. On real-world applicability: We agree that evaluation on experimental data is essential. We will update our paper to include demonstrations of MR-LFADS on multi-region electrophysiology from 5 simultaneously recorded Neuropixel probes in mice (Chen et al., Cell, 2024) performing a decision making task (Fig. 2a). We show that MR-LFADS predicts experimentally perturbed firing rates across multiple brain regions following photoinhibition of the anterior lateral motor (ALM) cortex (Fig. c-e). Critically, MR-LFADS was not trained on photoinhibition trials. Thus, MR-LFADS' ability to predict photoinhibition effects suggests identification of an accurate model of inter-region communication. Additionally, we show that MR-LFADS(R) yields more consistent solutions across random initializations (seeds) on these experimental data, as compared to the under-constrained MR-LFADS(G) (Fig. 2f,h).\n\n* Photoinhibition study:\nWe applied MR-LFADS-R to a $5$-region network spanning ALM, thalamus(A), thalamus(O), midbrain reticular nucleus (MRN), and superior colliculus (SC). Thalamus(A) includes the ventral medial (VM) and ventral anterior lateral (VAL) nuclei, which are strongly reciprocally connected with ALM (Guo et al., Nature, 2017), while Thalamus(O) includes other thalamic sub-areas, such as anterior ventral (AV) and lateral dorsal (LD) nuclei. In a subset of trials, ALM is briefly photoinhibited (Fig. 2a,c). We only fit MR-LFADS to unperturbed \"control\" trials and reserve the photoinhibition trials for post-training validation.\n\n    We first confirm that MR-LFADS fits held-out control trials well (Fig. 2b). We then simulate photoinhibition in MR-LFADS by setting the MR-LFADS ALM-region rates to be the trial-average time-varying firing rates from the photoinhibition trials from the experiment (Fig. 2d), denoted $r^\\text{inh}_*$.\n\n    These perturbed firing rates then propagate through MR-LFADS via inferred messages, ultimately predicting photoinhibited firing rates $\\hat{r}_*(\\hat{m}^\\text{inh})$. We find that the predicted ordering of ALM photoinhibition influence aligns with the experimental data (Fig. 2e).\n\n* Consistency across seeds: We fit MR-LFADS-R and MR-LFADS-G to a $3$-region network (ALM, Thalamus, MRN). Both models produced consistent effectomes across different random seeds (Fig. 2f), but MR-LFADS-R exhibited significantly greater consistency in message content (Fig. 2g,h).\n\n3. On notational issues and typos:\nWe greatly appreciate the reviewer’s attention to detail. We have corrected the identified issues and revised the notation for the priors over inferred inputs and initial conditions in the updated manuscript (Fig. 4). Typographical errors on lines 28 and 633 have also been fixed. Additionally, the reviewer was correct to point out the inconsistency in the KL divergence description on line 211 -- we have updated this to reflect the correct direction of the penalty: \"Therefore, we carefully set the KL penalty weight for inferred inputs ($\\lambda_\\text{KL}^{\\hat{u}}$) to be higher than that for messages ($\\lambda_\\text{KL}^{\\hat{m}}$).\"\n\n4. On clarifications for incorrect inferences due to incomplete recordings: We appreciate this thoughtful question from the reviewer. The interpretation offered — “I suppose you mean that if you were to add more complete recordings, the inferred connectivity would change” — captures our intended meaning precisely. In practice, we often cannot guarantee that all relevant areas are recorded. As a result, inputs from unobserved regions must be accounted for (as latent variable inferred inputs) to minimize the deleterious effects of unobserved nuisance variables. We will update our paper to clarify this point.", "claims_and_evidence": "The authors claim that their model infers region-specific inputs from unobserved areas, outperforming approaches that rely on manual input specification. They also claim that their model constrains inferred inputs and inter-region communication more accurately than existing approaches. The authors evaluate their model on synthetic data only and show that it outperforms previous communication models across various neuroscience related synthetic tasks. While they explore a wide range of synthetic tasks, there is no real data analysis to support their claims. Moreover, while they propose a model variant for Poisson output, they only consider Gaussian output in their experiments.", "methods_and_evaluation_criteria": "The proposed MR-LFADS model extends the LFADS model by replicating a version of the LFADS model for each region and adding constrained communication channels between the regions. They propose variants of the single-region LFADS models for Poisson outputs and for Gaussian outputs. Communication messages are multivariate Gaussians variables with mean and covariance depending on the time varying average activities of other regions. The inference scheme follows the typical VAE variational inference. The evaluation metrics include the similarity between inferred pairwise message norms and ground truth quantified using Jenson-Shannon divergence, and the correspondence between inferred messages and ground truth quantified using linear regression to predict ground truth messages from inferred messages. These are reasonable evaluation metrics for the proposed model.", "theoretical_claims": "There are no theoretical claims in the paper.", "experimental_designs_or_analyses": "The authors evaluate their model on three experiments with synthetic data with known ground truth network connectivity, external inputs and cross-region interactions. They compare their model against reduced-rank regression, multi-population sticky recurrent SLDS and multi-region switching dynamical systems. The first experiments involves a synthetic memory network trained to recall private stimuli and received messages to assess the impact of unobserved inputs. The second experiment involves a synthetic pass decision network to assess the effect of direct constraint of messages via observed data. The third experiment involves synthetic randomized networks trained on randomly selected tasks to assess flexibility in coordination. The experiments are designed to highlight common failure modes of identifying region-specific inputs and inter-region communication. Overall, the authors succeed in showing that their model outperforms existing approaches in the synthetic tasks they consider. However, for the third experiment, the results are not as clear cut as in the previous experiments. While they do not present a systematic ablation study, they consider different variants of the model in the experiments to assess the impact of e.g. the KL penalty weight for messages.", "relation_to_broader_scientific_literature": "The LFADS model is a key ingredient in the proposed MR-LFADS model. The authors go beyond the LFADS model by considering multiple regions and communication between them. This generalizes the original framework and makes it useful for learning communication across multiple interacting neural populations. The authors show that this approach has advantages over existing state-of-the-art methods for communication analysis in neuroscience.", "essential_references_not_discussed": "To the best of my knowledge, the authors have discussed all essential references.", "other_strengths_and_weaknesses": "The paper is very well written and easy to follow. The authors provide a detailed description of the model and the inference scheme. The experiments are well designed and the results are clearly presented. In my view, the main weaknesses are two-fold: the lack of real data analysis and in the evaluation the focus on Gaussian output only. I don't think these are fatal flaws, but the provided analysis did not fully overcome concerns about the validity of the claims in more realistic settings.", "other_comments_or_suggestions": "Line 28, right column: \"should should\" -> \"should\"\nFig 1: g_0 and u hat are multivariate but the sigmas in the Gaussians suggest univariate variables. The priors on the other hand use the identity matrix in the covariances. This looks inconsistent.\nLine 633, left column: \"doen\" -> \"done\"", "questions_for_authors": "1. Line 81, right column: You claim that existing methods risk incorrect inferences about connectivity and communication due to incomplete recordings. Could you elaborate more on the functional aspects of such incorrect inferences due to incomplete recordings? For models providing purely functional descriptions of neural activity, there is usually no claim about the physiological connectivity of the regions, nor about causality. What does it even mean for an inference to be incorrect in the context of the functional model? I suppose you mean that if you were to add more complete recordings, the inferred connectivity would change. Could you clarify this point?\n2. Line 211, left column: You said that you set the KL penalty weight for messages higher than for inferred inputs to encourage the model to use messages. But shouldn't it be the other way around, i.e. the KL penalty weight for messages should be lower? This is also suggested by the more detailed description in the appendix.", "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": "DaE1DGeEtY", "reviewer": "Reviewer_yYC3", "summary": "The paper introduces Multi-Region Latent Factor Analysis via Dynamical Systems (MR-LFADS), a novel model for accurately identifying inter-region neural communication in multi-region recordings. MR-LFADS extends the existing LFADS framework by jointly inferring region-specific inputs from unrecorded brain areas and explicitly constraining inter-region communication signals based on observed neural activity, effectively mitigating biases from incomplete recordings and under-constrained latent representations. The key conceptual contribution is the explicit inference of region-specific inputs and rate-constrained inter-region communication messages, reducing ambiguity and misidentification inherent in previous methods. By better constraining inferred communication signals through observable neural activity, MR-LFADS provides a robust tool for analyzing multi-region neural recordings and elucidating principles of brain-wide information processing.", "questions": "", "limitations": "", "rating": 3, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for their thoughtful and constructive feedback. We’re encouraged by the reviewer’s comments highlighting MR-LFADS as a novel extension of LFADS, and by the positive assessments on our synthetic dataset. Below, we respond to the reviewer's concerns on the lack of comparison to experimental data. Rebuttal figures referenced throughout our response can be accessed at the following link: https://drive.google.com/drive/folders/1R94up1vl04bkkE12tpT80YeHo39vEbqw?usp=share_link\n\nWe agree that evaluation on experimental data is essential. We will update our paper to include demonstrations of MR-LFADS on multi-region electrophysiology from 5 simultaneously recorded Neuropixel probes in mice (Chen et al., Cell, 2024) performing a decision making task (Fig. 2a in our rebuttal figures PDF). We show that MR-LFADS predicts experimentally perturbed firing rates across multiple brain regions following photoinhibition of the anterior lateral motor (ALM) cortex (Fig. c-e). Critically, MR-LFADS was not trained on photoinhibition trials. Thus, MR-LFADS' ability to predict photoinhibition effects suggests identification of an accurate model of inter-region communication. Additionally, we show that MR-LFADS(R) yields more consistent solutions across random initializations (seeds) on these experimental data, as compared to the under-constrained MR-LFADS(G) (Fig. 2f,h).\n\n* Photoinhibition study:\nWe applied MR-LFADS-R to a $5$-region network spanning ALM, thalamus(A), thalamus(O), midbrain reticular nucleus (MRN), and superior colliculus (SC). Thalamus(A) includes the ventral medial (VM) and ventral anterior lateral (VAL) nuclei, which are strongly reciprocally connected with ALM (Guo et al., Nature, 2017), while Thalamus(O) includes other thalamic sub-areas, such as anterior ventral (AV) and lateral dorsal (LD) nuclei. In a subset of trials, ALM is briefly photoinhibited (Fig. 2a,c). We only fit MR-LFADS to unperturbed \"control\" trials and reserve the photoinhibition trials for post-training validation.\n\n    We first confirm that MR-LFADS fits held-out control trials well (Fig. 2b). We then simulate photoinhibition in MR-LFADS by setting the MR-LFADS ALM-region rates to be the trial-average time-varying firing rates from the photoinhibition trials from the experiment (Fig. 2d), denoted $r^\\text{inh}_*$.\n\n    These perturbed firing rates then propagate through MR-LFADS via inferred messages, ultimately predicting photoinhibited firing rates $\\hat{r}_*(\\hat{m}^\\text{inh})$. We find that the predicted ordering of ALM photoinhibition influence aligns with the experimental data (Fig. 2e).\n\n* Consistency across seeds: We fit MR-LFADS-R and MR-LFADS-G to a $3$-region network (ALM, Thalamus, MRN). Both models produced consistent effectomes across different random seeds (Fig. 2f), but MR-LFADS-R exhibited significantly greater consistency in message content (Fig. 2g,h).", "claims_and_evidence": "Yes, the claims made in the paper are generally supported by clear and convincing evidence.", "methods_and_evaluation_criteria": "Yes, the methods and evaluation criteria chosen by the authors are sensible and well-aligned with their stated goals. The synthetic datasets are thoughtfully designed to capture specific challenges faced by communication models, such as the presence of unobserved inputs, nonlinear dynamics, and temporal structure in neural interactions.", "theoretical_claims": "Yes. This paper does not contain explicit theoretical proofs or formal theoretical claims requiring validation.", "experimental_designs_or_analyses": "While the authors present extensive validation using synthetic datasets, which clearly demonstrate the model's capabilities in controlled scenarios, but the evaluation is solely based on synthetic data. Although synthetic datasets have the advantage of known ground truths, allowing precise assessment of the model’s accuracy, they inherently lack the full complexity and variability present in actual neural recordings. Evaluating MR-LFADS on neural recordings may scientifically get meaningful insights into brain-wide neural communication that are currently unexplored or poorly understood.", "relation_to_broader_scientific_literature": "The paper introduces MR-LFADS, extending LFADS to infer communication across multiple neural regions by explicitly modeling unobserved inputs and constraining inter-region messages through observed neural activity. It situates itself relative to existing methods and demonstrates improved performance on synthetic benchmarks that emulate cognitive neuroscience tasks", "essential_references_not_discussed": "No", "other_strengths_and_weaknesses": "Strengths:\n- Creatively adapts the established LFADS framework to explicitly model communication among multiple interacting neural regions.\n\nWeakness:\n- The lack of validation on actual neural recordings somewhat limits the impact. Demonstrating applicability to neural data could substantially strengthen claims regarding scientific insight.", "other_comments_or_suggestions": "N/A", "questions_for_authors": "N/A", "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:

Meta-review: This study addresses the challenge of identifying connectivity patterns between brain regions from neural recordings. It introduces Multi-Region LFADS, an extension of the Latent Factor Analysis via Dynamical Systems for modeling multiple interacting neural populations. The key innovation is the introduction of constrained communication channels, ensuring that each region can be activated not only by its own inferred inputs but also by message signals from other recorded regions. The model is evaluated using synthetic data assessing the impact of unobserved inputs, the effectiveness of communication learning, and the ability to infer random connectivity patterns. The author demonstrate that their new method performs better compared to a static linear model and two recent dynamic state-space models. The reviewers appreciated that the paper makes a case for being accepted in ICML but also raised suggestions both of citation of literature and of additional analyses and results to be added, which are expected to be implement in a revision.
Decision: Accept (poster)