Learning Time-Varying Multi-Region Brain Communications  via Scalable Markovian Gaussian Processes
Authors: Weihan Li
Yule Wang
Chengrui Li
Anqi Wu
Venue: ICML 2025 oral (2025)
Keywords: Multiple Brain Region Communications; Markovian Gaussian Processes; State Space Model
Primary area: applications->neuroscience_cognitive_science
TL;DR: We developed a scalable method using Markovian Gaussian Processes (State Space Model) to track how multiple brain regions communicate dynamically over time.

Abstract: Understanding and constructing brain communications that capture dynamic communications across multiple regions is fundamental to modern system neuroscience, yet current methods struggle to find time-varying region-level communications or scale to large neural datasets with long recording durations. We present a novel framework using Markovian Gaussian Processes to learn brain communications with time-varying temporal delays from multi-region neural recordings, named Adaptive Delay Model (ADM). Our method combines Gaussian Processes with State Space Models and employs parallel scan inference algorithms, enabling efficient scaling to large datasets while identifying concurrent  communication patterns that evolve over time. This time-varying approach captures how brain region interactions shift dynamically during cognitive processes. Validated on synthetic and multi-region neural recordings datasets, our approach discovers both the directionality and temporal dynamics of neural communication. This work advances our understanding of distributed neural computation and provides a scalable tool for analyzing dynamic brain networks.

Learning Time- Varying Multi-Region Brain Communications via Scalable Markovian Gaussian Processes
Weihan Li!
Yule Wang'
Chengrui Li! Angi Wu!
Abstract
enable simultaneous measurement of neural activity across numerous brain areas with unprecedented temporal and spa- tial resolution (Steinmetz et al., 2021; Siegle et al., 2021; Li et al., 2024a; Nejatbakhsh et al., 2024), These large-scale neural recordings necessitate computational methods capa- ble of uncovering and characterizing the dynamic nature of neural communication within the brain.
Understanding and constructing brain commu- nications that capture dynamic communications across multiple regions is fundamental to mod- ern system neuroscience, yet current methods struggle to find time-varying region-level com- munications or scale to large neural datasets with long recording durations. We present a novel framework using Markovian Gaussian Pro- cesses to learn brain communications with time- varying temporal delays from multi-region neu- ral recordings, named Adaptive Delay Model (ADM). Our method combines Gaussian Pro- cesses with State Space Models and employs parallel sean inference algorithms, enabling ef- ficient scaling to large datasets while identifying concurrent communication patterns that evolve over time. This time-varying approach captures how brain region interactions shift dynamically during cognitive processes. Validated on syn- thetic and multi-region neural recordings datasets, our approach discovers both the directionality and temporal dynamics of neural communication. This work advances our understanding of dis- tributed neural computation and provides scal- able tool for analyzing dynamic brain networks. Code is available at https: : //gi thub. com/ Delay -Model.
Latent representations offer a promising approach build- ing time-varying multi-region communications (Wang et al., 2023; 2024; Zhang et al., 2024). In the brain, communi- cation patterns between brain regions manifest at varying temporal scales: some regions exhibit fast, synchronous interactions with short delays, indicative of strong func- tional coupling, while others display slower interactions with longer delays, reflecting more indirect relationships. Such communications provide framework for understand- ing how these diverse communication patterns evolve over time, capturing both feedforward and feedback pathways that may shift during different cognitive states (Lillicrap et al., 2020; Liu et al., 2025).
Current computational approaches for modeling multi- region communications can be broadly categorized into non- delay models and delay models. Non-delay models, such as mp-srSLDS (Glaser et al., 2020) and MR-SDS (Karniol- Tambour et al.), do not explicitly incorporate temporal de- lays when capturing latent dynamics across regions. As a result, they can only learn the directional information of communications. Delay models, including DLAG (Gokcen et al., 2022), m-DLAG (Gokcen et al., 2024a), their ap- proximated versions (Gokeen et al., 2024b), and MRM-GP (Li et al., 2024b), introduce mechanisms to learn temporal delays between pairs of communications. Learning delays provides not only directional information but also insights into relative communication speeds, which can help assess the strength of functional coupling between brain regions.
BRAINML- MRAINML-GT/Adappe GT/AdaatiAeeeeeeeeeea sesettiiieeelllyyy
1. Introduction
Modern system neuroscience faces a significant challenge in discovering communication patterns that capture dynamic interactions across multiple regions. Understanding these time-varying communications has become increasingly crit- ical with the advent of advanced recording technologies that
Although these methods have shown potential in capturing inter-region brain communications, delay models, such as DLAG and its variants, do not account for the dynamic na- ture of multi-region communication. They fail to model time-varying temporal dependencies and communication patterns, which are crucial for understanding neural process- ing. Non-delay models, such as MR-SDS and mp- mp-srSLDS.
School of Computational Science & Engineering, Georgia Institute of Technology, Atlanta, USA. Correspondence to: Anqi Wu <anqiwuu <angiwu@@aatcc..duss gatech.edu>.
Proceedings of the nrnd International Conference Machine Learning. Vancouver, Canada. PMLR 267, 2025. Copyright 2025 by the author(s).

---

Learning Time-Varying Multi-Region Brain Communications via Scalable Markovian Gaussian Processes
introduce dynamic message flow among regions but don't model delays between communications and assume sin- gle communication subspace, meaning they cannot capture concurrent communications over different subspaces.
icantly accelerate computation and reducing complexity to Ollog T)- This approach enables efficient analysis of long-duration recordings while capturing dynamic commu- nication patterns.
MRM-GP is a specific case that integrates a Gaussian Pro- cess (GP) with a State Space Model (SSM), making it a delay model with discrete changing of phase delays. MRM- GP has two key limitations: (1) it only supports kernels that separable across spatial and temporal domains, restrict- ing its applicability in the frequency domain by learning phase delays, whereas temporal delays are more generaliz- able for neuroscience applications: (2) it assumes that the delays across all communication subspaces across regions change synchronously with hidden state transition. How- ever, this assumption does not reflect the true brain mech- anism, where different communication pathways between regions can operate asynchronously.
We validate our model using neural recordings from multiple regions of the brain during visual processing tasks (Semedo et al., 2019; Siegle et al., 2021). Our results demonstrate the method's capability to uncover how information flow patterns dynamically change across multi-region networks, offering new insights into the temporal organization of large- scale neural circuits and advancing our understanding of distributed neural computation.
In summary, the key contributions of our work are:
We establish a universal connection between arbitrary temporally stationary GPs and SSMs, which has broader im- plications for other domains where computational efficiency is priority.
In terms of computational efficiency, DLAG and mDLAG, which are GP-based methods, incur an ) computa- tional cost, where T is the number of time samples, mak- ing them challenging to scale. Approximated GP and SSM-based methods, such as MRM-GP, MR-SDS, and mp- srSLDS, reduce the cost to OT)). However, they rely sequential inference, which remains inefficient for large neu- ral datasets with long recording durations. Moreover, the reliance on discrete hidden states, e.g., MRM-GP, introduces inefficiencies during inference.
Our model discovers time-varying multi- commu- nications from latent representations of neural recordings without introducing additional discrete hidden states in the SSM.
We propose a scalable method for analyzing multi-region communications in large-scale neural data with O T) complexity.
2. Method
To address these limitations, we propose the Adaptive De- lay Model (ADM). It falls within the family of Markovian Gaussian process-based delay models (Li et al., 2024b), but incorporates a continuous time-varying temporal delay. Thus unlike static communication models commonly used in neuroscience, ADM can accommodate communication patterns with varying temporal characteristics, offering a more flexible and biologically relevant framework.
We hegin by demonstrating how a Gaussian Process (GP) with Factor Analysis (FA) model can be used to capture static brain communication across regions (Section 2.1). Next, we establish universal connection between GP and a State Space Model (SSM), referred to as the Marko- vian Gaussian Process (Section 2.2). Finally, we illustrate how time-varying brain communication can be modeled (Section 2.3).
Markovian Gaussian Process (Markovian GP) models inte- grate the expressive power of GPs with the computational efficiency of SSMs, facilitating scalable analysis of large- scale neural recordings while discovering multiple evoly- ing communication patterns. However, existing Markovian GP approaches either rely on single-output kernels or re- quire multi-output kernels be separable across spatial and temporal domains (Solin et al., 2016; Loper al., 2021; Dowling et al., 2021; 2023; Li et al., 2024b). In this pa- per, we introduce novel, universal connection between arbitrary temporally stationary GPs and SSMs, making the framework highly flexible and broadly applicable to various neuroscience problems.
2.1. Gaussian Process Factor Analysis for Brain Communications
The Gaussian Process Factor Analysis (GPFA) for modeling brain communication employs the multi-output Squared Exponential (MOSE) kernel (Gokcen et al., 2022);
(T + 0,g) (T;))
(T)(T) exp
(I)
2(2
where i, represent two brain regions, 7 = t- t' is the time difference, O:ş is the temporal delay between regions and j. and I is the length scale shared across all regions. This MOSE kernel allows the identification of temporal delays that characterize communication dynamics across regions.
Additionally, we apply an advanced inference method for Markovian GP, leveraging parallel scan algorithms (Blel- 2020) to signif-
loch, 1990; Särkkä &
Our goal is to learn M N latent variables, æ E RMN
2

---

Learning Time-Varying Multi-Region Brain Communications via Scalable Markovian Gaussian Processes
from neural recordings y RDxT across N regions. Each region is associated with M latent variables.
modeled as a GP with MOSE kernel:
K(0)
K(-1) ((0)
((-T+1)
A typical assumption for z is to decompose it into two com- ponents: across-region variables and within-region variables (Gokcen et al., 2022; Li et al., 2024b). The across-region variables, x" E I Nx,", capture shared neural activity that reflects communication between brain regions. These variables exhibit the similar dynamies across regions, dif- fering only in temporal delays. The within-region variables, x E R "et xT represent neural activity unique to individ- ual regions and are independent of other regions. Together, these components form the latent representation of neural recordings, x = [>=", ",""" "", where ma + my = M. The re- lationship between y and x is then modeled using Factor Analysis (FA) model:
K(1)
K ((T++2)
GP(O,
) (4)
K(T - 1) K(T - 2)
K(0)
where each K(T) E RNXN is MOSE kernel in Eq. with interval over N brain regions. Our goal is to find state- space approximation of a which follows a Multi-Order SSM structure:
Q m,
Ap m,t-p + qe q: ~ ~2,, Q), (5)
where P represents the number of orders, A Ap R-NxN are the transition matrices, and Q E RRNN is the process noise matrix.
y y-Ca Ca Cd+
(2)
where C E DDMN is a block-diagonal matrix C = diag ...., CN with each C' representing the mapping from region i's latent variable to its neural recordings. Addi- tionally, d TR. Dx1 is a bias term, and NO, V) repre- sents Gaussian noise with diagonal covariance V e DXD
Determining an State Space Model using Kernels. To estimate transition matrices and measurement using a we can consider the SSM in Eq. 5 as a regression model (Neumaier & Schneider, 2001):
The across-region variables x" are designed to capture com- munication patterns among regions. For the m-th group of across-region variables, xn E gNx,, the activity from each region exhibits spatial correlations with the N V-1o other regions. These variables are modeled Gaussian Process (GP) with the MOSE kernel, which captures the temporal delay characteristics of across-region communication.
m,t Gu, + 4o q. ~o0 Q).
(6)
where G E NxNP is the regression coefficient and de E RNPX1 is the predictor:
G = [Ap A A A1), T.T
On the other hand, the m-th group of within-region vari- ables ah E RN representing region-specton activity, is modeled independently across regions using a single-output Squared Exponential (SE) kernel:
U, c m
x T.T m.
P41
Our ultimate goal is to use K(T) to represent G and Q First, we can represent G and Q as functions of x Con- cretely, given T samples, Im.1..... dictor matrix as V [UP+1, UT!
K (rn = exp
(3)
Dm we define pre-
2/2
UT] E R P) and target observation matrix as W [xt Pmm ................................ m,rl (T-P))
Additionally, independence is assumed across different groups of both the across-region and within-region variables, with each group m governed by distinct kernel parameters.
Then, can represent the regression model in Eq. 6 in the matrix form: W GV + R, where R is the residual matrix. By doing so, we can estimate coefficient matrix G and the process noise matrix Q by least squares estimation:
By explicitly separating the across-region and within-region latent variables, this framework offers a clear representation of across-region communication and region-specific dynam- ics, enabling more interpretable analysis of multi- multi-region neural recordings.
RR
G C=WVT(VVT)
Q Q= T
(8)
2.2. Connect Gaussian Process with State Space Model
where R W GV denotes an estimate of the residual matrix, and its covariance is an estimate process noise matrix Q. Now, if we can represent WVT and VVT with K(7), we will achieve the ultimate goal.
We develop novel universal connection between arbitrary temporally stationary GPs and SSMs, enabling us to effi- ciently model both across- and within-region dynamics.
Gaussian Process and State-Space Approximation. The m-th group of across-region variables, 2 E NxT are
r in V and W is modeled sample in the GP in Eq. 4. We can represent VVT NNPXNP
Since each sample ait,,
3

---

Learning Time-Varying Multi-Region Brain Communications via Scalable Markovian Gaussian Processes
and WVT E RRNNNP using (7) as (full derivations see
2.3. Time- Mirr-ng Brain Communications
Appendix A):
Since the SSM in Eq. 10 follows a discrete structure, we can extend it to incorporate time-varying transition and process noise matrices follows: x Fm,t
((O) K(1)
K(-1)
( PK-P+ + K (
[(-P+1] 1)
((0)
VV &
A. Too
q. No, Qt),
(12)
K(P-1) 1) K(P- 2)
K(0)
I m,
H
WV & K(P) -(P-1). 1)... K(I)]
This formulation introduces a time-varying MOSE kernel, where the temporal delay parameter Ozg. evolves over time. In other words, at each time step f, we construct a Marko- vian GP (or SSM) as described Eq. 10, conditioned on the MOSE kernel specific that time t. Additionally, the length scale parameter / is held constant over time, which limits the flexibility of each Å, By sharing ( across time points, we constrain the temporal evolution of the delay parame- ters, promoting smoother dynamics and reducing the risk of misatrributing variability in the messages to fluctuations in delays.
Notably, each K(T) E RNNN, where E1 depends only on the number of brain regions N and can be efficiently computed using the stationary kernel func- tion employed in the GP. Furthermore, K(T) can represent any stationary temporal kernel, establishing a universal con- nection between GPs and SSMs. We apply this universal conversion to various kernels in GP regression task, see Appendix D for details.
This approach enables the model to learn tar-vngying tem- poral delays, effectively capturing the dynamics of multi- region brain communications. Compared to modeling time- varying phase delays using hidden discrete states (Li et al., 2024b), our method is more flexible, as it does not assume that each group of across ceross-region communications, Im un- dergoes simultaneous delay changes during state transitions.
Markovian Across-region Communcations. Now, the transition matrices and the measurement matrix in Eq. are uniquely determined by the kerne functions of GP by Eq. 8 and Eq. 9. Moreover, we can rewrite the SSM in Eq. 5 to an SSM with a Markovian structure, resulting Markovian Gaussian Process (Markovian GP) (Zhao, 2021):
E Fm.t
A 1-1-1 1-1 + Tc: le 2 (O Q),
Importantly, this computation can be efficiently performed in vectorized form across all T time steps, ensuring mini- mal impact on overall computational efficiency. In the FA model, the projection matrix C E D the bias term d E Dx11, and the observation Gaussian noise e E IR Dx1 remain innariitaii
(10)
am, H m,
where H E NNNP denotes mask matrix, A E PNPXNP is structured as a controllable canonical form (Grewal & Andrews, 2014) and small constants are added to Q RN NNXNP so it matches the shape of A:
3. Inference
A
A2
A A
A P
IN
0
0
Now, having established the connection between the m-th group of across-region communications and within-region neural activity to the Markovian GP (or SSM), as described in Eq. 12, the next step is to efficiently learn the latent variables and model parameters.
...
O
IN
0
0
0
0
0
(11)
0
0
IN
0
Q = 0
Q
xIN(P-1)
0
H = [IN 0].
Our model, ADM, formulated as an SSM, offers a signifi- cant advantage: it can learn the parameters using either sequential estimation method with complexity O(T) or parallel computation method with complexity O(log 7). On modern hardware, the parallel approach is consistently faster due to its efficient utilization of computational resources.
with E NNN denoting the identity matrix and o a small constant added for numerical stability. Notably, although we use Markovian structure to represent stationary GP, our method still incorporates information from multiple orders.
Markovian Within-region Neural Activity. Similarly, the state-space approximation of the m-th group of within- region variables, x E IR NxT can be seen as a specific case of across-region variables. In this case, each dimension, "... E T is independently modeled as a Markovian GP with a scalar single-output SE kernel (Fq.3) 3).
Parameter Settings The model parameters, denoted as O, include the kernel parameters Of and 1² from each la- tent dimension k. which define the transition matrix A, and the process noise covariance matrix Q, for each across- or within-region latent group. Additionally, the Factor Anal- ysis (FA) parameters include the projection matrix C, the

---

Learning Time-Varying Multi-Region Brain Communications via Scalable Markovian Gaussian Processes
bias term d, and the diagonal matrix V. The model also has hyperparameter P. representing the order of the autore- gressive process. To better understand the effect of different P values, we generate samples from our model for various P values, as shown in Appendix C.
by both fast and slow features.
Two Brain Regions (Semedo et al., 2019; Zandvakili & Kohn, 2019): Simultaneous spike train recordings from a monkey's 's primary visual area (VI) and secondary visual cortex (V2), with a 6Hz drifting grating as the external stimulus.
Parallel Scan Kalman EM Algorithm Given neural recordings y E TRDxT, our goal is to estimate the latent brain communications a RM along with the model parameters + To achieve this, we use the parallel scan- based Kalman Sopectation-Maxximize algorithm (Särkkä & García-Fernánder 2020), which in- troduces a parallel scan version of the Kalman Filter and Smoother: Specifically, in the E-step, we apply the parallel Kalman Filter and Smoother to infer the latent variables and expected log-likelinge likelihood. In the M-step, we update the kernel parameters using gradient descent and optimize the Factor Analysis (FA) parameters through closed-fome linear regression. See Appendix B for details.
Five Brain Regions (Siegle al.. 2021): Simultane- ous spike train recordings from a mouse's primary visual cortex (VISp), rostrolateral area (VISrl), anterolateral area (VISal), posteromedial area (VISpm), and anteromedial area (VISam), with a 4Hz drifting grating as the external stimu- lus.
(EM) inference
Baselines.
We compare our model with three methods:
DLAG (Gokcen et al., 2022): A Gaussian Process Factor Analysis model with a MOSE kernel, designed for neu- ral recordings from two brain regions. It can be used to learn both across-region and within-region latent communi- cations.
The objective of the Kalman Filter is to compute the pos- terior density p(as (o:)), given the neural data up to time step I, while the Kalman Smoother computes the posterior density p( (Borr) for all time steps. Traditionally, both filtering and smoothing are computed in (T) time using se- quential updates. However, sequential computation is often inefficient compared to parallel computation, particularly on modern hardware architectures (Chen et al., 2024).
mDLAG (Gokcen et al., 2024a): An extension of DLAG that supports more than two brain regions with a different inference approach. Unlike DLAG, it assumes all latent variables correspond across-region communications and does not explicitly model within-region dynamics.
MRM-GP (Li et al., 2024b): An approximation of a Gaussian Process with a Cross-Spectral Mixture (CSM) kernel (Ulrich et al., 2015), formulated as an SSM with O(T) complexity. designed to learn frequency- communications between two brain regions and can capture both across-region and within-region latent dynamics.
To address this inefficiency, Särkkä & Garcfa-Fernándes (2020) demonstrates that the sequential updates of the Kalman Filter and Smoother can be reformulated as an associative operator, enabling the use of the parallel sean algorithm (Blelloch, 1990). Consequently, the time com- plexity and memory cost of our model are given by:
Evaluation. We evaluate our model and baseline mod- els by randomly splitting the data into training, validation, and testing sets with ratio of 0.8, 0.1, and 0.1, respec- tively. Since all models assume a linear/Gaussian relation- ship between the latent variables and observed data, we assess their performance by computing the observation test ligeliiel:oood: L.L Brea), where Fes represents the inferred test latent variables, and Viest denotes the test neural recordings. To mitigate randomness, we report the average test log- likeliioddo over five different random seeds.
Time complexity: O Log
(13)
Memory complexity: O(MN2P21),
(14)
where N is the number of regions, M is the group number of across- and within-region latent dynamics, and P is the SSM order in Eq. 5. The cubic cost arises from Eq. 8. How- ever, we will show in the experimental section, the order parameter P and N can be significantly smaller compared to T' while still achieving strong generative and inference performance. Thus, the cubic cost does not pose a major computational bottleneck.
4.1. Synthetic Data
In this section, we simulate common phenomenon in neu- rosceecce where brain region communications are dynamic (Parra & Tobar, 2017). Our goal is to evaluate our model's s ability to recover time-varying temporal delays and latent variables.
4. Experiment
Datasets.
We evaluate our model on three datasets.
Synthetic Data: We generate synthetic data that incorpo- rate both across-region communications and within-region neural activities, along with time- varying temporal delays, to simulate dynamic brain communications characterized
Experimental setup. We generate 120 independent trials for two brain regions (N = 2) with an order of P = 5
5

---

Learning Time-Van Varring Multi-Region Brain Communications via Scalable Markovian Gaussian Processes
(A) time oin
Time bin
Our model effectively captures time-varying communica- tions for both latent dynamics and delay. See Appendix E for within-region neural activities.
scross-regionn temporal communications delays
esfimated deley
delay 4.
Figure ((B) shows the test log- likelihood summed over trials and time bins. The results indicate that performance remains relatively stable for P E [2,5], except for P = 1, which yields the lowest performance.
region 1
true latent
rogion 2
true latent
10 Bins
Figure I(C) compares the time costs of the sequential and parallel scan-based Kalman EM algorithms with GPU par- allelization. We generate synthetic data with up to T = 600 time bins. The results demonstrate that the parallel version is significantly faster than the sequential update.
(C)
test test logglllkkklhodd
time in minules
16
sequential
paratel
400k-
4
350k
Finally, Appendix G presents a more complex case with five-region synthetic dataset and evaluates the model's per- formance as function of the number of regions, latent variables, and data length.
128
256
512
Figure/ 1. Evaluation of the ADM on synthetic brain communica- tion data. (A) Estimated across-region communications, temporal delays, and ground truth over T bins for P = 5. The shaded area represents the variance of the estimated delay across five dif- ferent runs. (B) Test log- likelihood summed over trials and time bins, showing that performance remains stable for different P. (C) Comparison of the runtime for sequential and parallel sean- based Kalman EM across different values of T, showing that the parallel version is significantly faster, especially as T increases.
4.2. Two Brain Regions
In this section, we investigate the interactions between the mouse's 's primary visual area (V1) and secondary visual cor- tex (V2) in response to 6Hz drifting grating. Additionally, we compare our model's performance and inference time with MRM-GP and DLAG.
Experimental setup. We use smoothed multi-region spike trains from session (06r001p26 with an orientation of 0 This dataset consists of 400 trials, each containing 64 time bins (20 ms per bin), with 72 VI neurons and 22 V2 neu- rons. The monkey begins receiving the visual stimulus (drift- ing gratings) at the first time bin, and the stimulus persists throughout all 64 time bins. The number of across-region and within-region latent dynamics follows previous works (Gokcen et al., 2022; Li et al., 2024b). where m, = 2 and my = 2. The order P = 4 is selected based on performance evaluation on the validation dataset.
and T = 200 time bins. Each region contains 50 neurons, with m. = 2 groups of across-region communications and my = group of wethin-region variables. For across-region communications, the first group represents forward commu- nication, characterized by a larger positive delay of 5 bins and a smaller positive delay of 1 bin during time bins 30 to 70. The second group represents feedback communica- tion, with a larger negative delay of -5 bins and a smaller negative delay of -1 bin during time bins 130 to 170. For across-region and within-region dynamics, the length scales are set to = 5 and 2.5, respectively.
Note that large delays in this context are intended to repre- sent slow communication in the brain, whereas extremely large delay would imply an absence of communication be- tween regions. A small delay indicates fast communica- tion. Therefore, our data simulation is designed to reflect a scenario where region A initially has minimal effective communication with region B (delay of 5), then suddenly transmits a signal (delay of 1), followed by another period of ineffective communication (delay of 5). A delay of -5 represents communication in the opposite direction. During fitting, we set ma = 2, mh = 1, and P= 5.
Results. Figure 2(A) shows the estimated across-region communications and aame-varying delays from the test dataset (ten trials shown; within-region dynamics are pro- vided in Appendix E). The first group of communications shows a shift from slower feedback (with larger absolute delays) to faster feedback (with smaller absolute delays) starting around 200 ms after stimulus onset. In contrast, the second group exhibits a periodic pattern driven by the external drifting grating stimulus, along with a change in communication direction immediately following stimulus onset. The tire-narying delays suggest the following inter- pretation: shortly after stimulus presentation, V2 generates strong feedback signal from V2 to VI, potentially reflect- ing the emergence of surprise or prediction error (Rao & Ballard, 1999). As time progresses, both regions transition into more synchronized oscillations. Our findings on V1-V2
Results. Figure I(A) presents the estimated and truth across-region communications, temporal delays, and ground truth over T bins for P = 5. For the estimated delay, the shaded area represents the variance across five different runs.
6

---

Learning Time-Varying Multi-Region Brain Communications via Scalable Markovian Gaussian Processes
tamporal delays
ma
(A)
ma
(B)
lest log-likelhhodd
30k
(C)
30k
tog-likelihood
validation
-51
25k-
20k
estimateed delay
201
zero line
100ms
ADM
MRM-GP
DLAG
V1
(D)
16-
ADM DLAG
V2--1
V2-01
V1--22
time minutes
v2
across-region communcation
acrosa-region communcation
50
100
150
200
Figure 2. Evaluation of the ADM spike trains from V1 and V2. (A) Estimated across-region communications and time-varying temporal delays from the test dataset (ten trials are shown), showing a time-varying feedback and forward communications between V1 and V2. (B) Test log-likelihood comparison, where ADM outperforms MRM-GP and DLAG by capturing continuously tine-varying temporal delays. (C) Validation log-likelihood across different P values. (D) Computational time comparison across different length of T, demonstrating ADM's efficiency through parallel computation, outperforming both MRM-GP and DLAGG
interactions are similar to previous studies (Gokcen al., 2022; Li et al., 2024b; Gokcen et al., 2024a).
time bins (each 10 ms), and a total of 202 neurons, with ex- ternal visual stimuli comprising 4 Hz drifting gratings. Fol- lowing the approach in (Gokcen al., 2022), we first apply Factor Analysis to estimate the total number of across-region and within-region latent dynamics, determining M = 4 (see Appendix F for details). We then conduct a grid search with 5-fold cross-validation to refine the number of scross-region and within-region latent dynamics and the model order P.
Figure 2(B) compares the test log- likellhodd summed over trials and time bins, showing that our model, ADM, outper- forms MRM-GP and DLAG under the same ma and m settings. This improvement is attributed to ADM's abil- ity to capture continuously time-varying temporal delays. Figure 2(C) presents the validation log- likellhoooo across different values, with P achieving the highest value. Combined with the insights from Figures 2(A-B), this sug- gests that small P value can effectively estimate model parameters and latent variables.
Results. Figure 3(A) presents the ten estimated pairwise temporal delays from one group of across-region communi- cations. See Appendix for the estimated latent dynamics. Our results reveal consistent forward communication from VISp downstream visual areas, such as VISrl, VISal, and VISpm, aligning with the known anatomical hierarchy of the mouse visual cortex (Siegle et al., 2021). Additionally, these forward communications exhibit time-varying dynam- ics. For instance, communication between VISp and VISrI transitions from slow fast, indicating an enhanced inter- action that gradually becomes more synchronous following the initial surprise response the visual stimulus onset. In contrast, the communication between VISp and VISal shifts from fast to slow, suggesting inhibition induced by the external stimulus. Furthermore, our findings indicate that all communications involving VISam are feedback signals. This expected, as VISam is positioned at the end of the anatomical hierarchy of the mouse visual system, consistent with the anatomical hierarchy scores reported in (Siegle et al., 2021).
Figure 2(D) compares the computational time of our model with MRM-GP and DLAG on spike trains of varying T, obtained by concatenating trials. The use of parallel compu- tation significantly improves efficiency, outperforming the linear model (MRM-GP) and the cubic model (DLAG).
4.3. Five Brain Regions
In this section, we scale up our model larger neu- ral recording spanning five regions with increased time resolution. Our objective is to investigate across-region communications and identify the time-varying meso- meso-scale brain network, defined as the dynamic network spanning sAppendixub-brait regions, e.g., regions in visual cortex.
Experimental setup. We use smoothed multi-region spike train data from the Visual Coding Neuropixels project by the Allen Institute (Siegle et al., 2021), specifically from ses- sion 750749662. This dataset includes spike trains recorded from VISp. VISrI. VISal, VISpm, and VISam VISam-sub-areaa sub-areas of the mouse visual cortex. consists of 120 trials, T 200
Figure 3(B) depicts the meso-scale brain network corre- sponding to the across-region communications presented in Figure 3(A). Each node represents a region in the vi- sual system, while directed edges indicate the directional
7

---

Learning Time-Vare Vime-Varying Multi-Region Brain Communications via Scalable Markovian Gaussian Processes
ViSpm
ViSpm
ViSp
VISr
VISp
VISal
VISp
VISpm
VISp
VISam
ViSri
ViSal
ViSam
VISp
viSam
VISp
VISal
0-
temporal delays
-20
25
ViSal
ViSr
viSri
meso-scale brain network meso-scale brain network
ViSri
ViSpm
VISi
ViSam
VISal
ViSpm
ViSal
ViSam
viSpm
ViSam
when 1=3
when [=50
(C)
(D)6e5-
highest
2.1e6
22
50
20
.50
Seseoottin loa
1.966
estimated delay
zero line
100ms
20 hyperparameter combinations
ADM mDLAG
Figure 3. Evaluation of the ADM on spike trains from five visual brain regions. (A) Estimated pairwise temporal delays from one group of across-region communications, showing time-varying dynamics and alignment with the anatomical hierarchy of the mouse visual cortex. (B) Meso-scale brain networks at and 50 time bins, derived from the communications (A), illustrating changes in communication speed and direction. (C) Cross- validation results for hyperparameter selection, with m. = 3, m. 1, and P = 5 achieving the highest validation log-likelihood. (D) Test log-likelihood comparison between ADM and mDLAG, demonstrating ADM's superior fit due to its ability model time-varying communications.
communications. The length of each edge is determined by the estimated delays, reflecting the speed of communi- cation. The figure presents two meso-scale brain networks at t = 3 and = 50 time bins. The primary differences be- tween these networks include changes in the speed of certain forward communications and a direction change in the com- munication direction between VISrl and VISal, suggesting the emergence of stimulus presentation.
demonstrate that the Adaptive Delay Model (ADM) effec- tively captures these dynamics while maintaining compu- tational efficiency. Existing methods for studying across- region neural interactions can be broadly categorized into non-delay models and delay models. While non-delay mod- els provide directional communication patterns, they fail to capture temporal delays, limiting their ability to infer the communication speed. Conversely, delay models, such as DLAG and MRM-GP, introduce delay estimation but assume static discretely changing delays, which do not reflect the continuously evolving nature of brain commu- nications. Our results show that ADM overcomes these limitations by incorporating a flexible, time-varying delay mechanism, enabling more biologically relevant represen- tation of neural interactions.
Figure 3(C) presents the cross-validation results for twenty hyperparameter combinations. We first determine M = 4 using Factor Analysis (see Appendix F for details), then conduct a grid search over all combinations of ma E [0, 4],
m. E [0, 4], and P E [2, 5). The highest validation log- likelihood is achieved with m. 3. m = 1, and P 5.
Figure 3(D) compares the test log- likelihoodod summed over trials and time bins, between our model (ADM) and mD- LAG with m. = 3 latent communication channels, where mDLAG is an extension of DLAG that supports more than two brain regions using variational inference. The results indicate that ADM provides a better to the data, attributed to its ability to model varying communications. We do not compare MRM-GP and DLAG since they are lim- ited to two brain regions. Additionally, we skip time cost comparison with mDLAG because it is implemented only in MATLAB, which is significantly slower than our GPU-optimized implementation.
Neuroscience Implications. Our results from large-scale neural recordings show that across-region communication delays are not static but change over time. Notably, we observe transitions from slow feedback and forward com- munication to fast forward interactions in both datasets (Sec- tion 4.2 and Section 4.3), aligning with adaptive sensory processing in the visual cortex. These findings show the im- portance of time-varying models in capturing the dynamic nature of brain communications.
Computational Advancements. Beyond its neuroscien- tific implications, our model contributes to the broader field of computational modeling by bridging temporally station- ary GPs with SSMs. Traditional GP-SSM connections often rely on separability assumptions in spatial and temporal kernels, limiting their flexibility. Our proposed universal
5. Discussion
Summary. Our findings highlight the importance of mod- cling time-varying multi-region neural communications and
8

---

Learning Time-Varying Multi-Region Brain Communications via Scalable Markovian Gaussian Processes
connection between arbitrary temporally stationary GPs and SSMs removes this restriction. Furthermore, by leveraging parallel scan algorithms, ADM achieves impressive com- putational complexity of Ollog T), significantly improving scalability compared to existing methods.
contrastive learning problem is secretly distribution
alignment problem. Advances in Neural Information Processing Systems, 37:91597-91617, 2024.
Dowling, M., Sokót, P., and Park, I. M. Hida-mat''"em
kernel. arXiv preprint arXiv:2107 07098, 2021.
Limitations and Future Directions. Our model has cubic time complexity with respect to the number of brain regions N and the SSM order P. Although these values are typically much smaller than T, they can still become computational bottlenecks for specific cases. A potential so- lution may involve leveraging frequency domain techniques. Parnichkun et al. (2024) proposed state-free SSM with controllable canonical transition matrix, similar to ours in Eq. 11, and utilized the Fast Fourier Transform to achieve linear scaling in latent size. Similarly, Gokcen et al. (2024b) approximated the GP kernel in the frequency domain to reduce the computational cost to linear in latent size.
Dowling, M., Zhao, Y., and Park, I. M. Linear time gps for inferring latent trajectories from neural spike trains. arXiv preprint arXiv:2306.0 01800, 2023.
Glaser, J., Whiteway, M., Cunningham, P., Paninski, L., and Linderman, S. Recurrent switching dynamical sys- tems models for multiple interacting neural populations. Advances in neural information processing systems, 33: 14867-14878, 2020.
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-525, 2022.
Acknowledgement
Gokcen, E., Jasper, A., Xu, A., Kohn, A., Machens, C. K., and Yu, B. M. Uncovering motifs of concurrent signaling across multiple neuronal populations. Advances in Neural Information Processing Systems, 36, 2024a.
This work is supported by National Institutes of Health BRAIN initiative Doriiieee(100 (IU01NS131810).
Impact Statement
Gokeen, E., Jasper, A. I., Kohn, A., Machens, C. K., and Yu, B. M. Fast multi-group gaussian process factor models. arXiv preprint 2.. 2024b.
Our work presents the Adaptive Delay Model (ADM), scalable and biologically relevant framework for uncover- ing time-varying multi-region neural communications. By bridging Gaussian Processes (GPs) with State Space Models (SSMs), ADM enhances our ability to analyze large-scale neural recordings, offering insights into sensory process- ing, cognitive flexibility, and neural disorders. Beyond neuroscience, its adaptability to fields such as robotics, au- tonomous systems, and signal processing broadens its soci- etal impact. Ethically, as machine learning models increas- ingly shape neuroscience research, it is crucial to ensure their responsible application, avoiding overinterpretation of inferred neural dynamics and considering the broader impli- cations for brain-computer interfaces and neurotechnology.
Grewal, M. S. and Andrews, A. Kalman filtering: Theory and Practice with MATLAB. John Wiley & Sons, 2014.
Karniol-Tambour, O., Zoltowski, D. M., Diamanti, E. M., Pinto, I.., 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.
Li, C., Li, W. Wang. Y., and Wu, A. A differen- tiable partially observable generalized linear model with Borward-backwart message passing. arXiv preprint arXiv:2402.0126 2024a.
References
Li, W., Li, C., Wang, Y., and Wu, A. Multi-region marko- vian gaussian process: An efficient method to discover directional communications across multiple brain regions. In Forty-first International Conference on Machine Learn- ing, 2024b.
Balzani, E., Noel, P., Herrero-Vidaa Vidal, P., Angelaki, D. E., and Savin, C. A probabilistic framework for task-aligned intra-and inter-area neural manifold estimation. arXiv preprint wwxXiv:2009.281120 2022.
Lillicrap, T. P., Santoro, A., Marris, L.., Akerman, C. J., and Hinton, G. Backpropagation and the brain. Nature Reviews Neuroscience, 21(():3355344 346, 2020.
Blelloch, G. E. Prefix sums and their applications. 1990.
Boots, B. Learning stable linear dynamical systems. On- line). Avail.: htlp:://www: ml. cmu. edu/reseurch/dap- papers/dap_bowd boots. pdf [Accessed 30 05 2016), 2009.
Liu, Y. A., Nong, Y., Feng, J., Li, G., Sajda, P., Li, Y., and Wang, Q. Phase synchrony between prefrontal noradren- ergic and cholinergic signals indexes inhibitory control. bioRxiv, pp. 2024 2025.
Chen, Z., Lin, C.-H., Liu, R., Xiao, J., and Dyer, E. Your

---

Learning Time-Varylln Multi-Region Brain Communications via Scalable Markovian Gaussian Processes
Loper, J., Blei, D., Cunningham, J. P., and Paninski, L. A general linear-time inference method for gaussian pro- cesses one dimension. Journal of Machine Learning Research, 22(234):1 2021.
Ulrich, K. R., Carlson, D. E., Dzirasa, K., and Carin, L. Gp kernels for cross-spectrum analysis. Advances in neural information processing systems, 28, 2015.
Wang. Y., Wu, Z.,, Li, C., and Wu, Extraction and re-
Nejatbakhsh, A., Geadah, V., Williams, A. H., and Lip- shutz, D. Comparing noisy neural population dynam- ics using optimal transport distances. arXiv preprint arXiv:2412.144211 2024.
covery of spatio-temporal structure in latent dynamics alignment with diffusion models. Advances in Neural Information Processing Systems, 36:38988-39000 2023.
Wang. Y., Li, C., Li, W., and Wu, Exploring behavior- relevant and disentangled neural dynamies with genera- tive diffusion models. Advances in Neural Information Processing Systems, 111-4722-34736, 2024.
Neumaier, and Schneider, Estimation of parameters and eigenmodes of multivariate autoregressive models. ACM Transactions Mathematical Software (TOMS), 27(1):27-57, 2001.
Williams, A. H., Kunz, E., Kornblith, S., and Linderman, S. Generalized shape metries on neural representations. Advances in Neural Information Processing Systems, 34: 2021.
Parnichkun, R. N., Massaroli, S., Moro, A., Smith, J. T., Hasani, R., Lechner, M., An, Q., Ré, C., Asama, H., Ermon, S., et al. State-free inference of state-space models: The transfer function approach. arXiv preprint mrXiv:2405.06147, 2024.
Wilson, A. and Adams, R. Gaussian process kernels for pattern discovery and extrapolation. International conference on machine learning, pp. 1067-1075. PMLR, 2013.
Parra, G. and Tobar, F. Spectral mixture kernels for multi- output gaussian processes. Advances in Neural Informa- lion Processing Systems, 30, 2017.
Zandvakili, A. and Kohn, A. Paired vI-v2 neuronal spiking responses in anesthetized macaque monkey. CRCNS, org. 2019.
Rao, R. P. and Ballard, H. Predictive coding in the visual cortex: a functional interpretation some extra-classical ficeptive-field effects. Nature neuroscience, 2(1):79-87, 1999.
Zhang, Y., Wang, Y., Jiménez-Benetóó, D., Wang, Z., Az- abou, M., Richards, B., Tung, R., Winter, O., Dyer, E., Paninski, L., et al. Towards a" universal translator" for neural dynamics at single-cell, single-spike resolution. Advances in Neural Information Processing Systems, 37: 80495-80521, 2024.
Safaie, M., Chang, J. C., Park, J., Miller, L. E., Dudman, T., Perich. M. G., and Gallego, J. Preserved neural dynamics across animals performing similar behaviour. Nature, 623(7988):765-7777 2023.
Särkkä, S. and FeréádzFer,ández, Á. F. Temporal paral- lelization of bayesian smoothers. IEEE Transactions on Automatic Control, 66(1):299-306, 2020.
Zhao, Z. State-space deep gaussian processes with applica- tions. arXiv preprint arXiv:2111. 12604, 2021.
2020.
Semedo, J. D., Zandvakili, A., Machens, C. K., Byron, M. Y., and Kohn, A. Cortical areas interact through communication subspace. Neuron, 102(1):249-255 2019.
Siegle, J. H., Jia, X., Durand, S., Gale, S., Bennett, C., Graddis, N., Heller, G., Ramirez, T. K., Choi, H., Luviano, A., al. Survey of spiking in the mouse visual system reveals functional hierarchy. Nature, 592(7852):86- 2021.
Solin, A. et al. Stochastic differential equation methods for spatio-temporal gaussian process regression. 2016.
Steinmetz, N. A., Aydin, C., Lebedeva, A., Okun, M., Pa- chitariu, M., Bauza, M., Beau, M., Bhagat, J., Böhm, C., Broux, M., et al. Neuropixels 2.0: A miniaturized high-density probe for stable, long-term brain recordings. Science, 3726559))eahf4588 2021.
10

---

Learning Time-Varying Multi-Region Brain Communications via Scalable Markovian Gaussian Processes
A. Derivation for VVT, WV and WWT
Let's inspect VVT RNNPNN first. To simplify the notation, we use æ to represent o Eq. 7. have:
22
-T-P
a 2
22
43
TT-
32
3g
2 P+1
VVT
aP
3P+1
UT
TT
2T
T--P1
TT-1
(15)
+ ... + 2T -T- -T--PTT- PUT-P
310p 1 + 22 p+1
.... + WT- PT-PPT
pppr ..... aT -T-1
pppp + 1p410 + 2T 2T-13T-1
where the first element I o In other words, since æ is modeled as a stationary GP, the elements æ 1 TT- - are all equivalent and correspond to the diagonal elements ((0) in Eq. 4. Similarly, the elements a ..... 2T - 1- represent cross-covariances with time interval P - 1, which correspond to the off-diagonal elements K(P 1). Therefore, we can further write Eq. 15 as:
NxN represents the auto covarricv of æ1, which is essentially the kernel ( (Eq. 4),
((0) + ... K(0)
K(-P+1)) K(-P-t 11]
VV
K(P - 1) + K(P-1) -
K(0) + K(0)
K(0)
K(-1)
K(-P+D 1)
(16)
K(1)
K(0)
K(- - P 2)
x
K(P-1)
K(P- -2) 2)
K(0)
Following the same way, we can also represent WVT NXNP and WWT E NNN using K:
WVT [(P) K(P-1) WWT & K(0).
K(1)]
(17)
If the computation of G in Eq.8 leads to numerical issues because VVT has singular values that are nearly zero, more numerically stable approach to rewrite VVT by Cholesky factorization:
VVT
VWT
L
D
= LLT, T
L=
VV = LL La
(18)
WVT
WWT
L2
where D E [N((HI)XN(P-)) WV L.L. WWT K(0). and L1 NPXNP, L2 E NXMP, L3 E R the sub-matrices of L. In practice, Eq. 18 factorizes D + öI with small postive number å to ensure the positive definite of D. Then, the estimation for G can be cast in the form of L and the measurement matrix Q is the residual covariance of residual R:
G (VV LgLi'
Q ( GV) GVV
Lg
(19)
T-P-1
T-PP1

---

Learning Time-Varying Multi-Region Brain Communications via Scalable Markovian Gaussian Processes
B. Details for Kalman EM
consider the ShnearGGaussint state-space model
4 A + Wt,
wi ~ M( Qt),
y: = Hæx + ve.
Vi ~ ~M0, V),
where the transition matrix A, and process-noise covariance Q change with time. Let = {AITT QU. H, V) denote the full parameter set.
E-step
With parameters fixed at O² the expected complete-data log-likelihood
Q(+ )*) = Ex 0n logglops p(x, ())
2 yi V y: V + TH( (VV.))
Ti(Q cc)) - A. + THA, AA,CI.)))
T
T
log
VI
Σ
log Q
(20)
2
2
t=1
Here
x,- E(a, y, O"), C. C( y, O*,, 1.--1-Eaaaa__ I ly, O*).
which are obtained with Kalman filter followed by a Rauch-Tung- SStiieel smoother run with the time-varying parameters (Boots, 2009). The individual expectations that appear in (20) expand to
B[x, H V
Hæ:
x H
H + V
(21)
Biz Q z1 = x Q x + Qr C.).
(22)
E| B|z Q A A = Q, A, + TIQ A, A,C,,,
(23)
E| Bi Qi A - =I A Q A, + TT(A Q, - A,C,.-)).
(24)
M-step
The M-step maximises (20) with respect to O
@ arg max Q(e 6*),
where the paramerer sets θ is updated either via gradient descent or in closed form, depending on the specific component.
12

---

Learning T'ime-Vayying Multi-Region Brain Communications via Scalable Markovian Gaussian Processes
C. Generation Samples
To better understand the effect of different values, we generate samples with T' 200 time bins from our model using various values. Figure 4 shows that when is very small (e.g., P 1), the generated samples appear unsmooth. However, for P 2, the generated samples exhibit noticeable visual differences.
M
W
P=1
P=2
Pegion region
rightr region
P=3
P=4
Figure Generated samples from our model with MOSE kernel when P 1,2,
13

---

Learning T'me-Varrinng Multi-Region Brain Communications via Scalable Markovian Gaussian Processes
Table 1. MSE for GP regression with single-output kernels.
Reg-MSE/ 10-1
Exp
Matern 3/2
SE
RQ
SM
GP
5.7 = 0.1
5.9 - 0.2
3.1 - 0.1
3.0 0.1
3.0 + 0.2
SSM-Approx
5.9 0.1
6.2 0.1
3,3 0.1
3.4 = 0.1
3,3 + 0.2
Table 2. MSE for GP regression with multi-output kernels.
Reg-MSE/
10-1
MOSE
MOSM
CSM
LMC
GP
7.4 + 0.02
7.5 + 0.02
8.2 + 0.05
0.66 + 0.02
SSM-Approx
7.6 + 0.04
7.9 0.08
7.7 + 0.09
0.72 + 0.02
D. Gaussian Process Regression
To verify the universal connection between arbitrary temporally stationary Gaussian Processes (GPs) and State Space Models (SSMs), compare GP regression performance using our SSM approximation and the standard GP. We generate samples of 300 points from GP with added Gaussian noise regression data. The samples then randomly split into training (Arainn Mwain) and testing (heat Mhea) sets, with 60% used for training and 40% for testing.
The kernels we evaluated are:
Exponential (Exp): outplt with K(t,t') = o² exp
(-1b(1)
Matern 3/2 (Matern): Single-output with ((I,'") = o² (1+ vālt-11
V3|t-1
exp
Squared Exponential (SE): Singlc-output with K(t,t') = =² exp
=eexxp
(
20
Rational Quadratic (RQ): Single-output with K(ttt') = 0² (I+ (1-t"))
2ntt
Spectral Mixture (SM) (Wilson & Adams, 2013): Single-output with
k(t,t')=
Σg og exp
(t-2") cos (t- ((ff-t)))
Multi-Output Squared Exponential (MOSE) (Gokcen et al., 2022): Multi-output with Kij(1, ()
=
(
exp
2
Multi-Output Spectral Mixture (MOSM) (Parra & Tobar, 2017): Multi-output with tt,) = Eg oijn exp (t-t'+55,.)) cos (Wij.g(t - t') + Φij.g)-
212
Cross-Spectral Mixture (CSM) (Ulrich et al., 2015): Multi-output with
(iy(t,t')= = I T exp
(t-+')²)
cos (Wig.n(t - 1') + P.
Linear Model of Coregionalization (LMC): Multi-output with K(t,t') Eg Bg & (t,t'))) where B, is coregion-
alization matrix and kg (1,1") (') is single-output kernel.
The number of orders P for each kernels are as follows:
Exponential (Exp): P = 1.
Matern 3/2 (Matern): P= 2.
Squared Exponential (SE): P=
Rational Quadratic (RQ): 4.
Spectral Mixture (SM): 2
Multi-Oulput Squared Exponential (MOSE): / = 2.
Multi-Output Spectral Mixture (MOSM): 2.
Cross-Spectral Mixture (CSM): P 4.
14

---

Learning Time-Varying Multi-Region Brain Communications via Scalable Markovian Gaussian Processes
Linear Model of Coregionalization (LMC): When (ig (') is SE kernel, P P 2.
The results shown in Table and Table 2, where our SSM approximation demonstrates regression performance comparable to GP in terms of MSE.
15

---

Learning Time-Varying Multi-Region Brain Communications via Scalable Markovian Gaussian Processes
E. Additional Across- and Within-Region Latent Variables
Figure presents within-region neural activity for both synthetic data and V1-V2 neural spike trains.
(A)
region ... true
region 2 ... true
within-region activity for synthetic data
(B)
V1
V2
within-region activity
within-region activity 2
(C)
VISp
VISrl
VISal
VISpm
VISam
across-region activity 1,2,3
within-region activity
Figure 5. (A) Within-region neural activity for synthetic data. (B) Within-region neural activity for V1-V2 neural spike trains (ten trials are shown). (C) The across- within-region latent variables for the five-region dataset (ten trials are shown).
16

---

Learning Time-Varying Multi-Region Brain Communications via Scalable Markovian Gaussian Processes
F. Factor Analysis for Five Regions Spike Trains
Figure presents the Factor Analysis results for neural spike trains from five regions, and we select the latent size to be largest optimal latent size across five regions, which is
ViSp
-18.8
-20A
-19.0
28.5
-003
L05 0p Beeteehooo ioodd
28.5
-1.0
-994
-20
-99.6
41.2
41.3
-100.2
-29
10
10
10
Istent sze
Istent s2e
Iatent S2e
Vi5pm
ViSam
-53.1
53.2
-59.0
59.2
-335
5994
536
59.8
-598
-339
10
2
5
6
7
8
9
10
Iatent 926
latent 52c
Figure 6. The figure presents Factor Analysis results for neural spike trains from five regions. The optimal latent size, determined as the maximum across all regions, selected yielding highest test log likelihood.
17

---

Learning Time-Varying Multi-Region Brain Communications via Scalable Markovian Gaussian Processes
G. More Synthetie Data Experiments
(A)
delays between region 1 and 2
delays between region and 3
delays between region and 4
delays betwwenn region
delays between region 2 and 3
temporal delays (time bins
25
delays between region 2 and 4
delays between region 2 and 5
delays between region and
delays between region and 5
delays between region 4 and 5
M0
estimated delays ground truth
Figure 7. Estimated Pairwise Temporal Delays for Synthetic Data with Five Regions. Dashed lines indicate the ground truth delays, red lines show the estimated delays, and shaded areas represent the variance across different runs. The results demonstrate that the model accurately recovers the true delays. The increased MSE observed in Figure 8(A) is attributed to amplitude variability, which is not practical concern, as temporal delay patterns are well preserved.
(B)
0.04
1.0
0.04
1.0
0.9
0.9
0.03
Pearson's CC
0.03
Pearson à CC
MSE
0.02
0.8
MSE
0.02
0.8
0.7
0.7
0.01
0.6
0.01
0.6
2 regions
regions
regions
regions
0.5
regions
3 regions
4 regions
5 regions
0.0
lstents
latents
6 latents
10 latents
latents
latents
latents
10 latents
(C)
0.04
(D)
500k
0.9
0.03
Pearson's DC
400k
test log-llkellhhce
MSE
300k
0.02
200k
0.01
100k
0.0
200
400
800
200
400
600
800
under specified under-speciled
correctly specified
over- veerapecifled
time points
time points
Figure 8. Additional Results for Synthetic Dataset. (A) Model evaluation using MSE and Pearson's correlation coefficient (CC) between estimated ground truth delays across different numbers of regions. While MSE increases to greater amplitude variability with more regions, CC remains stable, indicating reliable recovery of temporal delay patterns. See also Figure 7 for a visualization of estimated delays with five regions. (B) Model evaluation across varying numbers of latent variables, showing stable MSE and CC performance. (C) Model evaluation under different data lengths. Longer sequences yield lower MSE, likely due improved estimation of underlying delays with more data. (D) Test log- likeliieldood comparison when the number latent variables is under-specified, correctly specified, or over-specified. Both under- and over-specification result lower log- likelihood and higher variance.
18

---

Learning Time-Vore Varying Multi-Region Brain Communications via Scalable Markovian Gaussian Processes
(A)
500k
lest "gglllellhood
450k-
400k
350k
ADM
MRM-GP
DLAG
mp-srSLOS
Figure 9. Test LokeLikelihood Comparison of ADM, DLAG, MRM-GP, and mp-srSLDS on the Synthetic Dataset used in Section 4.1. ADM consistently outperforms other methods by effectively capturing continuously time-varying temporal delays.
19

---

Learning Time-Varying Multi-Region Brain Communications via Scalable Markovian Gaussian Processes
H. Additional Results for Decoding Visual Stimuli in the Five-Region Brain Dataset
(A)
0.8
test disssification accuracy
0.6
0.4
0.2
all data
latent variable
across
latent variable
within
Figure 10. Decoding Visual Stimulus Orientation from the Visual Rostrolateral Area (VISrl). We analyze neural data from five regions and extract the learned latent variables. A linear decoder then used classify the orientation of visual stimuli (0", 90", and 135%) presented to the mouse during data collection. We evaluate decoding performance using three inputs: (I) raw observed neural activity from Visual Rostrolateral Area, (2) across-region latent variables (representing the communication subspace) Visual Rostrolatcral Arca, and (3) within-region latent variables of Visual Rostrolateral Area. The results show that decoding directly from the observed neural data yields the highest classification accuracy. Among latent spaces, the communication subspace achieves higher accuracy than the within-region subspace. Notably, both the observed data and the communication subspace perform above random guessing. indicating that orientation information is preserved communication subspace of Visual Rostrolateral Area, which a region known to he involved in motion and spatial processing. Finally, the drop in accuracy from the observed data to the communication subspace is likely due to information loss from dimensionality reduction.
20

---

Learning Time-Varying Multi-Region Brain Communications via Scalable Markovian Gaussian Processes
I. Related Work on Neural Alignment and Task-Dependent Variability Modeling in Multi-Region Anslysis Interaction
Multi-region neural analyses increasingly rely on tools that can jointly characterise activity patterns within each area and the structured co- variability that links them. Williams et al. 2021 gives a geometric foundation for such comparisons by turning pprreentationallssmiaarit heuristics into proper metric spaces; they demonstrate that these "shape metrics" scale to surveys of 48 distinct visual areas in the Allen Brain Observatory, enabling distance-based clustering and embedding of whole-brain activity patterns. Safaie al. 2023 pushes the multi-region perspective across species, showing that low-dimensional trajectories extracted from motor cortex recordings in both monkeys and mice align common latent space and remain predictive of movement individual brains. Complementing these observational studies, Balzani et al., 2022 introduces TAME-GP, probabilistic manifold model that factorises population variability into within-region and across-region components, imposes Gaussian- process priors for temporal smoothing, and infers single-trial latent dynamics that reveal task-dependent communication subspaces among multiple brain areas
in the absence of overt behaviour, pointing conserved circuit-level dynamics that transcend
21


Reviews:
[{"review_id": "nqbN9UXp1C", "reviewer": "Reviewer_Wfpo", "summary": "This submission describes an approach for inferring latent factors underlying shared neural responses across brain areas. Notably, the approach enables inferring a continuous and time-varying delay factor that captures the temporal delay between two brain areas. The submission formulates this as both a Gaussian Process model and state-space model (SSM) and leverages the SSM setup to enable faster inference via parallel scans. The method is validated in a simulated dataset and demonstrated on two experimental applications for the analysis of visual data. Importantly, the method found time-varying delay responses and was applied to more than two brain areas.", "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": "Dear Reviewer Wfpo,\n\nThank you for the encouraging feedback and practical suggestions! We have made several clarifications according to your questions and comments. Hopefully these will resolve most of your concerns, and that they can be taken into account when deciding the final review score.\n\n> Figure 3(B) is challenging to understand. I suggest presenting these results in a different format, potentially via a matrix where each element represents a signed pairwise connection.\n\n- We thank the reviewer for pointing this out. Figure 3(B) is intended to highlight the time-varying nature of temporal delays across five brain regions. It visualizes the learned temporal delays from Figure 3(A) at time points $t=3$ and $t=50$. In the figure, the orange arrows indicate the direction of communication: a positive delay denotes forward communication, while a negative delay indicates feedback communication. Additionally, the length of each edge represents the absolute value of the corresponding temporal delay. \n\n- To provide a clearer and more quantitative view, we will add a matrix alongside the network diagram to explicitly show both the magnitude and sign of the learned temporal delays.\n\n> Can the authors provide additional details on the formulation and fitting of the time-varying delay? It appears to be one of the crucial innovations. I'm curious on how flexible this model is and whether it presents any issues in fitting. The paper states that the resulting $\\mathbf{\\hat{A}}_t$ are still constrained at each time point via a shared length scale. However, it is not obvious to me how strongly that constrains these parameters.\n\n- At each time step $t$, we construct a time-specific Markovian GP conditioned on the MOSE kernel corresponding to that time step. The state transition matrix $\\mathbf{\\hat{A}}_t$ is uniquely determined by the delay $\\theta_t$ and length scale $l$ as mentioned in Sec. 2.2.  While the delay $\\theta_t$ is allowed to change over time, the length scale $l$, which is shared across all time steps, limits the freedom of $\\mathbf{\\hat{A}}_t$. This shared length scale fully determines the temporal dynamics of the latent variables and thus constrains the evolution of dynamics from $\\mathbf{\\hat{A}}_t$ to follow the same smoothness. In practice, this allows the model to flexibly capture time-varying delays while maintaining similar temporal dynamics across time.\n\n- The blue and red latent variables shown in Figure 1(A) illustrate this behavior: they share overall temporal dynamics but exhibit different delays at each time point. Empirical results in Sec. 4.1 further confirm that this shared length scale leads to robust and smooth delay estimates without fitting issues if a reasonable initialization is used, e.g., setting a relatively large length scale to encourage initial smoothness. This initialization aligns with standard practice in GP regression, where initializing with a larger length scale helps guide optimization during early iterations.\n\n> Additionally, Appendix B does not appear to describe the time-varying dynamics, as $\\mathbf{A}$ is not indexed by in the EM equations. Is there a different objective for the time-varying formulation? If so, it would be very helpful and important to put that in the appendix as well.\n\n- We apologize for the confusion regarding Appendix B. In our time-varying extension, the underlying EM objective remains the same: maximizing the expected complete-data log-likelihood, but the state-space model is augmented to include time-specific transition matrices $\\mathbf{\\hat{A}}_t$ and process noise covariances $\\mathbf{\\hat{Q}}_t$. The derivations for the EM updates are analogous to those in Appendix B, with the key difference that the E-step now involves running a vectorized Kalman filter and smoother over the entire time sequence. As a result, the shapes of $\\mathbf{\\hat{A}}$ and $\\mathbf{\\hat{Q}}$ become $(T, NP, NP)$. We will revise Appendix B to explicitly incorporate these modifications and clearly describe the time-varying dynamics.\n\nAgain, we thank the reviewer for the comments, constructive feedback, and insightful suggestions, which significantly improved the quality of our paper. If you have further questions, we are happy to discuss them!", "claims_and_evidence": "The submission presents convincing evidence that the approach enables inference of time-varying delays in brain communication.", "methods_and_evaluation_criteria": "Yes, the datasets and methods make sense. The submission evaluates the method using standard approaches on both simulated datasets and two relevant neuroscience datasets.", "theoretical_claims": "N/A", "experimental_designs_or_analyses": "Yes, I checked the details for the simulated and two experimental data analyses.", "relation_to_broader_scientific_literature": "This submission is related to a line of important work developing methods for the analysis of neural recordings that span multiple brain areas. In particular, this submission builds on previous work on Gaussian process models for capturing shared and private variability across brain areas, potentially with time delays.", "essential_references_not_discussed": "N/A", "other_strengths_and_weaknesses": "Strength:\n\nThe submission nicely demonstrated the generality of their proposed approach, showing how under various GP kernels they could relatively accurately approximate a GP regression using their SSM approximation. \n\nWeakness:\n\nFigure 3(B) is challenging to understand. I suggest presenting these results in a different format, potentially via a matrix where each element represents a signed pairwise connection.", "other_comments_or_suggestions": "N/A", "questions_for_authors": "Can the authors provide additional details on the formulation and fitting of the time-varying delay? It appears to be one of the crucial innovations. I'm curious on how flexible this model is and whether it presents any issues in fitting. The paper states that the resulting $\\hat{A}_t$ are still constrained at each time point via a shared length scale. However, it is not obvious to me how strongly that constrains these parameters. \n\nAdditionally, Appendix B does not appear to describe the time-varying dynamics, as $A$ is not indexed by $t$ in the EM equations. Is there a different objective for the time-varying formulation? If so, it would be very helpful and important to put that in the appendix as well.", "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": "hOVWwNcH8O", "reviewer": "Reviewer_9Ksf", "summary": "Modeling neural activity across networks of populations of neurons across multiple brain regions is critical to understand neural computation and how information is processed. While the recent recoding technological advances have made it possible to acquire the data, modeling tools are limited in their ability to capture all the variability present in it. To address some of these limitations, the authors present a new model that extends latent space models to capture time varying variability across brain regions while allowing for temporal delays between them. The authors also introduced an inference procedure to reduce computational complexity. They test the model in a synthetic dataset and two neural datasets showing the feasibility of the approach.", "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": "Dear Reviewer 9Ksf,\n\nThank you for the constructive comments! We have made several improvement according to your questions and comments. We hope these sufficiently clarified your concerns, and that they can be taken into account when deciding the final review score.\n\nThe additional results are provided (https://anonymous.4open.science/r/rebuttal-figures-for-ICML-2025-42E6/rebuttal_figures.pdf), including:\n\n- A comparison of our model with existing methods on the synthetic dataset (Rebuttal-Figure 4(C)).   \n- The model performance vs. the number of regions, latents, and length on the synthetic dataset (Rebuttal-Figure 1(A-C)).  \n- Decoding visual stimulus-related information from the learned latents (Rebuttal-Figure 3).\n\n> The aforementioned references should be added (Williams et al. NeurIPS 2021, Safaie et al. Nature 2023, Balzani et al. ICLR 2023).\n\n- Williams et al. (NeurIPS 2021) introduce a metric framework for comparing static neural representations, Safaie et al. (Nature 2023) show preserved latent dynamics across subjects using alignment, and Balzani et al. (ICLR 2023) propose a task-aligned model capturing within- and across-region neural variability. In contrast, our work focuses on modeling delay-based brain communication across regions, aiming to capture dynamic, time-varying delays that mediate across-region interactions. This enables insights into the directionality and strength of functional coupling.\n\n- Therefore, we respectfully disagree that our work should be directly compared to alignment-based approaches such as Williams et al. (NeurIPS 2021) and Safaie et al. (Nature 2023), as our modeling goals differ. On the other hand, we greatly appreciate the novel perspective introduced by Balzani et al. (ICLR 2023) on task alignment within communication subspaces. We will include a discussion of their work in our revised manuscript to clarify how it complements and contrasts with our approach.\n\n> They minimally evaluate the proposed model in a synthetic dataset, but the authors fail to show comparative performance with other models in this data.\n\n- We thank the reviewer for raising this concern. In Rebuttal-Figure 4(C), we present the test observation log-likelihoods on the time-varying synthetic dataset used in Sec. 4.1. Our model outperforms existing multi-region modeling methods, demonstrating its effectiveness in capturing time-varying multi-region communications.\n\n> How does the model behave with respect to the latent dimensionality, or number of areas, or length of recording?. \n\n- We thank the reviewer for this valuable suggestion. In Rebuttal-Figure 1(A-C), we evaluate our model using MSE and Pearson's correlation coefficient (CC) between estimated and ground truth delays under varying numbers of regions, latent dimensions, and lengths. The model shows stable performance across conditions. While MSE increases due to amplitude variability with more regions, CC remains stable, indicating reliable recovery of temporal delay patterns. As shown in Rebuttal-Figure 2, this amplitude variability is not a practical concern because the temporal patterns, which are crucial for understanding brain communication, are well preserved.\n\n> The method uses vanilla FA to estimate the latent dimensionality, ..., why not use likelihood estimates from the model itself?\n\n- We follow the strategy from prior work (Gokcen et al., 2022), see Sec. 4.3, using FA to estimate total latent dimensionality and then performing a grid search with model log-likelihood to determine across- and within-region dimensions. This approach balances accuracy with computational efficiency, as a full grid search of latent dimensionality is time-consuming.\n\n> Axes ticks and labels are missing in multiple figures.\n\n- We apologize for the oversight and thank the reviewer for pointing it out. We will add axis ticks and labels in the revised version.\n\n> Possible extensions would be to include comparisons of the latent representations across models, or decoding performance with respect to relevant neural computation variables, such as visual stimulus class.\n\n- Thank you for the suggestion. We have visualized our model’s latent representations on both neural datasets in Figure 2(A) and Appendix E. In the revised version, we will include latent representations from other models for comparison. Additionally, Rebuttal-Figure 3 demonstrates that decoding from our latent variables shows stronger task-related information (e.g., grating orientations) in the communication subspaces compared to the within-region subspaces. We also view the discovery of further task-relevant information in these subspaces as a promising direction for future work.\n\n> Different dimensionality across different latent spaces?\n\n- Yes, we allow the number of within-region latents to vary over regions.\n\nAgain, we thank the reviewer for the insightful suggestions, which improved the quality of our paper. If you have further questions, we are happy to discuss them!", "claims_and_evidence": "The authors clearly frame their work and provide theoretical proof and empirical results to back them. The authors present results for one synthetic dataset and two neural datasets, showing the feasibility of their model. However, the authors motivate their work by stating that existing methods fail to capture neural representations that could provide additional insight into neural processing, but they fail to present novel neuroscience discovery beyond corroboration with results surfaced by existing models. Comparatively showing how this model can improve interpretability or drive new insight is critical to understand the impact of the contributed work.", "methods_and_evaluation_criteria": "The methods are well motivated and described. The authors compared their model to alternative existing models that explicitly model temporal dynamics and delays, but do not provide systematic comparisons to other static methods such as vanilla Procrustes alignment (Williams et al NeurIPS 2021, Safaie et al. Nature 2023). Including such comparisons would be critical to evaluate the tradeoffs between computational cost, data demands and explanation capabilities across models. \nThey minimally evaluate the proposed model in a synthetic dataset, but the authors fail to show comparative performance with other models in this data. Moreover, they could have further used this setting to assess the applicability and limitations of the model (e.g. how does the model behave with respect to the latent dimensionality, or number of areas, or length of recording?). The authors present results on two sets of neural data, where they prove that their method is more computationally efficient and provides efficient reconstruction of the neural responses and captures the temporal delays between regions. It would be relevant to assess not only the delays and directionality of the connections between regions, but also i) the neural representations in the communication space and ii) the task relevant information present between the communications subspaces. This would further highlight the potential of the model to drive neuroscientific discovery.", "theoretical_claims": "The theoretical claims are correct, and full proofs are shown in the supplementary material.", "experimental_designs_or_analyses": "As mentioned previously, the design and analysis address the feasibility of the approach but lacks to fully illustrate the advantages of the model with respect to alternative models or the ability to drive new scientific insight. Including comparisons to simpler methods, which are expected to fail to capture the neural variability, and a comparison beyond log likelihood estimates would further strengthen the presented results. Possible extensions would be to include comparisons of the latent representations across models, or decoding performance with respect to relevant neural computation variables, such as visual stimulus class.", "relation_to_broader_scientific_literature": "The authors adequately motivate their work and present alternative methods. However, they could include references to other simpler neural alignment methods (Williams et al NeurIPS 2021, Safaie et al. Nature 2023). And add a references to methods that capture communication spaces between brain regions and also model the task variability, such as Balzani et al. ICLR 2023.", "essential_references_not_discussed": "The aforementioned references should be added to illustrate alternative alignment methods for neural recordings without explicit modeling of temporal dynamics (Williams et al NeurIPS 2021, Safaie et al. Nature 2023) or that capture task and behavioral information (Balzani et al. ICLR 2023).", "other_strengths_and_weaknesses": "The method uses vanilla FA to estimate the latent dimensionality, and it is not clearly listed as a limitation or design choice. Why not use likelihood estimates from the model itself?", "other_comments_or_suggestions": "Axes ticks and labels are missing in multiple figures, including Fig. 1A, 4 and 5 both x- and y- axes, Fig1C, 2A,D , 3A, 6 x-axis.", "questions_for_authors": "Does the model support different dimensionality across different latent spaces?", "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": "w7iyTILP8Z", "reviewer": "Reviewer_MsHr", "summary": "This paper extends GP-based methods for modeling multi-region neural recordings to the case where temporal delays in communication can dynamically shift. This is achieved by a novel combination of GPs with state-space models. The authors evalute their work on both synthetic and real multi-region neural recording datasets.", "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": "Dear Reviewer MsHr,\n\nThank you for the encouraging feedback and suggestion! We have made several clarifications according to your questions and comments. We hope these sufficiently clarified your concerns, and that they can be taken into account when deciding the final review score.\n\nThe additional results are provided (https://anonymous.4open.science/r/rebuttal-figures-for-ICML-2025-42E6/rebuttal_figures.pdf), including:\n\n- Delay estimation on the synthetic dataset with five regions (Rebuttal-Figure 2).\n\n> It would be interesting to see how well the model reconstructs true delays in the setting where there are more than two regions though.\n\n- We thank the reviewer for this insightful suggestion. In Rebuttal-Figure 2, we present the estimated time-varying delays from synthetic data with five regions. Rebuttal-Figure 1(A-C) further shows model performance in terms of MSE and Pearson’s correlation coefficient (CC) between the estimated and ground truth delays as the number of regions increases. While MSE rises due to greater amplitude variability, CC remains stable, indicating reliable recovery of the temporal delay patterns. As shown in Rebuttal-Figure 2, this variability in amplitude is not a practical concern, as the temporal patterns, which are crucial for understanding brain communication, are well recovered. We will add these results in the revised version.\n\n> To what degree is the inference of delays vs. communications an ill-posed problem? Is there a degeneracy where the true message could have static delays in terms of the transmission channel but it looks like a time-varying delay due to the time-varying nature of the communications?\n\n- We thank the reviewer for raising this question. Separating the effects of time-varying delays from changes in the message content is inherently challenging. In principle, the inference problem is ill-posed because the same observed data could be explained either by static delays combined with time-varying messages or by genuinely time-varying delays. In other words, there is a potential degeneracy: even if the true communication channel has static delays, the variability in the messages may cause the inference process to interpret these as time-varying delays.\n\n- To mitigate this issue, our model incorporates a regularization strategy by using a shared length scale parameter $l$ across all time points. As described in Sec. 2.2, the state transition matrix $\\mathbf{\\hat{A}}_t$ is uniquely determined by the delay $\\theta_t$ and the length scale $l$. By sharing $l$ across time, we constrain the temporal evolution of the delay parameters, enforcing similar temporal dynamics and reducing the risk of misattributing variability in the messages to changes in delays. Furthermore, our synthetic experiments with known ground truth demonstrate that the model reliably recovers the true dynamic delays, indicating that the inferred time-varying delays reflect genuine changes in the communication pathways rather than artifacts arising from the message content.\n\n- However, we acknowledge that further research could explore additional constraints or independent measures to better address this problem, but our current results suggest that the degeneracy is minimal in the context of our application.\n\nAgain, we thank the reviewer for the comments, constructive feedback, and insightful suggestions, which improved the quality of our paper. If you have further questions, we are happy to discuss them!", "claims_and_evidence": "Overall, the paper did a good job supporting its claims with a thorough and extensive evaluation. They are able to reconstruct the true time-varying delays in the setting of synthetic data. And they show that they can fit real-world data, and the results make sense in the context of known details about the dynamics of the studied regions.", "methods_and_evaluation_criteria": "The proposed methods and evaluation criteria (both benchmarks and metrics) make sense for this problem setting. The real-world data analysis is more of a proof-of-concept, which is appropriate for the venue.", "theoretical_claims": "I looked over the derivations in Appendix A and nothing popped out to me as incorrect.", "experimental_designs_or_analyses": "I looked over the experimental design of the main results and found them to be sound. The synthetic data with known time-varying temporal delays is a great place to start as a sanity check. And the real-world datasets show the scalability of the method to more than two regions. It also demonstrates that it fits to real data, and illustrates the type of hypotheses it can help generate for further testing. It would be interesting to see how well the model reconstructs true delays in the setting where there are more than two regions though. This could be tested with a synthetic model similar to the proposed one, but with more than two regions. Additionally, the results show that the model is fairly consistent across seeds, which is a nice property for a modeling framework.", "relation_to_broader_scientific_literature": "This paper extends state-of-the-art GP-based multi-region brain modeling techniques which focus on modeling inter-region delays in communication. Prior work primarily focuses on static delays, but this method considers the case where they are influenced by cognitive states. Their developed method also combines GPs with state-space models in an interesting way which may be of indepdent interest outside of neuroscience.", "essential_references_not_discussed": "There are no essential references missing from this paper as far as I'm aware of the literature.", "other_strengths_and_weaknesses": "While not the primary concern of this paper, it is unclear how well this framework can account for the ill-posedness of inferring time-varying messages with static delays versus time-varying delays. It's possible that some of what appears to be time-varying in the delay of data is more due to the time-varying nature of the messages than something about the actual communication channel. But this comes down to the interpretation of the model, which is an important but different problem than the focus of the paper.", "other_comments_or_suggestions": "No additional comments.", "questions_for_authors": "- How well does the model reconstruct true time-varying delays in data generated by a synthetic model with more than two regions? \n- To what degree is the inference of delays vs. communications an ill-posed problem? Is there a degeneracy where the true message could have static delays in terms of the transmission channel but it looks like a time-varying delay due to the time-varying nature of the communications?", "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": "cjral8x2bj", "reviewer": "Reviewer_sX24", "summary": "This paper proposes a statistical model for estimating the time-varying delay of communication between multiple brain regions. This is achieved through low-dimensional latent variable modeling and incorporating multi-output Gaussian Process models. In addition to the modeling contribution, the paper exploits a connection between state-space models (SSM) and Gaussian Process (GP) with arbitrary temporally stationary kernels which allows the authors to develop fast and scalable inference models which scale logarithmically with the number of time points. Results are shown on a toy dataset as well as two neuroscience datasets with recordings from 2 and 5 regions respectively. The model discovers interesting interplay between feedforward and feedback flow of information across regions.", "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": "Dear Reviewer sX24,\n\nThank you for the suggestion! We hope these improvements clarified your concerns, and that they can be taken into account when deciding the final score.\n\nThe additional results are: (https://anonymous.4open.science/r/rebuttal-figures-for-ICML-2025-42E6/rebuttal_figures.pdf), including Rebuttal-Figures 1&4.\n\n> Claim anything causal from the delay?\n\n- We thank the reviewer for this thoughtful comment. We agree that delay or directional estimates from observational data alone do not prove causality. We do not make causal claims here; instead, we use these estimates to form hypotheses about across-region information flow. We acknowledge that only direct manipulations can establish true causal effects. Our model’s value is in highlighting putative, time-varying interactions, which can help pinpoint the circuits and time points most relevant for further causal testing. We will elaborate on this limitation in future revisions.\n\n> More complex sensory stimulus\n\n- Thanks for this suggestion. We agree that more complex stimuli offer exciting possibilities, and we see this as an important direction for future work.\n\n> Parallel scan access to O(T) GPUs.\n\n- Parallel scan algorithm relies on having enough parallel threads rather than literally requiring $T$ separate GPUs. We will clarify this in limitation part.\n\n> Improve the diagonal covariance model to a full covariance model?\n\n- In FA model, the latents approximates a full kernel GP, so the marginal covariance $\\mathbf{C}Cov(\\boldsymbol{x})\\mathbf{C}^\\top+\\mathbf{R} $ is full.\n\n> Where are you computing the complete LL?\n\n- In the E-step, we use a Kalman filter/smoother to compute the posterior of the latents and compute the expected complete LL. In the M-step, we maximize this expectation with respect to model parameters using gradient descent, as part of the EM framework rather than directly optimizing the marginal LL. We will provide a pseudo algorithm in future revisions.\n\n> Approximation gap between SSM and GP.\n\n- Although our model is linear in state evolution, it can capture the nonlinear covariance structure of GP with a sufficient order $P$ in Eq. 5. In practice, choosing an appropriate $P$ results in performance that closely matches that of a GP.\n\n> How do you deal with Poisson data?\n\n- We preprocess using Gaussian smoothing, following prior work setting (Li et al., 2024). We acknowledge this is not ideal and developing a Poisson Kalmen Parallel scan is a key direction for future work.\n\n> Why not comparing against SLDS?\n\n- SLDS is not designed to capture temporal delays between regions. But we also show log-likelihood comparisons with multi-region SLDS in Rebuttal-Figure 4.\n\n> What enforces orthogonality of the within vs. across variables?\n\n- We do not claim strict orthogonality; instead, the FA model encourages decorrelation between variables through a block-diagonal observation matrix that assigns latents exclusively to shared or independent activity. Across-region variables use separate kernels, encouraging different delay characteristics; however, if the data indicates similar delays, the latents can reflect that similarity.\n\n> Do you also perform cross-validation on the number of latents?\n\n- Yes, we perform cross-validation on all hyperparameters, including the number of latents, in Figure 3(C). Also, in Rebuttal Figure 1(D), we test cases with larger and smaller latent numbers than true value. Incorrect latent settings result in larger variance across runs and lower test log-likelihood. We will address this issue in future revisions.\n\n> How are the shaded errorbars computed? Degrees of freedom?\n\n- The shaded error bars represent the variance of the learned delay across runs. We initialize the projection matrix $\\mathbf{C}$ using CCA, which yields stable fits across random seeds. Our model incorporates a shared kernel length scale over time, which reduces the model’s degrees of freedom by determining the temporal dynamics and constraining the evolution of delays.\n\n> Can you include state space visualization?\n\n- We visualized our model’s communication subspace on both neural datasets in Figure 2(A) and Appendix E.\n\n> Performance vs number of regions, latents, and length?\n\n- We thank the reviewer for this insightful suggestion. As shown in Rebuttal Figure 1(A-C), the model demonstrates stable performance across different conditions.\n\n> The interplay between the kernel and the delay?\n\n- The estimated delay is an effective parameter that captures the time shift between signals, as defined by the chosen kernel. While different kernels may yield varying delay values, the directional information they convey remains consistent. The sign and relative ordering of delays are stable across reasonable kernel choices, as shown in prior work (Li et al., 2024), where a different kernel produced the same directional insights on the two regions data.\n\nAgain, we thank the reviewer for insightful suggestions, which improved the quality of our paper!", "claims_and_evidence": "The claims regarding the speed of inference and the accuracy of the estimation for ground truth experiments are indeed supported by the results. For results on neuroscience datasets, the biological relevance of the estimated delays is still an open question. While the method provides valuable insights into the dynamic communication across regions, further validation experiments are required to support whether the delays have biological or causal relevance. Addressing the causality aspect of the estimated delays is beyond the scope of this paper and requires further careful investigation.", "methods_and_evaluation_criteria": "Methods that jointly model the low-dimensional latent evolution of neural dynamics along with modeling the dynamic communication across regions did not exist before. Indeed this paper addresses an important and open question in neuroscience data analysis and the modeling framework used here is sensible and appropriate. In addition, the authors discover a deeper connection between SSM and GP with temporally stationary kernels with opens up new research direction to the broader ICML community.", "theoretical_claims": "The derivations seem correct to me. Apart from those the paper does not include theoretical claims. A few approximations are proposed in the paper. Addressing the quality of those approximations can be further investigated. I will explain this below in more detail.", "experimental_designs_or_analyses": "The datasets considered and the application of the model to those datasets is sounds and appropriate.", "relation_to_broader_scientific_literature": "The connection between SSM and GP can have potential applications that go beyond the results shown in the paper. In addition, the methodology is relevant to the computational and experimental neuroscience community and enables them to gain new insight into multi-region datasets.", "essential_references_not_discussed": "Nothing off the top of my mind.", "other_strengths_and_weaknesses": "**Strengths**\n\n- The paper is very well-formatted and well-written. I had a relatively easy job following the arguments and derivations.\n- The accompanied code package follows standard code development practices.\n- The connection between SSM and GP with temporally stationary kernels is of potential interest to the broader ML community.\n- The method opens up new interesting analyses to the neuroscience community enabling new insights about multi-region communication.\n\n**Areas of improvement**\n\n- In some figures your description of the results has a causal flavor to it (feedforward vs. feedback direction of information). Imagine a scenario in which you have two regions exhibiting sine and cosine waves that are phase shifted; essentially it’s impossible to tell which one is sending/receiving information to the other, unless one region is manipulated and the effect is seen in the other region (but not in the opposite direction). Therefore it sounds like it’s impossible to claim anything causal from the estimated delay values. While this is a fundamental challenge associated with all the modeling approaches that do not consider causal manipulations, it’d be great if the authors can discuss this in the limitation section of their work. That said, I think there’s still a lot of value with the development of methods like the one proposed in the paper and it can help build insights and hypotheses that can then be causally tested.\n\n- Both datasets presented in the paper use drifting grating stimuli which induce synchronous activity across brain regions; it’d be nice if the authors could show results on datasets with more complex sensory stimulus or datasets recorded from freely behaving animals.", "other_comments_or_suggestions": "**Minor issues**\n\n- Please reword the following, there are two many communications in one sentence: [Understanding and constructing brain communications that capture dynamic communications across …].\n- Please fix [Dicovering] -> discovering.\n- Please fix [Dimmention] -> dimension.", "questions_for_authors": "- Parallel scan $O(logT)$ that’s assuming you have access to O(T) GPUs, right? The arguments in the paper give the impression that $O(logT)$ comes for free. It'd be great if the authors could discuss this in the limitations section.\n- Can you improve the diagonal covariance model and incorporate a full covariance model (or low-rank + diagonal)?\n- I didn’t quite understand the inference, where are you computing the complete log likelihood? If you're running gradient descent isn't your algorithm essentially performing marginal likelihood optimization?\n- Samples from a GP are nonlinear functions, looks like everything in your modeling framework is linear though, can you characterize the approximation gap you incur by treating a GP as an SSM?\n- How do you deal with Poisson data? Can you expand the model to account for Poisson observations? In the current settings, do you apply any transformations to your spike counts to make them consumable by the model?\n- Why not comparing against SLDS? Can you include these comparisons?\n- What enforces orthogonality of the within vs. between variables? Also among “between” variables what makes them orthogonal? Is it not possible for all of the between variables to reflect the same temporal delay characteristics?\n- You perform cross-validation on the delay, do you also perform cross-validation on the number of latents? If not, how do you set it? Can you show cross-validation results to get an overall sense of the sensitivity of the method to these parameters? Also can you include experiments where you set the number of latents to larger or smaller than the true value and see the effect? Specifically, if you set to larger this should result in some non-identifiability and inconsistency across runs or random initializations; this would specifically be an issue since you’re using the same number of between area latents which is practical assumption but might not be the the right assumption.\n- How are the shaded errorbars computed? If you run the  method multiple times with different initializations, do you get the same result (same latents, same estimated delays)? How well defined is the problem? There seems to be many degrees of freedom in the model which suggests that multiple solutions can possibly co-exist (which is perhaps exacerbated by the existence of both positive and negative delays). \n- Can you include state space visualizations too (e.g. PCA, or projections to the communication subspace)?\n- Can you expand the simulation results to study the performance of the model (in terms of the recovery of the true delay parameter) as a function of number of regions, latents, delays, trials, smoothness of the kernel, etc.?\n- What's the interplay between the choice of the kernel and the estimated delay. The estimated delay parameter seems to be dependent on the choice of kernel. In other words, if we change the kernel, the estimated delay parameter will be different too. Given this non-identifiability, how should we interpret the estimated delay values from a physical or biological standpoint?", "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 paper addresses itself to the estimation of communication between different recording sites in multi-region brain activity data. This problem has received increasing attention in the literature as more and more experiments generate these data, and as a result, several methods have been proposed, all of which necessarily trade off computational efficiency and flexibility. 

This work proposes to approximate a flexible but inefficient Gaussian Process model with between-region time delays by a multi-lag state space model (SSM). While similar papers such as Li et al. (2024) have attempted this link, this paper provides an alternative approach that allows for both a) inferring time-varying delays between regions and b) computationally efficient parallel scan operations.

Reviewers found the paper clearly written and convincing, with well-chosen benchmarks and model comparisons. While the model does make some approximations (detailed by the authors), it clearly represents a new state of the art in analyzing communication between recording sites in multi-region brain data.
Decision: Accept (oral)