Uncovering motifs of concurrent signaling across multiple neuronal populations
Authors: Evren Gokcen
Anna Ivic Jasper
Alison Xu
Adam Kohn
Christian K. Machens
Byron M. Yu
Venue: NeurIPS 2023 spotlight (2023)
Keywords: neuroscience
multi-population neural recordings
dimensionality reduction
latent variable models
Gaussian processes
TL;DR: We developed a dimensionality reduction framework for characterizing the multi-dimensional, concurrent flow of signals across multiple neuronal populations.

Abstract: Modern recording techniques now allow us to record from distinct neuronal populations in different brain networks. However, especially as we consider multiple (more than two) populations, new conceptual and statistical frameworks are needed to characterize the multi-dimensional, concurrent flow of signals among these populations. Here, we develop a dimensionality reduction framework that determines (1) the subset of populations described by each latent dimension, (2) the direction of signal flow among those populations, and (3) how those signals evolve over time within and across experimental trials. We illustrate these features in simulation, and further validate the method by applying it to previously studied recordings from neuronal populations in macaque visual areas V1 and V2. Then we study interactions across select laminar compartments of areas V1, V2, and V3d, recorded simultaneously with multiple Neuropixels probes. Our approach uncovered signatures of selective communication across these three areas that related to their retinotopic alignment. This work advances the study of concurrent signaling across multiple neuronal populations.

Uncovering motifs of concurrent signaling across multiple neuronal populations
Evren Gokcen', Anna Jasper'. Alison Xu²,
Adam Kohn 223 Christian K. Machens Byron M. Yu
Deptt. of Electrical and Computer Engineering, Carnegie Mellon University
Dominick Purpura Dept. of Neuroscience, 'Dept. of Ophthalmology and Visual Sciences, Dept. of Systems and Computational Biology, Albert Einstein College of Medicine Chhampalimaud Neuroscience Programme, Champalimaud Foundation Dept. of Biomedical Engincering, Carnegie Mellon University egokcen@emu. edu, (anna. jasper, alison. xu, adam. moonh]@eistteimeed edu, christian machens@neuro. f ccampal imaud. org, byronyu@cmu.c .dd "Denotes equal contribution.
Abstract
Modern recording techniques now allow us to record from distinct neuronal popu- lations in different brain networks. However, especially as we consider multiple (more than two) populations, new conceptual and statistical frameworks are needed characterize the multi-dimensional, concurrent flow of signals among these populations. Here, we develop dimensionality reduction framework that deter- mines (1) the subset of populations described by each latent dimension, (2) the direction of signal flow among those populations, and (3) how those signals evolve over time within and across experimental trials. illustrate these features simulation, and further validate the method by applying it previously studied recordings from neuronal populations in macaque visual areas V1 and V2. Then study interactions across select laminar compartments of areas VI, V2, and V3d, recorded simultancously with multiple Neuropixels probes. Our approach uncovered signatures of selective communication across these three areas that related to their retinotopic alignment. This work advances the study of concurrent signaling across multiple neuronal populations.
Introduction
Cortical circuits functionally involve feedforward, feedback, and horizontal interactions between many neuronal populations that span distinct areas and layers. Recording techniques now allow us to record from many neurons across these populations ( (Fig. [T).. To capitalize the scientific opportunities presented by these recordings, however, new conceptual and statistical frameworks are needed, particularly as we consider communication across multiple (more than two) populations.
Characterizing interactions between just two populations is a challenging high-dimensionall problem. Dimensionality reduction techniques have therefore been increasingly used for this purpose 74 Methods like canonical correlation analysis (CCA) [8] and its probabilistic variants (9), in particular, identify www-dimensional set of latent variables that parsimoniously describe the interactions between two populations. This cross-population shared-latent model has inspired several extensions targeted toward neural recordings 140
Communication between two populations, however, occurs bidirectionally and likely concurrentl,, and disentangling this concurrent communication is substantial challenge in neuroscience. A
37th Conference on Neural Information Processing Systems (NcurIPS 2023).

---

a
Population A
Population B
Neurons Trials
Aa
Neural activity
B2
Latent activity
Time
Activity pattern
Neuron A
Time
Neuron B,
Two distinct interactions?
One global interaction?
B
Figure 1: Challenges in studying concurrent signaling across multiple populations of neurons. (a) High-dimensional, simultaneous recordings across multiple neuronal populations. (b) Disentangling concurrent signal flow. (c) Distinguishing network-level interactions.
recently developed dimensionality reduction approach, delayed latents across groups (DLAG) [5] addresses this challenge by leveraging two insights: (1) communication between two populations is not instantaneous, and (2) while concurrently relayed signals might be difficult to tease apart from the neural activity (Fig- Dp. black trajectories in each population space), they could be pinpointed along certain dimensions or axes (Fig. [[.. "latent activity" measured along the magenta and gray dimensions; note how activity in one population leads activity in the other population).
Pairwise methods are ill-equipped, however, to analyze communication across three or more recorded populations. Suppose, for example, that wish to study the interactions of three recorded popu- lations, A, B, and C (Fig. One could consider applying one of the aforementioned methods to each pair of populations. However, Suppose that populations A and B exhibit shared activity fluctuations, and populations A and C also exhibit shared fluctuations. Do populations A. B, and C all co-fluctuate together (Fig. De. right; black arrows: direction of influence)? Or do A and B co-fluctuate in a way that is uncorrel ated with the way in which A and C co-fluctuate (Fig. [T. left)? Only by analyzing all populations together can we differentiate these possibilities. Multi- multi-popplltion dimensionality reduction approaches, such as group factor analysis (GFA) [16). could be applied toward that end, but the challenge of disentangling concurrent, bidirectional signaling remains. New approaches are needed to jointly characterize the multi-dimensional, concurrent flow of signals among multiple populations.
would encounter the following interpretational ambiguity.
We therefore propose multi multi-poppulation DLAG (mDLAG), a dimensionality reduction framework that determines (1) the subset of populations described by each latent dimension, (2) the direction of signal flow among those populations, and (3) how those signals evolve over time within and across experimental trials. We illustrate these features in simulation, and further validate mDLAG by applying it to previously studied recordings from neuronal populations in macaque visual areas V1 and V2. Then we study interactions across select laminar compartments of areas V1, V2, and V3d, recorded simultancously with multiple Neuropixels probes. mDLAG uncovered signatures of selective communication across these three areas that related their retinotopic alignment. Throughout the analyses of simulated data and neural recordings, mDLAG provided improved interpretation and quantitative performance relative alternative approaches. This work advances the study of concurrent signaling across multiple neuronal populations and its role in brain function.
2 Delayed latents across multiple groups (mDLAG)
Observation model and automatic relevance determination For population m (comprising In neurons) at time i on trial we define a linear relationship between observed activity, ym E 19m.
and latent variables, m R2 (Fig. h):
n.t + m +e"
E" ~ N(O, (1m
(2)

---

b
GP auto- covariance
GP parameters
A
C
Automatic relevance determination
Ø Tr
B
C2
Dmn
Dmis
GP cross covariance
1
dm
Ym
C
m
C
N trials
M populations
Latent
B B- 1/2
0
A
Figure 2: Delayed latents across multiple groups (mDLAG) (a) mDLAG directed graphical model representation. Filled circles represent observed variables. Unfilled circles represent probabilistic latent variables and parameters. Black dots represent point estimates. Arrows indicate conditional dependence relationships between variables. (b) Example mDLAG loading matrix (here the loading matrices for individual populations, C² C² and C have been concatenated vertically). Each element of the matrix is represented by square: magnitude is represented by the square's 's area, and sign is represented by the square's 's color (red: positive; blue: negative). Note the population-wist sparsity pattern of each latent variable, which is estimated from the neural activity. (c) Gaussian process auto-covariance (top) and cross-covariance (bottom) functions in the mDLAG model. Here chose the squared exponential covariance function. The width of the auto- and cross-covariances corresponds to timescale parameter (Ty). The center of the cross-covariance between populations mi and m2 (e.g., population A and B, respectively) corresponds to the relative time delay Dmaadd -mmidt
where the loading matrix C e Rm XP and mean parameter d E R9m model parameters E a zero-mean Gaussian random variable with noise precision matrix 4 E 89m Xqm (8 Tm is the set of 4m X qm symmetric matrices). As will describe, at time point t, latent variables x m = 1..... M are coupled across populations, and thus each population has the same number of latents, p. Because we seek low-dimensional description of neural activity, the number of latent variables is less than the number neurons, i.e., < q, where = =m 9m-
A core goal of the mDLAG framework is to identify and distinguish multiple network-level inter- actions (Fig. [.). To do so requires identifying the number latent variables across all recorded populations, and which subset of populations each latent involves, in computationally tractable man- ner. We take a Bayesian approach to this problem (thereby avoiding computationally intensive grid search, see Discussion), and let d!" dT and C be probabilistic parameters with prior distributions.
The parameter d describes the mean firing rate of each neuron. We set a Gaussian prior over d"
p(d = nd 10
where 8 E Rz0 is a hyperparameter, and I is the 9m X 4m identity matrix. We constrain the precision matrix 4th diag(oT) ..... qm be diagonal to capture variance that is independent to each neuron. This constraint encourages the latent variables to explain as much of the shared variance among neurons as possible. We set the conjugate Gamma prior over each 6.n for each neuron i 1, ......9m Qmi
P( [(O ass bz)
where ag. by E Rso are hyperparameters.
The loading matrix CT linearly combines latent variables and maps them to observed neural activity. In particular, the ,th column of C"", c² E Th9 maps the jth latent variable r to population m. To determine which subset of populations is described by each latent, we employ automatic relevance determination (ARD), which has been used successfully in a variety of contexts [7 [8] Specifically, we define the following prior over the columns of each C
P(c (.)) = =N(c," 0, (m P(a.)) l(aa"" T
(5)
da. ba)
(6)
'We will define all variables they appear, but Supplementary Section S1 for more complete notation.

---

where aj" E Ryo is the ARD parameter for latent variable and population m, and aa ba E Ry0 are hyperparameters. As oz becomes large, the magnitude of cj becomes concentrated around and hence the jth latent variable " "n.j. will have a vanishing influence on population m. The ARD prior encourages population- wwise sparsity for each latent variable during model fitting (see below), where the loading matrix coefficients will be pushed toward zero for latent variables that explain an insignificant amount of shared variance within population, and remain non-zero otherwise (Fig. (F)
The parameter C also has an intuitive geometric interpretation. Each element of y the activity of each neuron in population m, can be represented an axis in high sigh-dimensional population activity space (Fig. IT). Then each column of CTMH defines a dimension in population m's m's activity space (Fig. Dp. "activity pattern")). Dimensions that appear in two more populations represent patterns of activity that are correlated across populations. Note that the columns of C.Th are linearly independent; but they are not, in general, orthogonal. The ordering of these columns, and of the corresponding latent variables, is arbitrary.
State model For each latent variable, seek to characterize the direction of signal flow among the involved populations (determined by ARD) and how those signals evolve over time within and across trials. We therefore employ Gaussian processes (GPs) (19]. and define GP over all time points f = =....., for each latent variable j = -1,, l,...,. as follows (Fig. 20):
K1.1
KLMJ
NN
0,
: KM.13
KM,M...
Under equation 7 latents are independent and identically distributed across trials. The diagonal blocks K1,1,j .... K..M.j E S describe the autocovariance of each latent, and each T-by-T off-diagonal block describes the cross-covariance that couples two populations.
To define these matrices, we introduce additional notation. Specifically, we indicate populations with two subscripts, mi = 1,..., M and m2 = 1,..., M. Then, define K Kmi RTxT to be either the auto- cross-covariance matrix between latent variable x E R in population mi and latent variable x E R population m2 on trial n. We choose to use the squared exponential function for GP covariances (Fig. 2 Therefore, element (t1,t2) of each Kmi,ma.j can be computed
as follows 15 20
(At)²
km man kmy (2) = 10:00011000 (1- exp
+ to²
Es
(8)
2T²
At = (t2 - Dmzid)
(t1 Dmn
(9)
where the characteristic timescale, Tj E Rxo, and the GP noise variance, of E (0, 1), are deterministic model parameters to be estimated from neural activity. dAt is the kronecker delta, which is 1 for Al = 0 and O otherwise.
We also introduce two new parameters: the time delay to population m., Dmni. R. and the time delay to population m2, Dmgd E R. Notice that, when computing the auto-covariance for population m (i.e., when my = 772- mi; Fig. 2c. top), the time delay parameters D D and D Dinaj are equal, and so At (equation 9 reduces simply to the time difference (fz t1). Time delays are therefore only relevant when computing the covarianee between distinct populations m, and m2. The time delay to population my. Dm..j and the time delay to population m2, Dmz.jt by themselves have no physically meaningful interpretation. Their difference Dmg.j Dm j. however, represents a well-defined, continuous-ve valuued time delay from population m. to population m2 (Fig. 2 bottom). The sign of the relative time delay indicates the directionality of the lead-lag relationship between populations captured by latent variable j (positive: population my leads population m2; negative: population m2 leads population m)), which we interpret description of signal flow. Note that time delays need not be integer multiples of the sampling period or spike count bin width of the neural activity.
Without loss of generality, we designate population m = as the reference area, and fix the delays for population 1 at (i.e., D1.5 = for all latent variables = .......P. p). We follow the same conventions as in 15 21 and fix of to small value (10-3). Furthermore, the GP is normalized so

---

that mm,, (2) = 1 if At = 0, thereby removing model redundancy in the scaling of the latent variables and loading matrices C
Posterior inference and fitting the mDLAG model Let Y and X be collections of all observed neural activity and latent variables, respectively, across all time points and trials. Similarly, let d. Φ. C, A. and D be collections of the mean parameters, noise precisions, loading matrices, ARD parameters, and time delays, respectively. From the neural activity, we seek to estimate posterior distributions over the probabilistic model components 0 = (X, d, Φ, C, A} and point estimates of the deterministic GP parameters & = {D. {T)})=}}-
We do so by employing variational inference scheme [6 22 in which maximize the evidence lower bound (ELBO), L(Q,9) ?), with respect to the approximate posterior distribution Q(f) and the deterministic parameters S2, where
log P(Y) 2 (Q) Eollog ([)) - Eo[log Q(0)
(10)
We constrain Q(f) so that it factorizes over the elements of 0:
Q(0)= - [000000000x000000000000000000000
(11)
This factorization enables closed-form updates during optimization. The ELBO can then be itera- tively maximized via coordinate ascent of the factors (0) and the deterministic parameters f2 Here all hyperparameters were fixed to a very small value, 8, ag, ba. Ga ba = 10-12, to produce noninformative priors [16). Throughout this work, we take estimates of the latent variables and model parameters to be the corresponding means of the posterior distributions comprising equation [1 Full
details are provided in Supplementary Section 52
Finally, we consider some special cases of the mDLAG model that illustrate mDLAG special cases
its relationship to other dimensionality reduction methods. In the case of two populations (M = 2),
mDLAG is equivalent to Bayesian formulation of DLAG [15]. In the case of one population
(M = 1), and when all time delays are fixed to zero (Dm.) 0), mDLAG becomes equivalent to a
Bayesian formulation of Gaussian process factor analysis 23 By removing temporal smoothing
(i.e., in the limit as all GP noise parameters o approach 1) mDAA becomes equivalent to GFA [16]
3 Validation in simulation and on neural recordings
Simulation 1: Uncovering directed interactions across multiple populations We first sought to demonstrate the ability of mDLAG to infer the structure of multi-population interactions. We therefore generated simulated neural activity from three populations (M 3) according to the mDLAG generative model (equations D For the sake of illustration, we set 10 neurons in each population (9m. 10). Importantly, we designed the loading matrices CTM so that all types of inter- population interactions were represented (Fig. 3R. left, ground truth): interactions shared globally. unique to each pair of populations, and local to one population. We also scaled the observation noise precision matrices d so that noise levels were representative of realistic neural activity. Specifically, activity due to single-neuron observation noise was 10 times stronger than activity due to latents (the "signal-to "signal-to-noise noise ratio" tr(C" CT ) )/tr( (r(($")=1) = 0.1 for each population). We selected Gaussian process timescales and (relative) time delays that ranged between 20 ms to 150 ms and between 5 ms and 50 ms, respectively. With all model parameters specified, we then generated N 100 independent and identically distributed trials. Each trial was 500 ms in length, comprising T = 25 time points with a sampling period of 20 ms, to mimick the 20 ms spike count time bins used to analyze the ncural recordings in Section 4
We then fit an mDLAG model to the simulated neural activity. We set the initial number of latents (p 10) to be greater than the ground truth number (p = 7), to verify that these additional latents would be pruned during fitting. In addition, to demonstrate the benefit of ARD in the mDLAG model, we also fit a modified mDLAG model that did not use ARD. Specifically, we fit (via an exact EM algorithm) a modified model with state model defined by equations 9 and observation model defined by equations] | and 2
but with no prior distributions over the parameters C dm and tm (i.e., we obtained point estimates). The number of latents for this modified model was selected via 4-fold vrossvallidation, considering model candidates with p = to p = 10 latents. Cross-validan valldaval is not needed for the mDLAG model with ARD, key computational advantage of the approach.
5

---

a
Ground truth
mDLAG (with ARD)
mDLAG ARD)
8.0 14.0 13.6
14.0
34.0
46.0
C
15.4
34.4
42.9
2.4
Latent
3
4
W
-2.4 Time
B
250 ms
A
2.4
W
Ground truth Estimate
Latent
-2.4
Figure 3: Automatic relevance determination encourages population-wise sparsity. (a) Loading matrix estimates. Left: Ground truth loading matrix. Center: mDLAG estimate. Right: mDLAG estimate where automatic relevance determination (ARD) was not used. i.e., no population-wise sparsity priors were placed on the loading matrix. Same conventions in Fig. 2p. Note that the sign and ordering of each loading matrix column is, in general, arbitrary. We have therefore reordered and flipped the signs of the columns of the estimates to facilitate comparison with the ground truth. (b) Single-trial latent time course estimates. Each panel corresponds to the ground truth and estimated time course of a single latent variable. For concision, we show only latents corresponding population (x..,.): the remaining latents are time-shifted versions of those shown here. Inset above each latent are the involved populations along with the signal flow and magnitude of time delays between populations (in ms). Magenta: mDLAG estimates; black: ground truth. a.u.: arbitrary units.
one
We were able to identify the total number of latent variables (p = 7) with both approaches. The mDLAG model with ARD additionally recovered the population-wiss sparsity structure with high accuracy (Fig-Bi Pp. center). The mDLAG model without ARD, however, produced an estimate of the loading matrix with mostly non-zero elements (Fig- [3,, right): had not known the ground truth in advance, it would be difficult to interpret which population subsets are involved in which interactions. The mDLAG model with ARD also estimated the latent variable time courses and time delays (along with their implied signal flow) with high accuracy (Fig. 3b; R22 between ground truth and estimated time courses: 0.936; mean delay error: 1.14 ms).
Simulation 2: Disentangling concurrent, bidirectional signaling We next considered simulated scenario that illustrates mDLAG's ability to distinguish concurrent, bidirectional signaling across populations, particularly in contrast with static multi-population methods like group factor analysis (GFA) We again simulated activity from three populations, each with 4m = 10 neurons. These populations interacted through two latent processes,
that propagated from population A to B to (Fig. 4h. Latent 1) and one that propagated in the opposite direction, from population to B to A (Fig. 4a. Latent 2). Each latent incurred a 15 ms time delay propagating from one population to the next. For the sake of visual demonstration, we generated latent signals with a simple transient peak of activity (Fig. 4h. black traces), and generated the same latent signals each trial. We then generated observed neural activity according to the mDLAG observation model (equations I and 2). Importantly, we scaled the columns of the loading matrices C (for each area m) that both columns had equal magnitude, and thus both latent interactions exhibited equal strength in each population: such communication scheme difficult disentangle (15]. Finally, isolate the consequences of model structure from issues like overfitting, we simulated low-noise, data-rich scenario by generating N 000 independent trials, each 500 ms length (comprising = 25 time points), with high signal-to-note noise ratios cr( t/tr( ($m) -1) = 10.0 for each population.
We fit both an mDLAG model and GFA model (see Supplementary Section S3) to this simulated neural activity, assuming for each model the correct number latent variables, = 2. mDLAG latent variable and time delay estimates accurately reflected the distinct signaling pathways across
"See also Supplementary Fig. S2 for performance comparisons between GFA and mDLAG on the neural recordings analyzed Section

---

a
mDLAG estimate
GFA estimate
Latent
Latent 2
14.7
15.0
14.8
14.9 B
B
-
C
A
Time 15 ms
15 ms
2.7
Latent 1
2
2.2
Latent 1
2
One trial
M
0
au.
0
15 ms
15 ms
-2.7
100 ms
-2.2
Time
Figure 4: mDLAG disentangles concurrent signaling where static methods do not. (a) Schematic of simulated ground truth. Populations A, B, and interact through two latents. Latent (left) propagates from A to B to C. Latent 2 (right) propagates from C to B to A. (b) mDLAG estimates disentangle the two interactions. Each panel corresponds to a single latent variable. Each black trace represents one trial (10 representative trials shown). For concision, we show only latents corresponding to one population (x..j..); the remaining latents are time-shifted versions of those shown here. Inset above each latent the involved populations along with the signal flow and magnitude of time delays between populations (in ms). a.u.: arbitrary units. (c) Group factor analysis (GFA) estimates represent mixtures of the two interactions. Same conventions as in (b).
the three populations (Fig. 40; R² between ground truth and estimated time courses: 0.989; mean delay error: 1.16 ms). Each latent estimated by GFA, however, notably reflected a mixture of both interactions (Fig. 4c. each latent time course exhibits two peaks). With no description of signal flow, static methods like GFA struggle to distinguish concurrent, bidirectional signaling.
Validating mDLAG on recordings from V1 and V2 To further validate mDLAG, we turned to neural recordings from two in the macaque visual cortex, V1 and V2 (24]. These recordings have been studied extensively [15] 27), and therefore provide an excellent testbed for our framework. We applied mDLAG the simultaneously recorded activity of neuronal populations in the superficial (output) layers of V1 (61 to 122 neurons; mean 86.3), and the middle (input) layers of V2 (15 32 neurons; mean 19.6) in three anesthetized animals (Supplementary Fig. (TT).. We analyzed neuronal responses measured during the 1.28 second presentation of drifting sinusoidal gratings of different orientations, and counted spikes 20 ms time bins. Because we were interested in V1-V2 interactions timescales within a trial, we subtracted the mean across time bins within each trial from each neuron. This step removed activity that fluctuated on slow timescales from one stimulus presentation to the next [28). In total, we analyzed separately 40 "datasets," corresponding to five recording sessions, each with eight different orientations. Each dataset included 400 trials: we allocated at random 300 trials performance (we employed a cave-group-out prediction metric; see Supplementary Section SS).
training set and 100 trials a test set which to measure model
Given prior work on these V1-V2 recordings, we reasoned that mDLAG ought to qualitatively recover several hallmarks of the V1-V2 activity: (1) a significant number of latent variables local to each area, indicative of selective V1-V2 communication (Supplementary Fig. [SIp. 'V1' and 'V2'; Supplementary Fig. ST latents 4 and 5), (2) the presence of periodic structure representative the drifting grating stimulus (Supplementary Fig. STc latents and 4), and (3) bidirectional interactions (Supplementary Fig. STC. latents 1-3). Indeed, mDLAG detected all of these hallmarks.
These recordings were also previously studied with DLAG [5], which was designed to study concurrent signaling between two neuronal populations. Since DLAG is special case of the mDLAG framework, we reasoned that mDLAG ought to exhibit similar quantitative performance and extract qualitatively similar features of V1-V2 communication. Indeed, the latent time courses estimated by mDLAG (Supplementary Fig. ( were qualitatively similar to those estimated by DLAG [29] (Supplementary (Td. Interestingly, mDLAG outperformed DLAG across all datasets (Supplementary Fig. S22, points above the diagonal), suggesting that ARD provides an improved method of model selection over the constrained grid search method used for DLAG, while also avoiding grid search's computational drawbacks (see Supplementary Section SS).. Collectively, the results across our simulations and the VI-V2 recordings indicate that mDLAG well-equipped to study concurrent signaling across multiple simultancously recorded neuronal populations.
7

---

Neuropixels
Neuron
probes
V1
2
V1
V2
V3d
log tring rate
ayer
M
M
D
V3d
V3d S/M
562
1,000
Time from stim. onset (ms)
d
2 ms ms
ms
10 ms
V2
V1
No. shared dims.
Unique
v3d
V V1
V1 V3d
V1 V3d
4.4
Latent
3
Global
0
2-
One trial
V3d
V1-V2 V1-Vad
V2-V3d
4.4 Time
333 ms
2
Figure 5: Applying mDLAG to Neuropixels recordings across VI, V2, and V3d. (a) Schematic showing sagittal section of occipital cortex and the recording setup. Laminar compartments of VI. V2, and V3d were recorded using multiple Neuropixels probes. (b) Temporally smoothed peristimulus time histograms during the stimulus presentation period, for an example session and stimulus condition. (c) Schematic showing the subset of laminar compartments included in the analyses (shaded compartments). Compartments with the same color were treated as one population during analysis. Arrows indicate prominent anatomical connections [30]. (d) Visual receptive field locations of V1, V2, and V3d. Each circle represents receptive field size and location across the channels of single Neuropixels probe. Scale hars indicate degrees of visual field. (e) Total number of latent dimensions shared between each pair of populations. Light gray: dimensions unique to pair. Dark gray: dimensions shared globally across all populations.
(f) Single-trial time course estimates of select shared latents (same dataset as in (b), (d), and (e)). Same conventions as in Fig.
4 Interactions across laminar compartments of V1, V2, and V3d
We next used mDLAG to study interactions across VI, V2, and V3d of anesthetized macaques recorded using multiple Neuropixels probes (Fig. 5p). These recordings included hundreds of neurons across different layers in each area (Fig-5p. superficial (S), middle (M), and deep (D) compartments) allowing us investigate longstanding questions about functional cortical organization.
(Fi
First, we sought to quantify the performance of mDLAG models relative to alternative approaches, and to further demonstrate the empirical benefit of mDLAG's ' model components. We therefore considered three types of models: group factor analysis (GFA); mDLAG, but with time delays removed ("mDLAG-0") ); and mDLAG. With each model type, we analyzed spike counts (in 20 ms bins) of neurons in select laminar compartments of each area, measured during the first 1,000 ms after the onset of drifting sinusoidal grating (Fig. (p,, trial-averaged responses for an example dataset). The selected laminar compartments (grouped into three sub-populations for analysis: Fig. 5). color shading) were those that likely directly interacted through combination of feedforward, feedback, and inter-l connections (30] (Fig. 55. black arrows). As we did for the V1-V2 recordings, above, we subtracted the mean across time bins within each trial from each neuron, to remove slow fluctuations beyond the timescale of a trial. We fit GFA, mDLAG-0, and mDLAG separately to 10 datasets, comprising five recording sessions (each with two grating orientations, 90" apart, presented 300 trials each) from two anesthetized animals. For each dataset, we allocated at random 225 trials as a training set and 75 trials as test set on which measure model performance. GFA, mDLAG-0, and mDLAG exhibited increasingly better leave- leave-group-on group-out prediction (Supplementary Fig. S2p: mDLAG- better than GFA; c: mDLAG better than mDLAG-0), demonstrating the performance benefit of (1) temporal smoothing and (2) time delays (see also Supplementary Fig a demonstration of mDLAG performance versus number of available training trials, and Supplementary Fig. S4 for characterization of GFA and mDAAG runtimes).
8

---

We then focused on one particular set of questions: Do VI, V2, and V3d communicate selectively, and if so, what aspects of the recorded populations might contribute this selectivity? One likely aspect is the alignment of visual receptive fields (RFs): retinotopically aligned neuronal populations might co-contribute to computations that non-aligned populations would not. This functional organization would appear as mDLAG latent variables shared exclusively between the aligned populations. We therefore studied a dataset from a recording session in which the retinotopy of the recorded populations allowed us to investigate this hypothesis (Fig. (t). RFs of the VI and V3d populations were largely overlapping (Fig. 5i. blue: VI, gold: V3d), whereas the RFs in V2 did not overlap with either VI or V3d (Fig. 5d, red: V2).
We quantified the total number of dimensions shared between each pair of populations (Fig. 5p: a dimension was considered significant in a population if it explained at least 2% of the shared vanance, Supplementary Section SS5] equation (552). mDLAA allowed us to distinguish dimensions that were unique to pair of populations (e.g., between VI and V3d, but not V2) from dimensions that also involved the third population. V1-V3d interactions included unique dimensions (Fig. 5e, 'V1-V3d' light gray) in addition to a dimension shared globally across all three populations (Fig. 55. 'VI-V3d', dark gray), whereas V1-V2 interactions and V2-V3d interactions could be attributed primarily to the global dimension (Fig. St. 'V1-V2' and 'V2-V3d', only dark gray). The single-trial time courses of the global latent evolved according to a long timescale (Fig. latent 1, 132 ms) relative to the unique V1-V3d latents (Fig. 5T. latents 2 and 3, 8 ms and 62 ms, respectively). The global latent was also associated with a direction of signal flow (V3d to V2 to V1) opposite to that of the VI-V3d latents (V1 to V3d). The time courses of the V V1-V3d latents exhibited features of the stimulus response, including a fast transient response (Fig. 5 latent 2) and periodic structure with the same period as the drifting grating stimulus (Fig. 5T. latent 3, 333 period). Results were largely insensitive the initialization of the mDLAG fitting procedure (Supplementary Fig. (S5). These findings are consistent with the hypothesis that visual cortical populations communicate selectively in retinotopic manner, as they perform computations on inputs from localized regions of visual space [31].
Discussion
Even for two populations, identifying which latents are shared between populations or local to population is a challenging computational problem. Existing approaches explicitly designate latents as shared local 15). and then rely on heuristics to avoid the computational cost of grid search. Scaling this type of approach to three or more populations (or external experimental variables) would be prohibitive. mDLAG is thus advance toward scaling to large-scale multi-populatee recordings, and could continue to be improved through the many approaches used to scale Gaussian process a
methods in neuroscience 13 23 32 36
mDLAG treats time delays as constant parameters across trials, time, and neurons. The delay estimated for each latent variable thus represents summary of inter nntur-population interaction throughout the course of an experiment. To assess the variability of delay estimates, one could fit mDLAG to subsets of trials, subsets of neurons, or to separate trial epochs. If additional subpopulation labels are available (for example, laminar compartments in our Neuropixels recordings, Fig. 5D), one could also specify these subpopulations to the mDLAG model to better incorporate variability of delays across neurons.
A latent variable shared across three or more populations is potentially consistent with multiple signaling schemes. For instance, Latent of Fig. 5 is consistent with signal relayed from V3d to V2 (with delay of 2 ms) and then from V2 to V1 (with a delay of 4 ms). It is also consistent, however, with a scheme in which V3d is a common input to V2 (with a delay of ms) and to V1 (with delay of 6 ms). A third scheme could involve common input to all three areas from unobserved source. Still, mDLAG narrows the set of populations that could be involved in any given interaction, and the sign and magnitude of mDLAG's time delays narrow the set of signaling schemes consistent with the data. This hypothesis set can be narrowed further by experimental interventions [37].
The mDLAG model includes assumptions of linearity and temporal smoothness; specifically, we have employed here the commonly used squared exponential function for mDLAG's GP covariances. It might be desirable to incorporate temporal structure more appropriate for certain signals (for example, the sinusoidal structure of Latent 3, Fig. 55)) Then, an alternative GP covariance function can be substituted in straightforward manner into the mDLAG model (equation g
[99
32
33
38
9

---

Nonlinearities, for example related to neuronal tuning or receptive fields, could also be incorporated
under GP framework (39 40).
For exploratory data analysis, mDLAG's GP-hased description of multi-population temporal structure is advantageous over alternative linear dynamical system (LDS)-based description [TO in two respects: (1) a GP can be useful for exploratory data analyses where an appropriate parametric dynamical model is unknown a priori, and (2) mDLAG's continuous-time model enables the discovery of wide-ranging delays with high precision, which, in contrast to discrete-time LDS approaches, are not limited to be integer multiples of the sampling period or spike count bin width of the neural activity. Ultimately, these approaches can be complementary: one can use mDLAG to generate data-driven hypotheses about motifs of concurrent signaling across populations, and then test these hypotheses with a dynamical system-based approach.
Finally, while we applied mDLAG to the spiking activity of neuronal populations, mDLAG is applicable to any nigh-dimensional time series data, including other neural recording modalities. In fact, the groups analyzed by mDLAG need not all be neuronal populations, but could include, for example, a collection of dynamic stimulus or behavioral variables. mDLAG is general framework that advances the study of concurrent signaling throughout the brain.
Code availability
A MATLAB (MakhWorks) implementation of mDLAG is available on GitHub at http:///ittull com/eggoccen/mLAAG and on Zenodo at hhttps://dii.org:: 10t.///eood///doodoooooo 10048163 [4]
Acknowledgments
We thank A. Zandvakili for providing the V1-V2 recordings in Supplementary Fig. S1 This work was supported by the Dowd Fellowship (E. ECE Benjamin Garver Lamme/Westine Weeteighouueeeu Fellowship (E.G.), Simons Collaboration the Global Brain 542999 (A.K.). 543009 (C.K.M.), 543065 (B.M.Y.), NIH RFI NS127107 (A.K., C.K.M., B.M. Y.). NIH ROI EY035896 (B.M.Y. Y., A.K., C.K.M.), NIH UO1 NS094288 (C.K.M.), NIH CRCNS R01 MH1 18929 (B.M.Y.), NSF NCS DRL 2124066 (B.M.Y.,, NIH R01 N8129584 (B.M.Y.).
References
[1] Ahrens, M. B. et al. Brain-wide neuronal dynamics during motor adaptation in zebrafish. Nature 485, 471-477 (2012).
[2] Yang, W. & Yuste, R. In vivo imaging of neural activity. Nature Methods 14, 349-359 (2017).
[3] Steinmetz, N. A. et al. Neuropixels 2.0: A miniaturized high-density probe for stable, long-term brain recordings. Science 372, eabf4588 (2021).
[4] Semedo, J. D., Gokcen, E., Machens, C. K., Kohn, A. & Yu, B. M. Statistical methods for dissecting interactions between brain areas. Current Opinion in Neurobiology 65, 59-69 (2020).
[5] Kang, B. Druckmann, S. Approaches to inferring multi-regional interactions from simulta- neous population recordings. Current Opinion in Neurobiology 65, 108-119 (2020).
[6] Keeley, S. L.. Zoltowski, D. M., Aoi, M. C. & Pillow, J. W. Modeling statistical dependencies in multi-region spike train data. Current Opinion in Neurobiology 65, 194-202 (2020). [7] Machado, T. A., Kauvar, I. V. & Deisseroth, K. Multiregion neuronal activity: the forest and the trees. Nature Reviews Neuroscience 23, 683-704 (2022).
[8] Zhuang, X., Yang. 7.. & Cordes, D. A technical review of canonical correlation analysis for neuroscience applications. Human Brain Mapping 41, 3807-3833 (2020).
[9] Bach, F. R. & Jordan, M. I. A Probabilistic Interpretation of Canonical Correlation Analysis. Technical Report No. 688, Department of Statistics, University of California, Berkeley (2005).
[10] Semedo, J., Zandvakili, A., Kohn, A., Machens, C. K. & Yu, B. M. Extracting Latent Structure From Multiple Interacting Neural Populations. Advances Neural Information Processing Systems 27, 2942-292229 2950 (2014).
10

---

[11] Glaser, J., Whiteway, M., Cunningham, J. P., Paninski, L. & Linderman, S. Recurrent Switching Dynamical Systems Models for Multiple Interacting Neural Populations. Advances Neural Information Processing Systems 33, 14867-14878 (2020).
[12] Bong, H. al. Latent Dynamic Factor Analysis of nigh-Dimenniioal Neural Recordings. Advances in Neural Information Processing Systems 33, 16446-16456 (2020).
[13] Keeley. S., Aoi, M., Yu, Y.. Smith, S. & Pillow, J. W. Identifying signal and noise struc- ture in neural population activity with Gaussian process factor models. Advances in Neural Information Processing Systems 33, 13795- (2020).
[14] Balzani, E., Noel, J.-P. G., Herrero-Videe P., Angelaki, D. E. & Savin, C. A probabilistic framework for task aligned intra- and inter-area neural manifold estimation. The Eleventh International Conference on Learning Representations (2023).
[15] Gokcen, E. et al. Disentangling the flow of signals between populations of neurons. Nature Computational Science 2, 512-525 (2022).
[16] Klami, A., Virtanen, S.. Leppäaho, E. & Kaski, S. Group Factor Analysis. IEEE Transactions on Neural Networks and Learning Systems 26, 2136-2147 (2015).
[17] MacKay, D. J. Bayesian Non-Linear Modeling for the Energy Prediction Competition. ASHRAE transactions 100, 1053-1062 (1994).
[18] Neal, R. M. Bayesian Learning for Neural Networks, PhD thesis (University of Toronto, 1995).
[19] Rasmussen, C. E. Williams, K. Gaussian processes for machine learning (MIT Press, Cambridge, MA, 2006).
[20] Lakshmanan, K. C., Sadtler, P. T., Tyler-Kabara, E. C., Batista, A. & Yu, B. M. Extracting
Low-Dimensional Latent Strueture from Time Series in the Presence of Delays. Neural Computation 27, 1855-1856 (2015).
[21] Yu, B. M. et al. Gaussian- Process Factor Analysis for Now-Dimensional Single- Trial Analysis of Neural Population Activity. Journal of Neurophysiology 102, 614-635 (2009).
[22] Bishop, C. Variational principal components. 9th International Conference Artificial Neural
Networks, 509-514 (1999).
[23] Jensen, K., Kao, T.-C., Stone, J. & Hennequin, G. Scalable Bayesian GPFA with automatic relevance determination and discrete noise models. Advances Neural Information Processing
Systems 34, 1061 1-613-10626 (2021).
[24] Zandvakili, A. & Kohn, A. Paired V1-V2 neuronal spiking responses in anesthetized macaque monkey. https:///oi.ooo wtt.oo..g/....... 6080KK0277SHN (2019).
[25]
Zandvakili, A. & Kohn, A. Coordinated Neuronal Activity Enhances Corticocortical Commu- nication. Neuron 87, 827-839 (2015).
[26] Semedo, J. D., Zandvakili, A., Machens, C. K., Yu, M. & Kohn, A. Cortical Areas Interact through Communication Subspace. Neuron 102, 249-255 (2019).
[27] Semedo, J. D. et al. Feedforward and feedback interactions between visual cortical areas use
diffcrent population activity patterns. Nature Communications 13, 1099 (2022).
[28] Cowley, B. R. et al. Slow Drift of Neural Activity as a Signature of Impulsivity in Macaque Visual and Prefrontal Cortex. Neuron 108, 551-567 (2020).
[29] Gokcen, E. DLAG (v1.0.0). Zenodo. https://doiico ort./1..2// 22281//enodo.666
(2022).
[30] Salin, P. A. & Bullier, J. Corticocortical connections in the visual system: structure and function.
Physiological Reviews 75, 107-154 (1995).
[31] Carandini, M. al. Do We Know What the Early Visual System Does? Journal of Neuroscience
25, 10577-10597 (2005).
[32] Ulrich, K. R., Carlson, D. E., Dzirasa, K. & Carin, L. GP Kernels for Cross-Spectrum Analysis. Advances in Neural Information Processing Systems 28, 1999-2007 (2015).
[33] Gallagher, N. et al. Cross-Spectral Factor Analysis.Advancee in Neural Information Processing
Systems 30, 0842-6852.06 (2017).
[34] Zhao, Y. & Park, M. Variational Latent Gaussian Process for Recovering Single- -rrill Dynamics from Population Spike Trains. Neural Computation 29, 1293-1316 (2017).
[35] Duncker, L. & Sahani, M. Temporal alignment and latent Gaussian process factor inference in population spike trains. Advances in Neural Information Processing Systems 31, 10445- (2018).
11

---

[36] Dowling, M., Zhao, Y. & Park, I. M. Linear Time GPs for Inferring Latent Trajectories from Neural Spike Trains, Proceedings of the 40th International Conference on Machine Learning (2023).
[37] Reid, A. T. al. Advancing functional connectivity research from association causation. Nature Neuroscience 22, 1751-1760 (2019).
[38] Wilson, A. & Adams, R. Gaussian Process Kernels for Pattern Discovery and Extrapolation. Proceedings of the 30th International Conference on Machine Learning. 1067-1075 (2013).
[39] Wu, A., Roy, N. A., Keeley, S. Pillow, J. Gaussian process based nonlinear latent structure
discovery in multivariate spike train data, Advances in Neural Information Processing Systems
30, 3496-3505 (2017).
[40] Duncker, L., Ruda, K. M., Field, G. D. Pillow, J. W. Scalable Variational Inference for Low-Rank Spatiotemporal Receptive Fields. Neural Computation 35, 995-1027 (2023).
[41] Gokcen, E. mDLAG (v1.0.0). Zenodo. http:://oo..agg///..222/////ooo0
(2023).
12


Reviews:
[{"review_id": "OHU0eSwLC5", "reviewer": "Reviewer_DK8X", "summary": "This paper extends the delayed latents across groups (DLAG) to a recurrent general form (mDLAG), which allows analyzing the contribution of each latent dimension for multiple observational groups. Besides, the newly proposed mDLAG is able to identify the complicated directions of the information flow between groups with the corresponding time delay estimations. Experiments on two synthetic and two empirical datasets show the effectiveness of the mDLAG in analyzing multiple-grouped neural data.", "questions": "* It is hard for me to understand Fig. 1b.\n* Eq. 3, theoretically mean firing rate $\\boldsymbol d^m$ cannot have a Gaussian prior with zero means. A typical way is to model the log or logit of the firing rate since firing rates > 0.\n* Is that $\\boldsymbol x_{n,t}^m$ should be $\\boldsymbol x_{n,t}$, since from Fig. 3b, it seems like you have a shared 7-d latent across all three neural populations A, B, and C?\n* Is the inference very time-consuming?\n* Is there any other models that can solve the same task with some tiny modification, so that the performance can be compared?", "limitations": "/", "rating": 7, "confidence": 4, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "* The model derivation is logically clear. The authors provide a very detailed clarification of the notation, parameter set, and generative procedure of this generative model\n* The experiments are exhaustive, including two synthetic and two real-world datasets with very comprehensive analyses.\n* The model is intuitive and in a general form. Its special cases can reduce the model to some famous generative model like the Gaussian process factor analysis (GPFA).", "weaknesses": "* Some math notation is a bit hard to understand, especially the $\\boldsymbol x$. From my understanding, it is a high-dimensional tensor, but the author does not define this tensor in a clear way, e.g., the order of the dimensions.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": ">**Weaknesses**\n\nWe realize we could be more explicit in defining our notation for the latents, $\\mathbf{x}$. We have added the following text throughout Supplementary Section S1:\n- \"we define latent variable $j$ (out of $p$) in population $m$ at time $t$ on trial $n$ as $x^m_{n,j,t} \\in \\mathbb{R}$.\"\n- \"we represent the collection of all $p$ latent variables in population $m$ at time $t$ on trial $n$ as the vector $\\mathbf{x}^m_{n,:,t} \\in \\mathbb{R}^p$.\"\n- \"the latent variables in population $m$ on trial $n$ can be grouped into the matrix $X^m_n = [\\mathbf{x}^m_{n,:,1} \\cdots \\mathbf{x}^m_{n,:,T}] \\in \\mathbb{R}^{p \\times T}$. We represent a row of $X^m_n$ (i.e., the values of a single latent variable $j$ at all time points on trial $n$) as $\\mathbf{x}^m_{n,j,:} \\in \\mathbb{R}^T$. Finally, we can form the three-dimensional array $X^m$ by concatenating the matrices $X^m_1,\\ldots,X^m_N$ across trials along a third dimension.\"\n\n>**Questions**\n\n>- It is hard for me to understand Fig. 1b.\n\nFig. 1b illustrates the geometric intuition behind mDLAG's ability to disentangle concurrent signal flow between a pair of populations. This task is difficult given measurements of raw neural activity, i.e., the activity along axes $A_1$ and $A_2$ in population A's neural state space (Fig. 1b, left), and the activity along axes $B_1$ and $B_2$ in population B's neural state space (Fig. 1b, right). However, if it were possible to measure activity along alternative dimensions in each population's state space (Fig. 1b, magenta and gray dimensions), then it could be possible to tease apart signals concurrently relayed in opposite directions. For example, if one were to measure activity along the magenta dimension in both population A and B (Fig. 1b, middle, top black time courses), then it is apparent that the activity in population A (top trace) leads the activity in population B (bottom trace). Similarly, if one were to measure activity along the gray dimension in both population A and B (Fig. 1b, middle, bottom black time courses), then it is apparent that the activity in population B (bottom trace) leads the activity in population A (top trace).\n\nTo improve clarity, we have added the following to Fig. 1b: (1) labels for \"Population A\" (left) and \"Population B\" (right) to indicate which population corresponds to which set of axes, (2) explicit labels to indicate that axes $A_1$, $A_2$, $B_1$, and $B_2$ correspond to individual neurons, and (3) arrowheads on the gray and magenta dimensions to indicate the positive direction of the corresponding latent activity.\n\nIf the reviewer has any other suggestions to help clarify the panel, please let us know.\n\n>- Eq. 3...\n\nWe thank the reviewer for this clarification question. Please see the General Response to Reviewers, \"Positivity constraints on the mean parameter.\"\n\n>- Is that $\\mathbf{x}^m_{n,t}$...\n\nThanks to the reviewer's comment, we realized that we could have been clearer in indicating how we chose to display the latent variables. $\\mathbf{x}^m_{n,t}$ is still correct. For Fig. 3b, for example, the complete output would be 21 latent time courses (7 latents, time-delayed for each of the 3 populations). For concision, in Fig. 3b and throughout the paper, we display the set of latents for one of the populations. We have added the following text to the legends of Figs. 3 and 4:\n- \"For concision, we show only latents corresponding to one population ($\\mathbf{x}^m_{n,j,:}$); the remaining latents are time-shifted versions of those shown here.\"\n\n>- Is the inference...\n\nAcross all datasets, the average runtime per iteration was 2.5 seconds (Supplementary Fig. S3). We have also added the following text to Supplementary Fig. S3: \"We conservatively ran each mDLAG model for 50,000 iterations, resulting in an average total runtime of 34 hours. Had we used a less conservative, but still reasonable, stopping tolerance of $10^{-6}$, the average number of iterations required for convergence would have been 17,000, for an average total runtime of 11.5 hours.\" For further discussion, please see the General Response to Reviewers, \"Computational demands.\"\n\n>- Is there any other models...\n\nWe include the following qualitative and quantitative comparisons to alternative models throughout the text:\n- Section 3, Simulation 1, Fig. 3a: We compare an mDLAG model with automatic relevance determination (ARD) to an mDLAG model without ARD. We conclude that \"The mDLAG model with ARD … recovered the population-wise sparsity structure with high accuracy (Fig. 3a, center). The mDLAG model without ARD, however, produced an estimate of the loading matrix with mostly non-zero elements (Fig. 3a, right).\"\n- Section 3, Simulation 2, Fig. 4: We compare mDLAG to group factor analysis (GFA). We conclude that \"mDLAG latent variable and time delay estimates accurately reflected the distinct signaling pathways across the three populations (Fig. 4b; $R^2$ between ground truth and estimated time courses: 0.989; mean delay error: 0.16 ms). Each latent estimated by GFA, however, notably reflected a mixture of both interactions (Fig. 4c, each latent time course exhibits two peaks).\"\n- Section 3, V1-V2 recordings, Supplementary Figs. S1, S2a: We compare mDLAG to DLAG. We conclude that \"mDLAG outperformed DLAG across all datasets (Supplementary Fig. S2a, points above the diagonal)...\"\n- Section 4, Supplementary Figs. S2bc, S3: We compare mDLAG to GFA and a lesioned version of mDLAG in which time delays were fixed at zero ('mDLAG-0'). We conclude that \"GFA, mDLAG-0, and mDLAG exhibited increasingly better leave-group-out prediction (Supplementary Fig. S2b: mDLAG-0 better than GFA; c: mDLAG better than mDLAG-0)...\"\n- See also Reviewer Figs. R2 and R3, in which we show that mDLAG outperforms multiset CCA.\n\nTo our knowledge, these are the most comparable methods to mDLAG. If the reviewer is aware of other related methods that we should cite or compare to, please let us know."}, {"review_id": "M35I0pZs8z", "reviewer": "Reviewer_DL1P", "summary": "The recent developments in neural recording technologies allow recording from large populations of neurons from multiple brain regions simultaneously. Latent space models are often used to analyze these datasets, but they are generally limited to the study of single or two populations of neurons. This work expands on existing dimensionality reduction methods and introduces a new probabilistic method to characterize interactions across multiple (more than two) populations of neurons. The proposed model can capture local and shared variability across multiple populations as well as the direction of the signal flow between areas and their temporal evolution with trial resolution. The authors validated the method on simulated data and neural data, showing that the model outperforms an existing method without temporal smoothing, GFA, and lesion versions of the model.", "questions": "The simulated data has predetermined parameters, some set to match characteristics of the neural recordings and some others not. For example, the number of neurons in the recordings is in the order of dozens, while in one of the simulations, it is set to 10. What prompted the specific parameter choices? Or more interestingly, how robust is the model to potential changes in the different parameters: number of neurons, noise parameters, temporal delays, and so on? Another additional application listed for future work is the analysis of other signals like behavior. It would be relevant to know in these cases if the model requires similar dynamics or temporal binning for it to work.", "limitations": "The authors address the limitations of linearity and temporal smoothness but fail to address the potential limitations in computational cost and data demands.", "rating": 7, "confidence": 4, "soundness": 3, "presentation": 4, "contribution": 3, "strengths": "The paper is clearly presented and technically sound. The method was tested and shown to work well in simulation and neural data. This method builds on existing latent space models to enable multi-area signal characterization. Understanding inter-area population activity is an interesting problem in neuroscience, and this work proposes a new method that can help explain these emerging datasets. Importantly, even when tested in cases with only two neural populations, for which there are multiple existing methods, the proposed model still outperforms alternatives. This suggests that this tool could be broadly applied to single, two, or more than two areas population recordings. Additionally, when applied to neural data, they demonstrated the power of the approach by reporting new discoveries.", "weaknesses": "While the authors show the promise of their method to understand multi-area signals, they overlooked other existing methods that can also capture this information, such as multiset CCA or extensions of Procrustes alignment to multiple datasets. It would be interesting to add a performance comparison to these models. The model assumes linearity and Gaussian observations, but neural activity is often captured in spiking activity, which is better captured by a Poisson model. The authors show no indication of this limitation. Lastly, the computational cost of using a GP prior for temporal smoothing is relatively high compared to alternative methods, which could limit the potential applications. The method seems to be a direct extension of DLAG to multiple areas. If so, an explicit comparison between them could help assess the significance of the work.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "> **Weaknesses**\n\n> While the authors show the promise of their method to understand multi-area signals, they overlooked other existing methods...\n\nWe thank the reviewer for pointing out these alternative methods. We agree that they are relevant, and we cite several review papers that mention such methods (Semedo et al., 2020; Kang et al., 2020; Keeley et al., 2020; Zhuang et al., 2020; Machado et al., 2022). In summary, these methods fall short in addressing two key challenges in the study of multi-population recordings (Fig. 1): (1) distinguishing network-level interactions, and (2) disentangling concurrent signal flow. Below and in Reviewer Figs. R2 and R3, we demonstrate these points empirically for multiset CCA. We believe that our empirical comparisons of mDLAG to group factor analysis (GFA; Fig. 4, Supplementary Fig. S2, Supplementary Fig. S3) are a representative demonstration of the advantages mDLAG has over the broad class of static methods (i.e., methods that do not consider the flow of time), which includes multiset CCA, Procrustes alignment, and GFA.\n\nIn greater detail, we applied multiset CCA to both of our simulated datasets (Reviewer Fig. R2) and to the Neuropixels recordings (Reviewer Fig. R3):\n- When applied to Simulation 1 (Fig. 3), the loading matrix estimated via multiset CCA (Reviewer Fig. R2a, right) did not clearly show the group structure in the ground truth loading matrix (Reviewer Fig. R2a, left). Indeed, typical formulations of multiset CCA do not aim to distinguish these network-level interactions, unlike mDLAG or GFA.\n- When applied to Simulation 2 (Fig. 4, Reviewer Fig. R2b), multiset CCA's latent estimates represented a mixture of the two directed interactions (Reviewer Fig. R2c). GFA, a static method like multiset CCA, exhibited the same shortcoming (Fig. 4c). mDLAG successfully disentangled the two interactions (Fig. 4b).\n- When applied to the Neuropixels recordings (Fig. 5), multiset CCA was outperformed by both mDLAG (Reviewer Fig. R3a) and GFA (Reviewer Fig. R3b). These results are consistent with prior work in which GFA was shown to outperform multiset CCA on multi-area fMRI data (Klami et al., 2015), and DLAG was shown to outperform CCA on electrophysiological recordings from two brain areas (Gokcen et al., 2022).\n\n> The model assumes linearity and Gaussian observations...\n\nWe thank the reviewer for making this point. Please see the General Response to Reviewers, \"Gaussian observation model applied to spike counts.\"\n\n> Lastly, the computational cost...\n\nWe thank the reviewer for making this point. Please see the General Response to Reviewers, \"Computational demands.\"\n\n> The method seems to be a direct extension of DLAG...\n\nThe reviewer is correct that mDLAG is an extension of DLAG to multiple (more than two) neuronal populations. We directly compare mDLAG to DLAG throughout the text:\n1. In Section 2, under \"mDLAG special cases,\" we write that \"In the case of two populations (M = 2), mDLAG is equivalent to a Bayesian formulation of DLAG.\"\n2. In Section 3, under \"Validating mDLAG on recordings from V1 and V2,\" along with Supplementary Figures S1 and S2a, we make an empirical comparison between mDLAG and DLAG:\n\"mDLAG outperformed DLAG across all datasets (Supplementary Fig. S2a, points above the diagonal), suggesting that ARD provides an improved method of model selection over the constrained grid search method used for DLAG, while also avoiding grid search’s computational drawbacks (see Supplementary Section S5).\"\n3. In Supplementary Section S5, we outline the heuristic model selection approach employed for DLAG (Gokcen et al., 2022), and in the Discussion, we note that \"Scaling this type of approach to three or more populations (or external experimental variables) would be difficult. mDLAG is thus an advance toward scaling to large-scale multi-population recordings…\"\n\n> **Questions**\n\n> What prompted the specific parameter choices?\n\nWe thank the reviewer for this clarification question. Particularly in Simulation 1 (Fig. 3), we aimed to match all dataset characteristics to typical neural data. However, we chose to set 10 neurons per population to facilitate visual demonstration of the results (specifically, the loading matrices in Fig. 3a). The results do not meaningfully change if we set population sizes closer to those of the neural recordings. For example, if we re-run Simulation 1, but set 100 neurons per population instead of 10, then the accuracy of estimates (Fig. 3a,b) looks qualitatively similar. In fact, performance slightly improves: $R^2$ between ground truth and estimated latents is 0.957 (compared to 0.936 originally), and mean delay error is 0.58 ms (compared to 1.14 ms originally). We note that we further validated mDLAG on previously studied neural recordings (Section 3, \"Validating mDLAG on recordings from V1 and V2) to test mDLAG in a truly realistic setting.\n\nPlease see also the General Response to Reviewers, \"Amount of data needed.\"\n\n> Another additional application...\n\nThere is precedent for using latent time series models to study both neural activity and behavioral signals. For example, Kao et al., *Nature Communications*, 2015 developed a latent dynamical model to study the relationship between spike counts of a single motor area population and arm kinematics. They used a common set of latents (and consequently the same latent dynamical model) to describe neural activity and arm kinematics. The arm kinematics were sampled to the same time resolution as the spike count bins. Using mDLAG, we could extend this idea to multiple brain areas and arm kinematics. \n\n> **Limitations**\n\nWe thank the reviewer for these points. Regarding computational cost, we note that we report mDLAG's runtime per EM iteration and compare to GFA's runtime in Supplementary Fig. S3. Please see also the General Response to Reviewers, \"Computational demands.\" Regarding data demands, please see the General Response to Reviewers, \"Amount of data needed.\"\n"}, {"review_id": "iyQFZh6CMX", "reviewer": "Reviewer_qoeh", "summary": "This paper develops a new approach to analyze multi-region neural populations recordings. This approach extends DLAG to analyze communication across more than 2 regions. The main technical novelty lies in tractably extending the model definition of DLAG to multiple regions through the incorporation of the ARD prior. I found the approach elegant and particularly relevant to the current state of neuroscience research where multi-region recordings are growing in number. DLAG and mDLAG like approaches offer the ability to temporally disentangle concurrent signals, allowing us to make sense of multi-region recordings. The experiments are thorough and all the figures are very well made and described. Overall, I very much enjoyed reading this!", "questions": "In addition to the questions in the weakness section, I have some other minor questions:\n1. Shouldn't $\\bf{d}^m$ have a prior that requires it to be positive? I'd be curious to hear how the authors deal with this?\n2. I am curious to understand the amount of data needed to reliably estimate the parameters of mDLAG, and so I wonder if the authors have any plots demonstrating accuracy vs # of data points for a given # of regions and population size.\n\n\n\n\n\n\n\n\n\n\n\n\n\n", "limitations": "I do not envision any societal implications of this work in the short term.", "rating": 8, "confidence": 4, "soundness": 4, "presentation": 4, "contribution": 4, "strengths": "1. The paper is extremely relevant for a computational neuroscience crowd, where analyzing multi-region recordings is an important challenge. This approach can in practice help drive neuroscience research and answer questions about signaling across regions in the brain.\n2. Technically, the novelty lies in the model description (which allows for easy disentangling of latents involved across region(s)) and in discovering a tractable inference scheme (I should say that I did not read details of the inference in the supplement).\n3. The simulations and real-world experiments are both very thorough and nicely done. Details in the supplement were also helpful in contextualizing some of the results. \n4. The paper is very well-written, and the figures are very clean!", "weaknesses": " 1. The main question that I found myself struggling with was the order of latent variables across regions. I realize that there in a GP connecting them, but that does not guarantee that the ordering of latents is preserved in anyway across regions, right? Concretely, I am not sure I understand that reading through the columns of $C^m$ can tell us which dimensions are shared and which are not unless there is a post-processing step that reorders the latents (and the columns of $C^m$ respectively) so that the ordering of latents across regions is preserved (example: if latent 1 in regions 2 is correlated with a latent in region 1, then does it occupies the first dimension in the latent vector corresponding to region 1 at any time point?). Perhaps this falls our of the bayesian structure of the model, in which case it would be helpful to clarify this. \n2. It would be nice to also include a literature review / related work section discussing the connections of mDLAG to multi-region LDS based models. Including this in the introduction is also fine, but I do think that it would be helpful to describe the pros and cons of both.\n3. Finally, discussions of the amount of time and data needed to obtain reliable estimates from mDLAG would help practitioners understand the merits of the approach. Currently, there is a figure in the supplement describing the clock time per iteration, but as a reader I do not know how many iterations are needed in practice to fit the model, so having a number describing the overall runtime would be useful.  ", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "> **Weaknesses**\n\n> 1. The main question that I found myself struggling with...\n\nWe thank the reviewer for the opportunity to clarify this point. The ordering of latents across regions is preserved by definition of the model. Let us write out the observation model (Eq. 1) more expansively:\n\t$$y^m_{n,t} = \\mathbf{c}^m_1 \\cdot x^m_{n,1,t} + … + \\mathbf{c}^m_p \\cdot x^m_{n,p,t}  + \\mathbf{d}^m + \\boldsymbol{\\varepsilon}^m$$\nFrom here, we can see that latent variable 1, $x^m_{n,1,t}$, always maps to population $m$'s neural activity through the 1st column of $C^m$ ($\\mathbf{c}^m_1$), latent variable 2, $x^m_{n,2,t}$, always maps to population $m$'s neural activity through the 2nd column of $C^m$ ($\\mathbf{c}^m_2$), and so on. This ordering of the latents and columns of $C^m$ is the same for every population, by definition. We have added the following clarification to Section 2:\n- \"In particular, the $j$th column of $C^m$ maps the $j$th latent variable $x^m_{n,j,t}$ to population $m$.\"\n\n> 2. It would be nice to also include a literature review...\n\nWe thank the reviewer for this suggestion. We have added the following text to the Discussion:\n\n\"For exploratory data analysis, mDLAG's GP-based description of multi-population temporal structure is advantageous over an alternative linear dynamical system (LDS)-based description (Semedo et al., 2014; Glaser et al., 2020) in two respects: (1) a GP can be useful for exploratory data analyses where an appropriate parametric dynamical model is unknown *a priori*, and (2) mDLAG's continuous-time model enables the discovery of wide-ranging delays with high precision, which, in contrast to discrete-time LDS approaches, are not limited to be integer multiples of the sampling period or spike count bin width of the neural activity. Ultimately, these approaches can be complementary: one can use mDLAG to generate data-driven hypotheses about motifs of concurrent signaling across populations, and then test these hypotheses with a dynamical system-based approach.\"\n\n> 3. Finally, discussions of the amount of time and data needed...\n\nRegarding the amount of data needed, please see the General Response to Reviewers, \"Amount of data needed.\"\n\nRegarding the amount of time needed, we have added the following text to Supplementary Fig. S3: \"We conservatively ran each mDLAG model for 50,000 iterations, resulting in an average total runtime of 34 hours. Had we used a less conservative, but still reasonable, stopping tolerance of $10^{-6}$, the average number of iterations required for convergence would have been 17,000 (with similar parameter estimates), for an average total runtime of 11.5 hours.\" For further discussion, please see the General Response to Reviewers, \"Computational demands.\"\n\n> **Questions**\n\n> 1. Shouldn't $\\mathbf{d}^m$ have a prior that requires it to be positive?\n\nWe thank the reviewer for this clarification question. Please see the General Response to Reviewers, \"Positivity constraints on the mean parameter.\"\n\n> 2. I am curious to understand the amount of data needed...\n\nWe thank the reviewer for the interesting question. Please see the General Response to Reviewers, \"Amount of data needed.\""}, {"review_id": "vV1szgvFiW", "reviewer": "Reviewer_EMi6", "summary": "This work extends an established line of previously proposed methods for Gaussian Process factor analysis-based models of neuroscience data. In comparison with previously published methods, it extends the factor model to include correlations and time delays between different (matched) latents underlying different populations. Inference is performed via mostly closed-form variational updates. This method is applied to data from paired recordings from three regions in visual cortex, providing evidence of directional influence between them in the context of overlapping receptive fields.\n\nThis is a valuable contribution that extends previous methods in ways that are likely to prove increasingly useful in an era of widespread multi-region recordings. While the algorithmic innovation is somewhat incremental, the validation, presentation, and contextualization to questions in neuroscience are very good.", "questions": "- The authors report several latents with nearly sinusoidal form. Are inferred delays in these cases well-defined? Is it possible to infer directionality in these cases, when a lag of $2\\pi(1  - D/\\tau)$ (with $\\tau$ the period) should also produce similar results? This seems to me clearer in cases of sparse signals or strong transients but more difficult when latents are densely active.\n\n- While the authors cite a number of related factor analysis models in neuroscience (especially refs 14, 15), there is a limited discussion of some other lines of multi-region work (e.g., ref 33, but also Gallagher et al. NeurIPS 2017 from the same group) that have previously used a GP-based factor analysis model and included relative delays between regions. What, precisely, is the difference between these models apart from their application to spikes as opposed to LFP?", "limitations": "- From (7), it appears as if the latents for each $j=1\\ldots p$ are independent: $\\mathrm{cov}[x^m_{njt}x^m_{nj't}] \\propto \\delta_{jj'}$. That is, there is potential coupling _across populations_ for each latent but not _across latents_ even within a population. This choice might warrant some additional discussion in the text.\n- The variational inference algorithm for the latents requires inverting an $MpT \\times MpT$ matrix for each trial. I suspect this is feasible for typical numbers ($M \\sim 3$, $p \\sim 10$, $T \\sim 10$) and can be parallelized across trials, but will scale poorly if any of the relevant parameters become large.\n- While the model is fit to spike count data, the observation model in (1) is Gaussian. This is likely to hinder inference when bins are small and/or counts are low.  ", "rating": 7, "confidence": 4, "soundness": 3, "presentation": 4, "contribution": 3, "strengths": "- This is a valuable contribution that extends previous methods in ways that are likely to prove increasingly useful in an era of widespread multi-region recordings. While the algorithmic innovation is somewhat incremental, the validation, presentation, and contextualization to questions in neuroscience are very good.\n- The manuscript is clearly and carefully written, including a helpful notation key in the supplement.\n- Rigorous model evaluation through recovery of parameters on synthetic data and careful analysis of a challenge data set.\n- Details of the model and its inference are carefully laid out in the main text and supplement.", "weaknesses": "- While the authors have done a commendable job of using a consistent and careful notation, the resulting formulas are index-dense, and it may be easy for some readers to miss the key conceptual setup, particularly which latents are independent of which other latents. For example, l. 94 in the supplement clearly notes that latent posterior estimates are independent for different trials, and this is implicit in (7), but this probably deserves a mention in the main text. Similarly, different $\\mathbf{x}_j$ for $j=1\\ldots p$ appear to also be independent, and this is a substantive assumption (see comment under limitations below), but I didn't see this discussed (though I may have missed it).", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "> **Weaknesses**\n\nThanks to the reviewer, we realize we could have been clearer about the independence structure of the latents *a priori* and *a posteriori*. We have added the following text:\n- Section 2: \"Under equation 7, latents are independent and identically distributed across trials.\"\n- Supp. Section S2.1.1: \"Note that, under the posterior distribution, latent variables $j = 1,\\ldots,p$ are no longer independent, as they are under the prior distribution (equations 7, S2).\"\n\nWe hope that these statements better complement Supp. Section S2.1.1, in which we preface each posterior update (eqs. S11-S20) with a statement about independence structure. Please see also our response under Limitations, below.\n\n> **Questions**\n\n> - Sinusoids\n\nThe reviewer's intuition regarding purely sinusoidal signals is correct. The central quantity is the latent's underlying covariance function. For example, a signal generated from a squared exponential function (Fig. 2c) would be \"densely active,\" but the correct time delay can be unambiguously identified via the peak of the cross-covariance function. A sinusoidal signal has a sinusoidal covariance function, which has translational symmetry (i.e., it is periodic). Hence two time delay estimates spaced one period apart present two local optima. \n\nThis potential ambiguity can be resolved by additional context. For example, latent 3 of Fig. 5f is nearly sinusoidal with period 333 ms (corresponding to the grating stimulus) and estimated time delay +10 ms. An initial transient, however, breaks its translational symmetry (the 1st cycle is different from the 2nd and 3rd cycle). Furthermore, potential alternative delay estimates, -323 ms and +343 ms, would be inconsistent with the response latencies of the neurons in V1 and V3d (Fig. 5b; both populations respond within 10s of ms of stimulus onset). We can thus conclude that +10 ms is a reasonable estimate consistent with additional observations. Finally, we note that the sinusoidal signals identified here are due to experimental design: future experiments could be designed to avoid them.\n\n> - Related models\n\nWe thank the reviewer for this reference, and now cite it. In brief, we see three axes that distinguish mDLAG from the cross-spectral factor analysis (CSFA) approach of Gallagher et al., 2017: (1) group structure, (2) parametrization of time delays, and (3) automatic relevance determination (ARD). The design choices along these axes give mDLAG interpretational and computational advantages in the study of multi-population recordings.\n\nIn detail:\n1. Group structure: While Gallagher et al. indeed analyze activity in different brain regions, each \"brain region\" corresponds to a single LFP channel or time series. Consequently, there is no explicit group structure (i.e., the capability that brain regions may be multi-variate) built into CSFA, either in the observation model (eqs. 4 and 5 of Gallagher et al.) or the state model (eqs. 1 and 6 of Gallagher et al.). CSFA appears to be similar to time-delay GPFA (Lakshmanan et al., 2015), which studied time delays between pairs of individual neurons. In contrast, mDLAG incorporates group structure into both the observation (eqs. 1-6) and state models (eqs. 7-9; time delays are shared across neurons in the same population). Group structure facilitates inferences about network-level interactions (Fig. 1c, Fig. 2b) and about signal flow between populations (Fig. 1b), rather than between pairs of individual neurons or LFP channels.\n\n2. Parametrization of time delays: Gallagher et al. use the cross-spectral mixture kernel (Ulrich et al., 2015). Accordingly, (1) the GP covariance function of each latent variable is a weighted sum of spectral Gaussian kernels, and (2) for each latent variable and between each pair of LFP channels, each spectral Gaussian component contributes a constant phase difference at its center frequency. Effectively, \"time delays\" between a pair of LFP channels are parametrized as a piecewise-constant phase function in the frequency domain. In contrast, (1) the GP covariance function of each mDLAG latent variable comprises one component (e.g., the squared exponential function), and (2) for each latent variable and between each pair of populations, mDLAG directly defines a time delay, or equivalently, a linear (as opposed to piecewise-constant) phase function in the frequency domain. These design choices enable a simpler description of cross-population temporal structure (fewer parameters) while remaining reasonably expressive.\n\n3. ARD: CSFA requires cross-validation over a space of four hyperparameters. mDLAG employs ARD to enable computationally tractable model selection for multi-population recordings.\n\n> **Limitations**\n\n> - Independence\n\nIndeed, under the prior distribution, eq. 7 (see also Supp. eq. S2), the latents $x_{n,j,:}$, $j =1,\\ldots,p$ are independent. Under the posterior distribution, however, this independence no longer holds (see Supp. eqs. S6, S11, and S12). The interaction between the structure of the prior and posterior distributions can be understood by inspection of the ELBO (Supp. Section S2.2, eq. S31). The KL-divergence term $KL(Q_x(X) || P(X|\\Omega))$ penalizes deviations of the posterior distribution $Q_x(X)$ from the prior distribution $P(X|\\Omega)$. The independence structure of the prior distribution therefore acts as a form of regularization, not a hard constraint.\n\n> - Variational inference\n\nWe thank the reviewer for the opportunity to clarify an important point: the latent posterior covariance (Supp. eq. S11) is identical for all trials of the same length. This computation can therefore be reused efficiently across trials. We have added this statement to Supp. Section S2.1.1. Please see also the General Response to Reviewers, \"Computational demands.\"\n\n> - Gaussian model\n\nWe thank the reviewer for making this point. Please see the General Response to Reviewers, \"Gaussian observation model applied to spike counts.\"\n"}]

Author rebuttal:
## General Response to Reviewers

We thank the reviewers for their constructive comments, which helped to strengthen our submission. Here, we provide responses to comments shared by multiple reviewers.

### Amount of data needed

Reviewers qoeh and DL1P inquired about the amount of data needed to fit an mDLAG model. In general, the answer to this question is highly data-dependent: the number of trials, number of neurons, number of populations, trial lengths, latent dimensionality, latent temporal structure, and noise levels all interact to influence mDLAG's performance. However, we can address this question with an empirical example (Reviewer Fig. R1). We re-fit mDLAG models to the Neuropixels dataset analyzed in Fig. 5, but we limited the number of experimental trials available in the training set (i.e., we fit mDLAG to training sets with 10, 25, 50, 100, 150, up to 225 trials—the full training set size). With as few as 100 training trials (less than half of the full training set size), mDLAG's test performance suffered by little more than 5% (Reviewer Fig. R1, black curve). Furthermore, with as few as 25 training trials, mDLAG still outperformed the group factor analysis (GFA) model fit to all 225 training trials (Reviewer Fig. R1, red dashed line).

This example demonstrates empirically that mDLAG performs well with trial counts typical of neurophysiological experiments. In principle, this data-efficiency is due to mDLAG's incorporation of (1) dimensionality reduction, (2) temporal smoothing, and (3) automatic relevance determination. These components act as forms of regularization that benefit model performance in data-limited regimes.

### Computational demands

Reviewers EMi6 and DL1P have pointed out that temporal smoothing via a GP prior is computationally demanding. We agree, and acknowledge in Supplementary Fig. S3 (where we compare mDLAG and GFA runtimes) that "Overall, mDLAG is more computationally intensive than GFA due to the incorporation of Gaussian processes. Each mDLAG fitting iteration requires the inversion of a $MpT \times MpT$ matrix (equation S11)."

We emphasize, however, two points from the Discussion: (1) mDLAG provides a tractable approach to model selection for neural recordings involving three or more populations, a problem that was previously computationally prohibitive; (2) The problem of scaling the latent inference step for this class of Gaussian process latent variable model has been studied extensively (Gallagher et al., 2017; Zhao et al., 2017; Duncker et al., 2018; Keeley et al., 2020; Jensen et al., 2021; Dowling et al., 2023). Approaches include diverse uses of Fourier approximations, inducing points, and additional variational approximations. mDLAG is compatible with these approaches, and we believe that they provide a promising avenue for continued computational development of the mDLAG framework.

### Gaussian observation model applied to spike counts

Reviewers EMi6 and DL1P have pointed out that mDLAG incorporates a Gaussian observation model, whereas the spike counts to which we apply mDLAG are generally better described by a Poisson observation model. While we agree with this point, we note that this limitation is not unique to mDLAG, and is shared by widely used methods such as factor analysis, Gaussian process factor analysis, and latent linear dynamical systems. Intuitively, one can think of the Gaussian model as capturing the first and second moment of the spike count data, rather than requiring that the data be Gaussian-distributed. In our experience, for neural recordings similar to those analyzed in this work, we have not encountered a scientific finding that could be seen with a Poisson or point process observation model that could not also be seen using a Gaussian model. Hence mDLAG is still widely applicable to spiking neural activity. We also emphasize, as we have in the Discussion, that mDLAG's Gaussian observation model enables wider applicability to other neural recording modalities. For spiking neural activity in which a Gaussian observation model is meaningfully limiting, a Poisson or point process observation model can be successfully incorporated, but at the expense of increased computational requirements (Zhao et al., 2017; Duncker et al., 2018; Keeley et al., 2020; Jensen et al., 2021; Dowling et al., 2023).

### Positivity constraints on the mean parameter

Reviewers qoeh and DK8X have inquired about potential positivity constraints on the mean parameter $\mathbf{d}^m$. We thank the reviewers for this question, as we realized we needed to further clarify our data preprocessing steps. Introducing a positivity constraint on $\mathbf{d}^m$ could make sense for applications to neural spike counts, especially if it is particularly important to interpret each element of $\mathbf{d}^m$ as the mean firing rate of each neuron across time and trials. However, a positivity constraint might not make sense for analyses of different recording modalities, or given certain preprocessing choices of spiking neural activity. For example, in this work, because we were interested in inter-population interactions on timescales within a trial, we subtracted the mean across time bins within each trial from each neuron. This step removed activity that fluctuated on slow timescales from one stimulus presentation to the next (Cowley et al., 2020). Consequently, the data input to mDLAG were no longer positive-valued spike counts. We have added the following text:
- Section 3: "Because we were interested in V1-V2 interactions on timescales within a trial, we subtracted the mean across time bins within each trial from each neuron. This step removed activity that fluctuated on slow timescales from one stimulus presentation to the next (Cowley et al., 2020)."
- Section 4: "As we did for the V1-V2 recordings, above, we subtracted the mean across time bins within each trial from each neuron, to remove slow fluctuations beyond the timescale of a trial."


Meta-review: This work provides a computational model for identifying the network interactions between multiple brain areas. The work builds off of a 2-area model, delayed latents across groups (DLAG), which finds low-rank representations with temporal delays between brain areas, representing the expected propagation time across neural populations. The model establishes a Bayesian fitting procedure that generalizes on the past DLAG method, using Gaussian Process latent states and ARD determination of neural subsets corresponding to each state. The method is shown on simulated data as well as real visual cortex data across three neural populations recorded with SOTA neuropixel probes.

The reviewers universally found this work interesting and timely with the increased capability to record from multiple brain areas simultaneously. The main remaining weakness is one of scaling (how will the method scale in computation and interpretation to 40+ brain areas?). Regardless this is a strong paper that I recommend for acceptance to NeurIPS.
Decision: Accept (spotlight)