Modeling Neural Activity with Conditionally Linear Dynamical Systems
Authors: Victor Geadah
Amin Nejatbakhsh
David Lipshutz
Jonathan W. Pillow
Alex H Williams
Venue: NeurIPS 2025 poster ()
Keywords: Linear Dynamical Systems
Gaussian Processes
Neural data analysis
Expectation-Maximization
Primary area: neuroscience_and_cognitive_science
TL;DR: We extended classical linear-Gaussian state space models of neural circuit dynamics to capture nonlinear dependencies on experimental conditions, while maintaining ease of fit and interpretability.

Abstract: Neural population activity exhibits complex, nonlinear dynamics, varying in time, over trials, and across experimental conditions. Here, we develop *Conditionally Linear Dynamical System* (CLDS) models as a general-purpose method to characterize these dynamics. These models use Gaussian Process priors to capture the nonlinear dependence of circuit dynamics on task and behavioral variables. Conditioned on these covariates, the data is modeled with linear dynamics. This allows for transparent interpretation and tractable Bayesian inference. We find that CLDS models can perform well even in severely data-limited regimes (e.g. one trial per condition) due to their Bayesian formulation and ability to share statistical power across nearby task conditions. In example applications, we apply CLDS to model thalamic neurons that nonlinearly encode heading direction and to model motor cortical neurons during a cued-reaching task.

Modeling Neural Activity with
Conditionally Linear Dynamical Systems
Victor Geadah
Amin Nejatbakhsh
David Lipshutz
Jonathan W. Pillow!.
Alex H. Williams
Program in Applied and Computational Mathematics, Princeton University, Princeton NJ
Center for Computational Neuroscience, Flatiron Institute, New York City NY
Department of Neuroscience, Baylor College of Medicine, Houston TX
*Princeton Neuroscience Institute, Princeton NJ CCenter for Neural Science, New York University, New York City NY
Abstract
Neural population activity exhibits complex, nonlinear dynamics, varying in time, over trials, and across experimental conditions. Here, we develop Conditionally Linear Dynamical System (CLDS) models general-purpose method to charac- terize these dynamics. These models use Gaussian Process (GP) priors to capture the nonlinear dependence of circuit dynamics on task and behavioral variables. Conditioned on these covariates, the data is modeled with linear dynamics. This allows for transparent interpretation and tractable Bayesian inference. We find that CLDS models can perform well even in severely data-limited regimes (e.g. one trial per condition) due to their Bayesian formulation and ability to share statistical power across nearby task conditions. In example applications, we apply CLDSS to model thalamic neurons that nonlinearly encode heading direction and model motor cortical neurons during a cued reaching task.
Introduction
A central problem in neuroscience is to capture how neural dynamics are affected by external sensory stimuli, task variables, and behavioral covariates. To address this, longstanding line of research has focused on characterizing neural dynamics through recurrent neural networks (RNNs) and their probabilistic counterparts, state-space models (SSMs; for reviews, see 2
Early work in this area utilized latent linear dynamical systems (LDS) with Gaussian observation noise. Although these assumptions are restrictive, they are beneficial in two respects. First, they simplify probabilistic inference by enabling Kalman smoothing and expectation maximization (EM)- two classical and highly effective methods (4) Second, they produce models that mathematically tractable to analyze with well-established tools from linear systems theory [ Indeed, many influential results in theoretical neuroscience have come from purely linear models
In reality, most neural circuits do not behave like time-invariant linear systems. Thus, more recent work from the machine learning community has cataloged variety of nonlinear models for neural dynamics. Although these new models often predict held out neural data more accurately than LDS models, they generally more difficult to fit and more difficult to understand. Thus, there has been a proliferation competing models for neural data analysis, such as RNNs ( transformers [10 [1] and -effussio-bbseed methods [12), well as competing inference methods, such generalized
39th Conference on Neural Information Processing Systems (NcurIPS 2025).

---

Trial-fevel conditions
Conditionally Linear Dynumical System
Ieall mnail
& MM
condicions
latent stutes
Condition space
Time- come--arrying conditiions
Afline
dyramio
Alline
relicom:
neural Cosa
Blacb-lloo
Cine
Figure 1: (a) Neural dataset consisting of spike trains collected over multiple trials, along with corresponding experimental conditions. (b) Conditionally Linear Dynamical Systems are linear in state-space dynamies and capture nonlinear dependencies over conditions. Shaded nodes observed, clear nodes are latent.
teacher forcing [3] amortized variational inference ( and sequential Monte Carlo [4] among these strategies and scientifically interpreting the outcomes is challenging.
Choosing
Here describe Conditionally Linear Dynamical Systems (CL.D)) framework to jointly capture some of the benefits of the classical (i.e. linear) and contemporary (i.e. nonlinear) approaches to modeling neural data. CLDS models parametrize a collection of LDS models that vary smoothly as function of an observed variable u, (e.g. measured sensory input or behavior at time t). Assuming the presence of U, often feature-not bug-of this approach. Indeed, common goal in neuroscience is to relate measured sensory or behavioral covariates to neural activity, Additional features of the CLDS framework include:
1. CLDS models are locally interpretable. Conditioned on u, the dynamics are linear and amenable number of classical analyses.
2. CLDS models are easy to fit ($2.) Gaussian noise is assumed, then exact latent variable inference (via Kalman smoothing) and fast optimization (via closed-form EM) is possible. Under more realistic noise models (c.g. Poisson), the posterior over latent state trajectories is still log concave and amenable to relatively fast and simple inference routines.
3. CLDS models are expressive. As uy changes, the parameters of the linear system are allowed to change non- non-linearly. Thus, CLDS can model complex dynamical structures such as ring attractors (3.2). that are impossible for vanilla LDS to capture.
4. CLDS models are data efficient. To the extent that LDS parameters change smoothly as function of ut, few trials per condition ($3.5). In fact, CLDS models can interpolate to make accurate predictions entirely unseen conditions.
the recover the parameters of the dynamical system with very
Finally, CLDS models have connections to several existing methods ($ For example, they can be viewed as dynamical extension of a Wishart process model [15] and extension of Gaussian Process Factor Analysis (GPFA; (16)) with a learnable kernel and readout function that can vary across time and conditions. They also similar to various forms of switching linear dynamical systems models [17 I8, 19]. The key difference is that the "switching" in CLDS models is governed by an observed covariate vector, u4, rather than by discrete latent process. This makes inference in the CLDS model much more straightforward, albeit at the price of not being fully unsupervised.
2 Methods
X' x -> R, to denote samples f : X - R. obtained from stacking independent Gaussian processes into an N d-dimensional vector. Also, Ip is the D D identity matrix.
Notation We use f(.) ~ GPNV GPN(GpPN)))m...m k(., -)) with mean m X - R and kernel :

---

2.1 Conditionally Linear Dynamical Systems
Consider an experiment with N neurons recorded over K trials of length T. Our dataset consists of the recorded neural firing-rate rate trajectories (yir)
LK with E R at time-step t E .1.....Ty T,, 1K along with the corresponding experimental conditions TuiTkk= with u, in the condition space U. All modeling done for single trials, and we drop superscripts & where clear henceforth. By experimental conditions, we refer to available neural data covariates, either experimentally set or collected as measurements. These conditions, see Fig. Th. can vary over time or remain constant (compare $3.4 vs. (..5). They can also be step function over time (e.g. animal moving vs. not moving), resulting in a switching-like mechanism between different dynamics.
We model responses y: as emissions from a latent time-varying linear dynamical system in X E IR with dynamics governed by the conditions uy. Specifically,
X++1 = A(,))x + b(u)) + c2
(1a)
y; = C(u,)x + d(u,) + we
(lb)
over time steps t E {1....,T) Th} from initial condition X₁ N (x1; (())) Q1), and where E, (0, Q) and wi ~ (0, It) are sources of noise, i.i.d. in time. We assume that the latent variables X₁ follow smooth dynamics defined by time-narying linear matrices A(u) E DDD from initial mean m m(u) (1) E RD, with bias terms b(u) E R., and emissions governed by C(u) E RXD and d(u) E R.V all dependent on E U. See graphical depiction in Fig. Ip.
The system in Ta Tbp parameterizes a family of linear systems indexed by a continuous variable u4 which importantly is observed. Our approach can be viewed as the linearization in x, of a fully nonlinear system x; and ut, under additive noise-we we detail this relationship in Appendix $A.2 Most methods treat conditions as additive inputs, influencing the dynamics in (Ta) via additive terms of the form Bu, for linear encoding matrix B E IR Dx|[1] In contrast, the mapping of experimental conditions onto linear dynamics, ト (A(u), b(u), C(u),d(u) m()), }. is allowed to be nonlinear and learnable. Specifically, we place an approximate Gaussian Process (GP) prior on each entry of the parameters through a finite expansion of basis functions, leveraging regular Fourier feature approximations 20). For any M E {A, b, C, d, m), we consider a prior
(0) (e()), M,y Σ We
(if) vid 1)0,
(Ic)
over each j,j-th entry, truncated at L E N basis functions. Each basis function Φi 4 + I is fixed, and the randomness in the prior purely comes from the weights, w which are drawn from standard normal distribution. When constructed appropriately, the prior in equation (Tc) converges to a nonparametric Gaussian process in the limit that L - oc, with kernel structure determined by the basis functions. Many choices are valid, and we elect to choose basis functions { with the goal of expressivity; they are designed to approximate a GP prior of the form (-) ~ GP(O, k.)) for the squared exponential kernel k. with variance o² and length-scale [21 eq. 13].
We denote F- = {A, b, C, d, m} as the set of random functions, and analogously the parameter set F(u) (u),, buu), C(u), d(u), m(u) for any condition u E U. The model distribution
p(yI.T, XI.T A, b, C, m, uI.T) = ply 1.7, X F, u..T)
describes a time-varying LDS, conditioned on parameter sequence set at experimental conditions. Therefore, we refer to the model as a Conditionally Linear Dynamical System (CLDS
2.2 CLDS modeling choices
Practitioners can adapt a CLDS model in several ways to suit different applications and modeling assumptions. First, the GP prior can be tuned to trade off model expressivity for interpretability and learnaillit.. In one extreme, as let function of u and become independent per u. In the other extreme, if one takes λ + co, then the LDS parameters become constant (do not change as a function of u) and we recover a time-invariant LDS model with autonomous dynamics. In this regime, we could also modify the GP prior over b(-)
0, the LDS parameters change rapidly, nonlinearly
OOur CLDS implementation is available at httppp////gthhb..c.co//o/essssss
3

---

to follow a linear kernel, k(u, u²) = -u u', resulting in time- invariant LDS with additive dependence on u, Thus, CLDS models capture classical linear models as a special case, Moreover, the model's prior can be tuned to capture progressively nonlinear dynamics.
A second source of flexibility is the encoding of experimental covariates, u. Recall notation from section 2.1 that u represents experimental covariates at time E {l,..., ,T and trial & E f1...., J. A simple, and broadly applicable, modeling approach would be to set u = t. This achieves time- vvrying LDS model in which the GP prior encodes smoothness over time. This is similar in concept to fitting a linear model to data over a sliding time window 22 23 However, the CLDS formulation of this idea is fully probabilístic, which has several advantages. For example, one can use a single pass of Kalman smoothing to infer the distribution over the latent state trajectory, x within each trial. It is comparatively non-trivial average latent state trajectories across multiple LDS models that are independently fit to data in overlapping time windows.
In section 3 we demonstrate more sophisticated examples where už is specified to track continu- ously measured behavioral variable (e.g. head direction or position of an animal) or follow a stepping or ramping function aligned to discrete task events (e.g. sensory "go cue" or movement onset). Section 4 discusses further connections between CLDS models and existing state space models.
2.3 Inference
As mentioned earlier, the conditional distribution in eq. (2) has the advantage of describing a latent LDS, or linear Gaussian state-space model, given estimates of F. As such, we can benefit from analytic tools like Kalman filtering to compute the filtering distributions p(xe Yise; F, u1:)) and marginal log-likelihood p(yI:T F, uI:T), and Kalman smoothing for p(x; YI:T, F, UI:T) and the latents posterior mode X1:r. We focus on performing maximum-a-pont a-posteriori (MAP) inference for these parameters. In principle, it would be a straightforward extension to use variational inference or Markov Chain Monte Carlo to approximate the full posterior over these parameters.
Conditionally Linear Regression As stepping stone towards our goal of performing MAP inference for {A, b, C, m), consider the model
(3) given data yn E RRP:, regressors x,, E ERD: conditions un E 44 repeats n E (1....,,N) N). noise
Yn = M((,))x, + En, Cn E(,,E,, M(.) ~ GDDixDs (0, ku),
covariance Σ > 0, and with an approximate GP prior M parametrized in (Ic). We refer to
the model (3) as conditionally linear regression, and our goal is to perform MAP inference for the
weights { {w/ Jijjk H.,jj, for M.
Our approximate basis representation in eq. implies that each entry is a dot product, M,y()) (6() (()), where ((-) = (0)()) Therefore, we have
(-), ..., (1)()) R is our vector of basis-functions evaluations.
M(u) = W (Φ(u) & X),
for any u E V and vector or matrix X E Pa... of appropriate dimension, with "ø" the Kronecker product, and with our weights aggregated into the matrix WD.e+f, = w (t)) W E Dx Derivations are provided in Appendix $A.I.I. With this, we can rewrite our regression problem as
Yn M(u,)x,, 6,, W 2n +
with z, ((u,)) & X, RDL Thus, we have reformulated our original model into Bayesian linear regression in an expanded feature space. Namely, the MAP estimate of the weights, WMAP, is given by
arg max logp(yi:N W, XI,N, uINN) ) + log P (W)
(5)
which is a linear regression problem with regularization from the prior logp W) (up to an additive constant) from (c]. We can analytically solve for the solution (derivations in Appendix $A.1.2). which yields that WMAP is the solution to the Sylvester equation
|
z + WE - Z Y
(6)
with Z E R NxD our matrix obtained by stacking {z, 1. and similarly for Y E NxD, We that if = for some o 0 then we obtain back the familiar looking penalized least squares estimate WMAP (Z Z + o -1 ZTY Y.

---

Lacasiaod demsanics
Nocliaca: Usrraniss
Pynsmicx A(9)
Heal frott
Blon b(6)
Single reconsiatiion
Compostee Symaniss
Lace-hratio
egewwakss
smoothing
Cous Fractice
Restsfoctios
Figure 2: Head direction synthetic experiment. (a) Schematic of latent dynamics and neural activity about E [0, 2x), ) the mouse HD, serving as conditions u = 0 in this task. (b) True nonlinear flow- field corresponding to the schematic in a, computed considering p(0 x) = 60). (c) Recovered composite dynamics by CLDS (see text). Grey scale indicates posterior X2 occupancy, and thus our ability to estimate the dynamics. The model fixed points (colored) as a function of form perfect ring, overlapping with the true fixed points. (d-e) Parameter recovery for the dynamics matrix A (Ø) (d) and the bias b(0) (e) as functions of head direction 0. (f) Recovered eigenvalues of A (() function of 0, true in dashed. (g) Co-smoothing reconstruction from the test-set. The firing rate of one neuron is held-out (bottom) while the rest (top) is observed, and reconstruct accurately the single-neurod firing rates for both the held-in (top) and held-out (bottom) neurons,
Inference with Expectation Maximization We can leverage the above to perform MAP inference for F {A, b, C, d, m} with the Aopectatioon-Maxiizaaaiee (EM) algorithm [24 In the E-step we obtain estimates of the moments of the latents with dalman-smoothing which then place us in setting akin to eq. (3) with sufficient statisties data and regressors. We can then perform closed-form M-steps with our updates in eq. (6). We provide in Appendix 1.3. example of the associated derivations with these E- and M-steps for the joint update for A and b(.). The resulting EM algorithm has several advantages: (1) all E- and M-steps are analytic, (2) the E-step provides us with exact (penalized) marginal log- likeeliooood calculations, and (3) the algorithm gives monotonic gradient ascent guarantees of the marginal log- aigliiklihood (resp. log posterior) objective.
We initialize the EM algorithm at samples from our GP priors for F. With the EM algorithm also learn the covariance parameters (Q1, Q, R}. The hyper-pane parameters {L, Ki o for the GP priors and the dimensionality of the latents D are determined through performance on held-out test sets from 80/20 trial splits on all experiments unless specified.
Extensions to non-Gaussian likelihoods The closed form EM updates are only applicable when the observation distribution of y, conditioned x; and u, (Ib) is Gaussian. This assumption is common to existing methods (e.g. [16 255). though to model spike count data directly, most work utilize Poisson (e.g. (26)) and COM-Poisson (27) likelihoods. Such emission likelihoods of interest have log concave distributions, for which inference in CLDS models would remain tractable. Indeed, conditioned on uI.r, the log posterior density associated with X1:7 equal to a sum of concave terms up to an additive constant [I]. Thus, we can use standard optimization routines to identify MAP estimate of x1.7 with theoretical guarantees, and leverage it for approximate EM (26).
3 Experiments
3.1 Setup
Metries For a given trajectory {y1,7; U1:7), we denote as data reconstruction the mean emission E [YI.T &1., F(uIr)) (C(u))xi,.... C(ur) from the posterior mode X computed with Kalman smoothing, given the observations yı.r and parameters F( As our primary metric, we use co-smoothing 228 to evaluate the ability of models predict held-out single-neuron
5

---

Taskschemaac
Turing curves C 16]
(d) eigannltuss
Composite dynamics
Single neuron reconstructon d
Fixed point dynamics
Figure 3: (a) Schematic of mouse foraging in an open environment. We have access 0₂ E [0, 2x) the mouse HD in time f, which we use as conditions u, just like Fig. 2 (b) Model reconstruction on the whole dataset recovers the true data. We plot single-neuron traces, averaged over 10s trials. (c) Model tuning curves over head direction 0, obtained as C(O) [x| 00, recover the empirical tuning curves. Plotted for the top three units in firing rate norm. (d) Dynamies around each fixed point in x-state space as a function of head direction 0, with the solid-line representing the complete set of
fixed points. (e) Eigenvalues and angles of eigenvectors of (0) as function of 6. (I) Composite dynamics in xp-spac,, with overlaid colored model fixed points as a function of O.
activity. Specifically, for the top 5 neurons with highest variance from the test set, we compute the coefficient of determination Ra between the true and reconstructed single-neuron firing rate, obtained by performing data reconstruction using only the other neurons.
Composite dynamics The latent dynamical system Ia) depends on the condition uy, which can make visualizations challenging. Building on the idea that CLDS models decompose nonlinear dynamical system into linearizations governed by u (see App. $A.2), we aim to approximate the general nonlinear system by marginalizing out uy conditioned on x,- That is,
(7) which we define as the composite dynamical system. Intuitively, we expect this to provide a good approximation to the underlying nonlinear dynamies when u, and x; tightly co-determine each other-i.e., when the encoding p(x, u,) and decoding p(u; x2) conditional distributions have low variance (see App. $A.2.2). In practice, we estimate the expectation (7) by averaging over the un associated with binned values X,, of the posterior mode given a trajectory {yI.r, U1.7 J. Thus our ability to accurately estimate the flow-field around a given point under this method depends on how many posterior samples pass by it, which we report alongside the composite dynamics plots.
X gxx)) + = Ep(()x)) [A(u)x; + b(u)] +G,
Model baselines For model comparison, more flexible gpSLDS [19] (implemented in way to include linear boundaries to encompass the iSLDS (29), both with additive inputs of the form Bu, in the latent dynamics, and (3) LFADS 9 model with controller inferred-input inputs as fully nonlinear alternative. For all, we consider Gaussian observation model fit directly to the firing rates. See Appendix [ implementation details.
use baselines (I) standard LDS model and (2)
3.2 Synthetic head-direction ring attractor
We start by considering a synthetic experiment of head direction neural dynamics. We conceptualize latent dynamics that capture the head direction (HD) of the animal, with attractor dynamics about a HD-dependent fixed point-see schematic in Fig. 2h. This synthetic experiment is designed to represent nonlinear system decomposed linear systems, per HD serving the condition. We plot in Fig. 2b what the resulting, "composite dynamics" (see $3.1p. nonlinear flow-field would be, assuming the latent state encodes the head direction exactly. The generative dynamics are an instance of CLDS model by construction, so this example allows us to explore recovery performance.
Concretely, let 0g [0, 2T) denote heading direction time step f, which treat our conditions u := 8,. To build ring attractor, we parametrize two orthogonal unit vectors
e1 (0) (cos(O) sin( (2) |- sin(O) cos(O))
(8)
that describe the position on the ring and the tangent vector respectively. We design (i.e. impose) that the system converges to stable fixed point at (0), and at head direction 0 we approximately

---

integrate speed input along the subspace spanned by e2(O). To do this, we define A to be a rank-one matrix that defines leaky line attractor, with attracting (i.e. contracting) dynamics along the orthogonal ei (0). For a hyperparameter 0 < < define:
A = (1 - ())) b(0) = e
Completing the model description, we assume that the firing rate of individual neurons is given by linear readout as Yr,i (O.) + w. for neuron at time t, with d(u,) 0 and C.., is sinusoidal bump tuning curve function (see App. $B.2). Finally, we sampled trials of length T = 100 with N 10 neurons, generating the evolution of the heading direction as a random walk, 0, ~ (i 0.5?), and initialize at the origin X1 ~ (0,
We report our recovery results in Fig. 2 fixing the decoding matrix C(-) to known value as avoid non-identifiability considerations. We refer to App. $B.2 for recovery plots of C when fitted. First, observed that can generally recover the nonlinear flow-field, plotting in Fig. 2 the composite dynamics obtained from the posterior trajectories. Second, for a more parameter- based account of the recovery, we plot in Fig. 2d-e the varying biases b(0) and dynamics matrices (0) functions of the head direction O-we recovered with high-fidelity the true parameters. This recovery translated into the properties of the dynamies such the recovered eigenvalues of (0) in Fig. 2 Third, we observed that the test data single-neuron reconstruction 2g) recovers the true observations, and
the model was able accurately (R² 0.86) reconstruct held out neuron from this test-set through so-ommoohhing. Finally, we report in Appendix şB.2]a a study on the impact of observational noise on our quality of fit-we find reliable recovery even as the signal-to-noise ratio degrades.
3.3 Neural circuit model of ring-attractor dynamies
We provide an alternate synthetic ring attractor experiment based on the neuroscience literature 33) It is meant to test the CLDS and its inference in synthetic data modality with known underlying computation but which, importantly, is not generated from a CLDS model.
We write model of continuous ring attractor dynamics with bumps of activity integrating angular velocity [3T] In this model, high-dimensional and non-normal attractor dynamies are implemented through -ippereferettial units arranged in "ring", with short-range excitation and Long-range inhibition forming a bump of activity for a specific unit in a specific HD, and with a HD velocity vtt) = @(t integrator causing a shift in the location of the bump. Let y(o, t) denote the network activity of cells with preferred cead-direction & [0, at time t. The dynamics follow the stochastic partial diffcrential equation
θ at y - (-+ -y+ (( -y+ + 06
+ +
(10)
for (Ø,, t) a white noise process on the ring. Here, f is a nonlinearity, w "Mexican-har hat" convolution filter, and 7 > 0 time constant. For implementation purposes, y is discretized to E R.. N = 32, for E (1,..., T). with each unit [y]; having preferred directions Φi regularly spanning the interval [0, 2mr. See model and implementation details in Appendix BB.3) We use our observations and head-direction as conditions again, generated the same way as in our previous synthetic HD model.
The underlying model captures nigh-dimensional ring attractor dynamics. When fitting a CLDS model, we expect the low-dimensional latent dynamics capture the ring-attractor structure and the nonlinearity of the dynamics, with a linear emission model bringing it back to the high-dimensionall activity. This precisely what we found when fitting the CLDS- see Appendix $B.3] The results make us confident that the CLDS can capture this nonlinear structure even under model mismatch.
3.4 Mouse head-direction dynamics
Next, we turned to the analysis of antero-dorsal thalamic nucleus (ADn) recordings from Ref. [32 of the mouse HD system in mice foraging in an open environment (Fig. 3a.. We considered neural activity from the "wakee" period, binned in 50ms time-bins, then processed to firing rates and separated into 1Os trials. As with the synthetic experiment of the previous section, we treat the recorded head-direction 0g [0, 2w) as conditions u, O,.
We recovered single- neuron firing rates with high accuracy (Fig. 3D) through data reconstruction. We further validated our fit by computing the empirical tuning curves, which our model recovered almost
7

---

b
Conditions (2D)
Delay
1200
1339
Augk
0110
Go cue
6:35
L466
other
6.75
Reach
Reach
6.00
Dallay Go
Reach
Reach
Inferred dynamics
points Dalay A Reach
Co smoomhhtiig reconstmuction
e
Co smoothing score
1.
0.15
wuni LIADS
w.2
200
200
-202
300
Projected latents
Figure 4: Macaque reaching experiment. (a) Task schematic (left) and oopulation-averagee firing rates per trial (right). (b) 2D conditions, with trial orientation 0, and reach variable 2 E {0, 1 switching ramp onset. (c) 3D projection of the 5 dimensional latents used, projected as to align best with condition decoding. We show the model fixed-points per reach angle 8 and reach condition 2, plotted over posterior mean trajectories per trial. (d) Co-smoothing reconstruction of single held-out neurons from the test-set. (e) Co-smoothing R² per model as a function the number of trials used per reach angle during training. Error bars indicated std. around the mean over initialization seeds.
exactly (Fig. 3c). The model tuning curves are given by
y 0 E ly: x, 0]] C.. (0) 0),
(11)
which follows from the law of total expectation. The later quantity [E[x 0) represents the expected encoding of the conditions 0, which we estimate by averaging posterior trajectories, obtained with Kalman smoothing, over (binned) given corresponding firing rates.
Finally, we analyzed the learned latent dynamics. Like the synthetic example, we identified ring attractor structure (Fig. 3d). Unlike the synthetic example, per the eigenvalues of A in Fig. 3e. this ring attractor is composed of HD-dependent fixed points as opposed to line attractors.
3.5 Macaque center-out reaching task
Finally, analyzed neural recordings of dorsal premotor cortex (PMd) in macaques performing center-out reaching task (Fig. 41) from Ref. [33]. In contrast to the previous experimental conditions, we consider here two-dimensional conditions u, (@+ (2)) where 0h [0, 2xt) is the instructed reach angle, constant per trial k; and 2, E (0. 1} indicates the task reach condition (see Fig. 44) set at O during the delay and 1 at 100ms past the go-cue, at the onset of the movement-related firing rate ramp Fig. (rightt. Discrete-valuce conditions, such as the reach onset 2 E {0, 1}. are considered supported continuous interval. The correlation between such discrete points is determined by the length-scale parameter k, which we've set to x = 0.5 from a hyperparameter search. More details be found in Appendix $BB Finally,
hyperparameters and prproprepooesssing emission matrix C' and let the latent dynamics capture the dependency on experimental conditions through (u)) and b(u).
fixed
We found the latent dynamics to encode the conditions through attracting fixed-points during both the delay and reach periods. We show in Fig. 4e the projection of the D 5 latent dimensions along the 3-dimensional subspace most aligned (i.e. best decoding) with the experimental conditions, following
common analyses [33] we observed clear aligned rings of fixed points from delay reach. In CLDS models, we obtain the fixed points by simply solving for (u() satisfying (I A(u))=* b(u)
for any u. in contrast to numerical fixed-point methods usually employed (34)
We performed co-smoothing (see $3.1) to evaluate the model. We recovered with good accuracy single held-out neurons from the validation set excluded from training (Fig. 4d). We then compared

---

the performance of the CLDS against the baseline models, exploring further how each fares in low-data regimes. We report in Fig. 4 the co-smoothing R² per model, computed as a function of the number of trials used in each reach-angle 0², averaged over 5 random seeds. We found that the CLDS outperformed other models, displaying the highest difference at training trial per condition. Introducing more training trials per condition, the gpSLDS model sometimes approached the performance of the CLDS. While the LFADS model showed progressively better performance that did not plateau yet, it nonetheless underperformed in these low data regimes.
4 Related work
Wishart process models CLDS models capture the dependence of neural responses y1:7 on continuous experimental conditions uIT Ref. [15] investigated very similar problem, focusing on single-trial responses Yı to continuous experimental conditions ug E U. They use a conditional Gaussian model for responses Y& given conditions uk
Yk UA ~ (y: E(ux))
(12a)
and they place Gaussian process and Wishart process
35 priors on the mean and covariance functions
m(.) ~ GpM (0, kw),
E(u) = Uu) + A(u)
(12b)
M (0, kx) for chosen kernel functions {ky Ks, h The nlyer-parameteer p E N determines the low-rank structure of Σ.
with U(-) ~ GpMxp(O, ky) and A(-) ~ GPM
For a single time step = T = and assuming a degenerate prior m(u;) = 0, the marginal distribution of y1 conditioned u1 in our system (I) reads
N (d(u:), Q(ui)) -R) -
(13)
which can be compared with 12. The models are analogous with u(u) = d(u) and with the CLDS emission matrix C(u) serving as the Wishart process prior decomposition matrix U(u), right-scaled by Q//² E DDDD This makes the parameter p = now bear meaning as the dimensionality of the latents x E R.,, akin to factor analysis. Thus, we can view CLDS models as direct extension of Wishart process models that capture condition-dependenl dynamics across multiple time steps.
such as GPFA 166, are widely used in neuroscience. One can view MAP inference in CLDS model as optimizing kernel that defines a latent GP prior. Taking this a step further, the stationary dynamics of an AR-1 process (i.e. linear dynamical system) can be expressed as draws from a GP [16] 38). Vice versa, all stationary, real-valued, and finitely differentiable GPs admit a representation in terms linear state space models 39
Markovian GPs Latent GP models 36 37
40 a relationship
We expand upon this direction by allowing the dynamics to vary not only across time but over conditions too. For fixed set of LDS parameters at experimental covariates, F F(,i)), the distribution of latent states in a CLDS are jointly Gaussian. Thus, a set of LDS parameters induces (generally non-stationary) GP prior on the latent trajectories. In this view, the GP prior we place over the parameters of the LDS can be seen as a hyperprior over the latent dynamical process.
that has proven useful for both modeling and aiding GP inference 440
41
Switching dynamical systems A second class of relevant models generalizing the LDS are Switch- ing LDS models (SLDS; 43 17).. SLDS models consist of a discrete latent state z₁, with
finite Markov chain dynamics, dictating the dynamics matrix A (z4) This switching behavior can be mimicked in our setting if the condition space is discrete (see, e.g., $3.5). We can take the relationship step further by embedding the discrete process in the continuous dynamics parameter space of A (44) E DDD In a similar line of thinking as with Markovian GPs, we show in Appendix $A.3 how dinediinall SLDS model of dynamics
plze+1 i 2 = -) Ptj, Xe+1 x = (1(( + Et,
is equivalent, up to the first two moments of the stochastic process a = a", to CLDS model with
a(-) GP ( a, T (ilt क
over time conditions u, = I, for a the vector of values taken by (=.) and . the stationary distribution of the z₁ discrete state process. In similar vein, previous work has considered the discrete states

---

dictating the dynamics to live on a continuous support [45], however without leveraging this continuity the parameters (-) themselves.
The recurrent SLDS (rSLDS) model 29 46 47 takes an important departure from the SLDS by
leveraging the continuous latent states x, guide the diserete state transitions. Recent work [48] use this dependency but turn to the linearization of nonlinear systems with x-space fixed points as guide for the linear dynamics. In contrast, we linearize based on observed conditions (see App. A.2
Smoothly varying dynamical systems models Switching LDS models can be contrasted with
models that smoothly interpolate between dynamical parameter regimes. The simplest example of this would be time-varying linear models (c.g. (22)): CLDS models are generalization of this idea that comes with several advantages (see $2.2). Recent work [49] relaxed switching dynamics to be subject to an approximately continuous and latent time warping factor. More similar to CLDS models
are smoothly varying latent linear dynamical models, such as the gpSLDS [19] which relax the
discrete switching in rSLDS models to allow smoothly varying soft mixtures of linear dynamics, and
the dLDS [50]. which learns a dictionary of linear dynamics combined with time-varying coefficients.
Again, the primary difference is that CLDS models achieve a similar effect but utilize observed experimental covariates to infer these dynamical transitions. This approach is closer to previous work [51 using neural network architectures to parametrize input-dependent linear recurrent dynamics.
5
Conclusion
In this work, we extend classical linear- Gaussian state space models of neural dynamics. Our results suggest that these models competitive with modern methods when the dynamical parameters vary smoothly as a function of available covariates. Like classical linear models, CL.DS models are easy to fit and interpret. Our main technical contribution was to introduce approximate GP prior over system parameters and show that this leads to closed-form inference under a Gaussian noise model.
While these results are promising, CLDS models have limitations that merit consideration. To start, we did not extend to non-Gaussian observations to model spike counts directly, nor leveraged the prior structure on the coefficients F to convey more complete parameter posteriors -w discuss both of these options in the text, and see them as essential avenues for future work. Furthermore, as their name implies, CLDS models assume conditionally linear latent dynamics. We detail in Appendix AA.2hor how our approach can be viewed as a first-order (in x) approximation to a full nonlinear system of the form
E[x++1 x,, u.] = (x,, u.), Ely, x,, u.l = fy, (x,, (x,,u,)
We derive upper bounds for the approximation error between our linearized dynamics and this nonlinear system, focusing of fa and showing that the quality of this approximation is guaranteed to be small if fa is well-behaved and the conditional covariance in x, given is small. These same guarantees help shed light onto limitations. First, these models rely on observing a time series of experimental covariates U and performance would suffer if the covariates are corrupted, such as during forecasting or with partial observations. Second, the model assumes linear dynamics conditioned on uy We believe this is good approximation in many settings of interest, particularly when there is strong tuning in x to sensory or behavioral variables u, and expect CLDS models to struggle when they are loosely correlated (e.g. cognitive tasks with long periods of internal deliberation). In these scenarios, we expect that more nonlinear approaches will outperform CLDS models when given access to large amounts of data. Nevertheless, neural recordings are often trial- limited in practice [52]. We therefore view CLDS models as a broadly applicable modeling tool for many neuroscience applications.
Future work could extend CLDS models to overcome these limitations, such as handling partially observed covariates, u. Since CLDS models can be viewed as a dynamical extension of Wishart process models (see $4 future work could also apply this method to infer across-time noise corre- lations [reviewed in [53] in addition to classical across- trial noise correlations. Recent work [54] shows how across-time correlations can be used to quantify similarity in dynamical systems- a topic that has recently attracted strong interest (55). CLDS models are potentially attractive framework for tackling the unsolved challenge of estimating this high-dimensional correlation structure in rrill-limited regimes.
10

---

Acknowledgments
VG was supported by the Porter Ogden Jacobus Fellowship at Princeton University and by doctoral scholarships from the Natural Sciences and Engineering Research Council of Canada (NSERC) and the Fonds de recherche du Québec Nature et technologies (FRQNT). JWP was supported by grants from the Simons Collaboration on the Global Brain (SCGB AWD543027), the NIH BRAIN initiative (9R0IDA056404-04) an NIH RO1 (NIH 1R01EY033064), and U19 NIH-NINDS BRAIN Initiative Award (U19NS AHW was supported by the NIH BRAIN initiative (IRFI MH133778))
References
[1] Liam Paninski, Yashar Ahmadian, Daniel Gil Ferreira, Shinsuke Koyama, Kamiar Rahnama Rad, Michael Vidne, Joshua Vogelstein, and Wei Wu. A new look state-space models for neural data. Journal of computational neuroscience, 29: 107-126, 2010.
[2] Lea Duncker and Maneesh Sahani. Dynamics on the manifold: Identifying computational dynamical activity from neural population recordings. Current opinion in neurobiology, 70: 170, 2021.
[3] Daniel Durstewitz, Georgia Koppe, and Max Ingo Thurm. Reconstructing computational system dynamics from neural data with recurrent neural networks. Nature Reviews Neuroscience, 24(1 1):69--710, Nov 2023.
[4] Zoubin Ghahramani and Geoffrey Hinton. Parameter estimation for linear dynamical systems. Technical report, University Toronto, 1996. Tech. Rep. CRG-TR-96-2. Available as fhttps
7 7/mlg. eng. cam ac uk//oobbinccuurs04/rrr 9.-996622200 pdf
[5] Thomas Kailath. Linear systems, volume 156. Prentice-Hall Englewood Cliffs, NJ, 1980.
[6] H Sebastian Seung. How the brain keeps the eyes still. Proceedings of the National Academy of Sciences, 93323):1333--1334, 1996.
[7] Mark Goldman. Memory without feedback in a neural network. Neuron, 61(4):621-634 2009.
[8] Brendan K Murphy and Kenneth D Miller. Balanced amplification: selective amplification neural activity patterns. Neuron, 61(4):635 648, 2009.
mechanism of
[9] Chethan Pandarinath, Daniel J. O' Shea, Jasmine Collins, Rafal Jozefowicz, Sergey D. Stavisky, Jonathan C. Kao, Eric M. Trautmann, Matthew T. Kaufman, Stephen I. Ryu, Leigh R. Hochberg, Jaimie M. Henderson, Krishna V. Shenoy, L. F. Abbott, and David Sussillo. Inferring single-trial neural population dynamics using sequential auto-encoders. Nature Methods, 15(50):005-815 September 2018.
[10] Joel Ye and Chethan Pandarinath. Representation learning for neural population activity with neural data transformers. arXiv preprini arXiv:2108.01210. 2021.
[11] Mehdi Azabou, Vinam Arora, Venkataramana Ganesh, Ximeng Mao, Santosh B Nachimuthu, Michael Jacob Mendelson, Blake Aaron Richards, Matthew G Perich, Guillaume Lajoie, and Eva L Dyer. A unified, scalable framework for neural population decoding. Thirty-seventh Conference on Neural Information Processing Systems, 2023.
[12] Jaivardhan Kapoor, Auguste Schulz, Julius Vetter, Felix C Pei, Richard Gao, and Jakob H. Macke. Latent diffusion for neural spiking data. In The Thirty-eighth Anmual Conference on Neural Information Processing Systems, 2024.
[13] Florian Hess, Zahra Monfared, Manuel Brenner, and Daniel Durstewitz. Generalized teacher forcing for learning chaotic dynamics. In Proceedings of the 40th International Conference on Machine Learning, 2023.
[14] Matthijs Pals. A Erdem Sağtekin, Felix C Pei, Manuel Gloeckler, and Jakob H. Macke. Inferring stochastic low-rank recurrent neural networks from neural data. In The Thirty-eighth Annual Conference Neural Information Processing Systems, 2024.

---

[15] Amin Nejatbakhsh, Isabel Garon, and Alex Williams. Estimating noise correlations across continuous conditions with wishart processes. In Advances in Neural Information Processing Systems, 2023.
[16] Byron M. Yu, John P. Cunningham, Gopal Santhanam, Stephen I. Ryu, Krishna V. Shenoy, and Maneesh Sahani. Gaussian-process factor analysis for low-dimensional single-trial analysis of neural population activity. Journal of Neuruphysiology, 102(1):614466) July 2009.
[17] Biljana Petreska, Byron M Yu, John P Cunningham, Gopal Santhanam, Stephen Ryu, Krishna V Shenoy, and Maneesh Sahani. Dynamical segmentation of single trials from population neural data. In Advances in Neural Information Processing Systems, 2011.
[18) Josue Nassar, Scott Linderman, Monica Bugallo, and II Memming Park. Tree-structured recur- rent switching linear dynamical systems for multi-scale modeling. In International Conference on Learning Representations, 2019.
[19] Amber Hu, David Zoltowski, Aditya Nair, David Anderson, Lea Duncker, and Scott Linderman. Modeling latent neural dynamics with gaussian process switching linear dynamical systems, 2024.
[20] James Hensman, Nicolas Durrande, and Amn Solin. Variational fourier features for gaussian processes. Journal of Machine Learning Research, 18(151):1-52, 2018.
[21] Viacheslav Borovitskiy, Alexander Terenin, Peter Mostowsky, and Marc Deisenroth. Matérn gaussian processes riemannian manifolds. In Advances in Neural Information Processing Systems, volume 33, 2020.
[22] Antonio C Costa, Tosif Ahamed, and Greg Stephens. Adaptive, locally linear models of complex dynamics. Proceedings of the National Academy of Sciences, 11665):1011-1100 2019.
[23] Aniruddh R Galgali, Maneesh Sahani, and Valerio Mante. Residual dynamics resolves recurrent contributions neural computation. Nature Neuroscience, 26(2):326-338, 2023. 20233
[24] A. P. Dempster, N. M. Laird, and D. B. Rubin. Maximum likelihood from incomplete data via the em algorithm. Journal of the Royal Statistical Society Series B: Statistical Methodology, 1-222,, September 1977.
[25] Alex H Williams, Tony Hyun Kim, Forea Wang, Saurabh Vyas, Stephen I Ryu, Krishna V Shenoy, Mark Schnitzer, Tamara G Kolda, and Surya Ganguli. Unsupervised discovery of demixed, low-dimensional neural dynamics across multiple timescales through tensor compo- nent analysis. Neuron, 98(6): 099-11191 2018.
[26] Jakob H Macke, Lars Buesing, John P Cunningham, Byron M Yu, Krishna V Shenoy, and Maneesh Sahani. Empirical models of spiking in neural populations. Advances neural information processing systems, 24, 2011.
[27] Ian H Stevenson. Flexible models for spike count data with both over-and under-dispersion. Journal of computational neuroscience, 41:29-43, 2016.
[28] Felix Pei, Joel Ye, David M. Zoltowski, Anqi Wu, Raced H. Chowdhury, Hansem Sohn, Joseph O'Doherty, Krishna V. Shenoy, Matthew T. Kaufman, Mark Churchland. Mehrdad Jazayeri, Lee E. Miller, Jonathan Pillow, Memming Park. Eva L. Dyer, and Chethan Pandarinath. Neural latents benchmark "21: Evaluating latent variable models of neural population activity. In Advances in Neural Information Processing Systems (NeurIPS), Track on Datasets and Benchmarks, 2021.
[29] Scott Linderman, Matthew Johnson, Andrew Miller, Ryan Adams, David Blei, and Liam Paninski. Bayesian Learning and Inference Recurrent Switching Linear Dynamical Systems. In Proceedings of the 20th International Conference on Artificial Intelligence and Statistics, 2017.
[30] Brad K. Hulse and Vivek Jayaraman. Mechanisms underlying the neural computation of head direction. Annual Review of Neuroscience, 43(1):31-54, July 2020.
12

---

[31] K Zhang. Representation of spatial orientation by the intrinsic dynamics of the head-direction cell ensemble: theory. The Journal of Neuroscience, 16(6):21122122) March 1996.
[32]
Adrien Peyrache, Marie M Lacroix, Peter Petersen, and György Buzsáki. Internally organized mechanisms of the head direction sense. Nature Neuroscience, 18(4):569 -5756,9955,,,0 March 2015.
[33] N. Even-Chen, B. Sheffer, Saurabh Vyas, Stephen I. Ryu, and Krishna V. Shenoy. Structure and variability of delay activity in premotor cortex. PLOS Computational Biology, 150))00000808 February 2019.
[34] David Sussillo and Omri Barak. Opening the black box: low-dimensional dynamics in high- dimensional recurrent neural networks. Neural computation, 25(3):626-649, 2013.
[35] Andrew Gordon Wilson and Zoubin Ghahramani. Generalised wishart processes. In Proceedings of the Twenty-Seventh Conference on Uncertainty in Artificial Intelligence, pages 736-744, 2011.
[36] Neil D. Lawrence. Learning for larger datasets with the gaussian process latent variable model. In Proceedings of the Eleventh International Conference on Artificial Intelligence and Statistics, volume of Proceedings of Machine Learning Research, pages 243-250. PMLR, 2007.
[37] Jack Wang, Aaron Hertzmann, and David Fleet. Gaussian process dynamical models. In Advances in Neural Information Processing Systems, volume 18. MIT Press, 2005.
[38] Richard Turner and Maneesh Sahani. A maximum-likelihood interpretation for slow feature analysis. Neural Computation, 19(4): 1022-1038, April 2007.
[39] Matthew Dowling, Piotr Sokót, and Il Memming Park. Hida-mat', 'ern kernel. arXiv preprint arXiv:2107.07098. 2021.
[40] Matthew Dowling, Yuan Zhao, and Il Memming Park. Linear time GPs for inferring latent trajectories from neural spike trains. In Proceedings of the 40th International Conference on Machine Learning, 2023.
[41] Jouni Hartikainen and Simo Sarkka. Kalman filtering and smoothing solutions temporal gaussian process regression models. In 2010 IEEE International Workshop on Machine Learning for Signal Processing, page 379-384. IEEE, August 2010.
[42] Weihan Li, Chengrui Li, Yule Wang, and Anqi Wu. Multi-region markovian gaussian process: an efficient method to discover directional communications across multiple brain regions. In P Proceedings of the 41st International Conference on Machine Learning, ICML'24. JMLR.org, 2024.
[43] Kevin Murphy. Learning switching kalman filter models. Tech Report 98-10, Compaq Cam- bridge Research Lab, Cambridge, MA, 1998.
[44] Vladimir Pavlovic, James M Rehg. and John MacCormick. Learning switching linear models of human motion. In T. Leen, T. Dietterich, and V. Tresp, editors, Advances in Neural Information Processing Systems, volume 13. MIT Press, 2000.
[45] Victor Geadah, International Brain Laboratory, and Jonathan W. Pillow. Parsing neural dynam- ics with infinite recurrent switching linear dynamical systems. In The Twelfth International Conference Learning Representations, 2024.
[46] David Zoltowski, Jonathan Pillow, and Scott Linderman. A general recurrent state space framework for modeling neural dynamics during decision-making, In Proceedings of the 37th International Conference Machine Learning, 2020.
[47] Manuel Brenner, Christoph Jürgen Hemmer, Zahra Monfared, and Daniel Durstewitz. Almost- linear RNNs yield highly interpretable symbolic codes in dynamical systems reconstruction. In The eiirth- Annual Conference on Neural Information Processing Systems, 2024.
[48] Jimmy Smith, Scott Linderman, and David Sussillo. Reverse engineering recurrent neural networks with jacobian switching linear dynamical systems. In Advances in Neural Information Processing Systems, 1021.
13

---

[49] Julia Costacurta, Lea Duncker, Blue Sheffer, Winthrop Gillis, Caleb Weinreb, Jeffrey Markowitz. Sandeep R Datta, Alex Williams, and Scott Linderman. Distinguishing discrete and continuous behavioral variability using warped autoregressive hmms. Advances in neural information processing systems, 015223838-388.5 2022.
[50] Noga Mudrik, Yenho Chen, Eva Yezerets, Christopher J. Rozell, and Adam S. Charles. Decom- posed linear dynamical systems (dlds) for learning the latent components of neural dynamics. Journal of Machine Learning Research, 25(59):1-44) 2024.
[51] Jakob N. Foerster, Justin Gilmer, Jascha Sohl-Dickstein. Jan Chorowski, and David Sussillo. Input switched affine networks: An RNN architecture designed for interpretability. In Doina Precup and Yee Whye Teh, editors, Proceedings of the 34th Internationa Conference on Machine Learning, volume 70 of Proceedings of Machine Learning Research, pages 1136-1145. PMLR, 06-11 Aug 2017.
[52] Alex H Williams and Scott W Linderman. Statistical neuroscience in the single trial limit. Current Opinion in Neurobiology, 70:193-205, 2021.
[53] Stefano Panzeri, Monica Moroni, Houman Safaai, and Christopher D Harvey. The structures and functions of correlations in neural population codes. Nature Reviews Neuroscience, 23(9):551-
567, 2022.
[54] Amin Nejatbakhsh, Victor Geadah, Alex H Williams, and David Lipshutz. Comparing noisy neural population dynamics using optimal transport distances. In The Thirteenth International Conference on Learning Representations, 2025.
[55] Mitchell Ostrow, Adam Eisen, Leo Kozachkov, and Ila Fiete. Beyond geometry: Comparing the temporal structure of computation in neural circuits with dynamical similarity analysis. In Advances in Neural Information Processing Systems, 2023.
(56) Philip Greengard, Manas Rachh, and Alex H. Barnett. Equispaced fourier representations for efficient gaussian process regression from a billion data points. SIAM/ASA Journal on Uncertainty Quantification, 13(1))63-89, January 2025.
14

---

A Modeling
A.1 Inference
A.1.1 Basis space form
We want to show
M(u) W W (o(u) 8 X)
for M defined as in 1c with basis functions (or) U, and X E RDxxss D2 ENN
For each j-t entry,
[M
M(u)a [x]
D2
We (1b) (2))
(Ø(u)) & Dsetks We (1k)
D2
Σ
Wi.D.et& (Ø(u) @ D) k=
W., (Ø(u) @ X
[ (Ø(u) X)
as desired, where we' 've defined W D.e+k, = W
(ik)
A.1.2 Least squares derivation for Conditionally Linear Regression Recall
M(u) = WT(((u) & X), W E D2LXD
for M(u) E D2 X E IR D2xM In particular, M(u,)x, = W 2,, What follows are standard
least-squares derivations for matrix coefficients with matrix regularization, which we include for completeness.
Our posterior objective reads as
log W YI:N, XI.N, uI.N) x logp(yi: W, Xi:N, uINN) + log p( W)
logp(yn W,x,,
Xn, un) + logr(
We have
logp(ya
Xn, ue)
(yn W 2n) TI- (yn WT 2n) -
- (y, W 2n) TD- - W kn)] -C
1 1Tr 2--(yy W zn)(yn W 2n) -C
1 (Tr Σ yny, 2Tr Σ W z,y
+ + (-T,,, W|)
with the normalizing constant 2 log 12
To optimize this expression with respect W, consider the zeros of the derivative
log p(yn
Xn, u.)
Tr z,y, ,,,, 2aw
Tr - 'W
2 2 W]
aw
aw
Σ 2 WW-1
Z, n
15

---

and
log W) = W
SW log (W) W
Taken together, we thus have that the stationary point of the posterior satisfies
(nNY nn wE- W 0.
Define Y E RNxD; Z E NNDxL by row-wise stacking and n, respectively, and note that En zn Z Y. We get
ZWE- Z'YE-1 -Z'Z---- - W =
zwzw + WE y
as desired.
A.1.3 Joint dynamics and bias EM update in weight-space
Here we detail how one EM update for the parameters governing {A, b) is carried out. Using the iunciionspacce weights W A E DDLD D and Wi IR LxD,, the dynamics read
X++1
A b(u)) + Et
W. (Ø(u)) & x1) + W o(u)) Ec
Define 2, = x((1) 2 X, E LLD. and note
Ex, [2e] = (w+)) 3 Ex, [1r]]
8x,x++1
[x,T11] x+11 ( & Exs [xxx1+1]
The quantity of interest for the M-step from the complete data log-likelihood is (up to an additive constant)
log p(x+111 X4, F, u))
E
[1+1 +
+ b(u)))) Q-1 (x++1 (A(u))x + b(u:))) + const.
E
2 x141 2x A(u,)x: -
+ x T A A A
b(uu) + const.
-
E E -Q A(u, 11 2Q bb
( ( b(u bb(u))
const.
Which with the finite basis expansion reads as
- Tr
Q E Ex - ((kkk) @ - (TT
2
1=1
-W W
E[x]]]
16

---

The partial derivatives satisfy, denoting Φ: = =(u)),
OC
Σ ( @ F(xxx - ( 2 [ ( 2 E[x,)] W.Q-
OWA
NAQ
- W
(14)
0C
Σ
o,E ( 8 E W QQ + 0.0, W WoQ-1
OWb
==1 = =x - ZW
(15)
where we've defined the matrices + ((T Z ( X (T obtained by
stacking ((z)), ((u) & E[x,] and [[x++1] respectively for t fI....,T 1}, and the sufficient statistics
F. & ][xxx]] NA Σ s(u)) @ E[xxx,+11+
(2)
which are all defined during our E-step. These statistics are contrast the "typical" sufficient stats without the weight-space parametrization
T Σ
M(
Σ E
MA
E(x,x
Incorporating the Gaussian prior on WA and WT and equating the two partial derivatives in obtain the stationary points, we obtain the system of equations
[4 15
A W
Q
+
N 7
T + T +
WA
A W
NA
(16)
-
X
which solving for WA A and W. jointly amounts to solving our M-step.
A.2 Nonlinear dynamics: linearization and composite dynamics
Consider input-driven nonlinear dynamics in x, E R.,,
X++1 f(x;, Us) +
(17)
governed by f RD x 4 + RDD, and where 6 is zero-mean Gaussian noise. We assume that f has continuous and bounded second-order partial derivatives in both x and u. Below we treat the state and input variables (x, u.) as random variables that are jointly drawn from some unspecified distribution.
A.2.1 CLDS conditional approximation error
Let fID denote the output dimensions of f: that is, f(x,, us)
[f (xe, u.) fD(x,, u.)] In a first step to relate CLDS dynamics to f. consider the first-order Taylor expansion for each output dimension in the first argument æ about E RD, For i DJ. this
fi(,, u²) fi(a, u.) Vi,i T a)
(18)
where V denotes the vector-valued gradient of with respect to it's first argument and is the residual of the Taylor approximation. The Lagrange remainder form of Taylor's theorem tells us that this residual can be expressed as:
& = TV2,a, ux)xx -a)
(19)
for some E RR, where V2.f, is the matrix-values Hessian of fa with respect to its first argument. We can upper bound the absolute value of this remainder using the Cauchy-Schward Schwartz inequality and standard operator norm inequality. Specifically, for E (1...., D), we have
| < f.lC, 41)l ||x; al/
(20)
17

---

where V2 u2)ll denotes the maximal singular value (operator norm) of the matrix Vifi u,) R DxD We assume that this operator norm is upper bounded globally by constant 7i>0 0,
[ u) S. L.
Vx,uER".
(21)
Intuitively, this assumption implies that the second-order derivatives of f with respect to are not too large, meaning that the accuracy of the first-order Taylor approximation degrades in proportion to the magnitude of curvature in f.
Returning equation (20), proceed by taking conditional expectations with respect x given u, on both sides of the inequality. This yields an upper bound on the expected approximation error,
E | u,] - L, E ||x,
a
2
This upper bound is minimized by choosing a [x Uz]- Plugging in this choice, observe that
E [[
u,] - T, Tr [oo[xx un]]
where Cov uz] = E [(x: - E[x; uz]) E[x; u.)) uu.] is the conditional covariance of X₁ given u, Finally, we can sum these upper bounds over = {I...., D} to bound the total expected approximation error
Σ EB[E ur] S L. Tr [Cov[x:
u.]]
(22)
where we have defined L= E; Li as a global constant bounding the second-order smoothness of f across all dimensions.
Returning to equation (8) and plugging in the optimal choice of a E(x; Uz). we obtain the following CLDS approximation to the nonlinear dynamics
h(x;, u+) f (IE(x u.]. u:) + Ve ff ue], ux)(x; E|xr uz]) A(u,)x + b(u))
(23)
where V is the matrix-valued Jacobian, fx u) E RDD, with respect to the first argument of f and we have re-arranged the terms and detined
A(u:) = vf
ur], ue), b(u,) = f[E[x; u.J, u,) - V v uz]. un 3[[x uz).
For each value of u, the quality of this approximation is guaranteed by equation (22 to be small if the second-order derivatives of f with respect to are small and if the conditional variance of x, given u, is small. We note that this analysis of approximation error and does not account for the additional estimation eтror we incur when learning the functions A(u) and b(u) from noisy and limited data. Nonetheless, this analysis tells us that we expect the CLDDS to perform well in circumstances where the underlying dynamics are smooth and the conditional distributions of x given u, have low variance.
A.2.2 Composite dynamics
When considering the composite dynamics in 7 (Section $3.1 we are interested in approximating the input-driven nonlinear dynamics given by equation with autonomous nonlinear dynamies, governed by some function g such that
f(x, u:) g(x:).
We can evaluate the quality of this approximation using the conditional expectation of the squared error
E [Ilf(x, u.) - 9 x
where the conditional expectation taken over u, given Xx. This approximation error is minimized by choosing
(xx)) (x, u.) x,]
E[x xx]
However, learning f from limited data is challenging, so we replace this with our CLDS model to achieve the composite dynamical system
g (xx) R EE [xx, u:) xs]
(24)
18

---

where h(xx u) = A(uxx + b(u) as in equation (23) Now we analyze the quality of this approxima- tion. Consider the residual along dimension [T.... D],
R₂ fx,,, u.) Eii (x,, ū.) xi]
where the expectation in the second term taken over ūg which is drawn from the conditional distribution of u, given x,. That is. R2 is a random variable that depends on joint sample of
(x,, u,) from the stationary distribution, and ü, is dummy variable that is integrated during the calculation of R, To proceed, we note that
R₂ = (x,, uz) - h4 (x,, ū.) x+]
[x uy) fo (xx, ū.) fx üi) hy (x,, ūi)) xt]
and, applying Jensen's inequality and the triangle inequality, we conclude
[Ra S Ea,xx. [xx,, uz) - fi( ü.) + fi(e; ūt) h1s (x,, ūt][]
< Eu, [xx,, u,) - [i( út) + Ea. (x,, ü.) - he(xe, ūx)]]
where we have introduced a minor change in notation, using Ea,ix, to denote the conditional expectation of ü, given X,. Now we take expectations on both sides of this inequality with respect to the remaining random variables, x, and u, which are sampled from some stationary distribution associated with the dynamical system. Using the law of total expectation, we obtain
E R|R.| < Ex; G. [ (x,, ue) fi(X,, úz] + Eu, [Fxx/uu [[ us) h4 (x+,,, ue)
(25)
On the right hand side, the first term takes the conditional expectation over u, and ü, followed by an expectation over xt. The second term reverses this order, taking the conditional expectation over x2, followed by an expectation over u, (since this is identically distributed to ú, we drop the tilde).
To upper bound the first term, we introduce an assumption that Sa is Lipschitz continuous in its second argument. That is, there exists a constant C. > such that
| u) fi(x,, u')l < C:llu u'1l2 In conjunction with Jensen's 's inequality, this Lipschitz assumption implies the following upper bound:
Yu, u' E Eu.
E, [x u,) fx úi) 1 C. Em, ax [ - u
1 CiV En 1x u, ūن l
= =C.v Tr Tr//2//Covv/uu x
remains to upper bound the second term in equation 25 A direct application of the results in 1.22.1 yields the bound
Ex, un [i u.) - h, hx uz)l < L. Tr|Cov[x, ur]
where Li, defined in (21) is constant bounding the second derivatives of fa with respect to X Putting these pieces together we conclude that
E|R C Ex, / TT x + La Ru, Tr(Cov[x; u.]]
And so an upper bound on the total absolute error of the composite dynamics is given by
Σ F.R. S C Ex, v 2TH((Covuu xx] + T. EE,, Tr[Covix, uz]]
where we have defined C = Σ; C and I. E.
In summary, we have shown that the approximation error of the composite dynamical system, defined in equation 77. is bounded by a sum two terms. The first term approaches zero the limit that the conditional covariance of u, given x, goes to zero, while the second term approaches zero in the limit that the conditional covariance of X; given u, goes to zero. Thus, the composite dynamics have the potential to provide accurate depiction of the true nonlinear dynamical system if x, and u, close to being in one-to-one correspondence with each other.
We note that in practice when computing the composite dynamics in 7 we make the simplifying assumption that p( [1 X, = æ²) does not depend on (i.e. we assume that this decoding distribution is stationary).
19

---

A.3 Correspondence between CLDS and Switching LDS
In this section, we explore the relationship between the linear time-variami dynamics of lap with Aij sid (O kz), and the dynamics of switching linear dynamical system. While the latter has parameters evolving over discrete set, we can nonetheless explore how to think of this discrete support embedded within R PxP, and seek to match the moments of these processes.
The first thing to note that by drawing the entries of Aij i.id.., we can gain insight by considering single process a(.) ~ GP(O,, k),, and a SL.DS with one-dimensional dynamics. Thus, we consider the SLDS model
p 74 = j) Pig
(26a)
X++1 (4) X - Et
(26b)
with discrete states z, governing dynamies in X2, with transition matrix P and dynamics a for R|2l for the vector of the |2| K values that can be taken by the stochastic process and similarly a (a(=1), ..... (=*x) T_ We consider the stationary distribution for the z process. Finally, note that the map 2n + (=*)) is deterministic, one-to-one and Let a == a(*), and
E Z. Denote 2 (Z1.....
onto, such that we think of the Markov chain in 74 having support denote a; (=)).
Let us now explore the moments the stochastic process a, determine its relationship with the CLDS. Assume has reached stationarity, such that p(z) = p(ar) = π. Then, first,
aa
(27)
and then we have the cross-correlation
E
[a
E
[E
a ac]]
K
E Jara +n a. aalp(aa = =a,)
ag (p(ar-n
a a,)) p(a: a,)
K K
a diag((a)
Denote II = diag(T).. We get the desired covariance
Cov(ar, ar4n) = E [ara++]] - EE [ar] E [actn] = =a" (["P11 nπ
yielding our kernel form
Cor(s(t)).a(tt)) ke(t,,t)) a (plii-filII- कक
(28)
Hence, in all, our approximation of the stochastic process over IR up to the first two moments is
a(-) GP a. (pit--tl11 कग
(29)
which in particular can be made zero-mean by consider values a such that म TT a = 0. This establishes the form of the GP prior over a(-) for which the corresponding CLDS best matches the SLDS.
B Experiments
B.1 Model implementations
We initialize observation matrices for all models (C in LDS, CLDS and gpSLDS models, log-rate decoder weights LFADS) as the PCA principal axes in y-space- that is, the top / right singular vectors of the data-for for each task.
20

---

CLDS In all experiments, we assume that d(uu) O, which forces the predicted firing rates, conditioned on uy lie in N-dimensional space spanned by the columns of C(u)).
gpSLDS We implemented the gpSLDS [9] using its accompanying code at httts:///ithhti com htam/lindermannnn gpslds We used latent dimension of D = in all experiments as encountered significant OOM errors with higher dimensions. We used exactly the same inputs U1:r those provided to the (C)LDS models. We J and the feature transformation
((x) =
21 222
22
(30)
as to allow both linear (as the rSLDS) and circular म (2) partitions. We experimented with providing the same basis functions for the feature transformation the ones provided to the CLDS but found that did not provide any advantage. Other pper-ppraeeters we set to default provided values.
LFADS We use the Jax implementation of LFADS available https://gittht: com/ googlerreeearee resalrrrrrraccc/ooottttttttt computation. thru- caaru-dynamics/// bllyyaamics//rreemmatteer///////
tutorial, which is under Apache-2.0 license. We choose the factor dimension be the same as the latent dimension D of the (C)LDS models and the inferred inputs to be of dimension [2/], both following the (C)LDS models on any given experiment. The other components of the architecture are held fixed across all experiments:
Encoder, controller and generator have hidden-states of dimension 32;
The inferred inputs are modeled as having an auto-correlation of 1.0 and a variance of 0.1;
We train for 1000 epochs, with an initial learning rate of 0.5 with exponential decay rate 0.995, along with a KI. warm-up coming in at 500 steps.
B.2 Synthetic experiment and parameter recovery
True
Tuning width = 0.1
Tuning width
001
104770
JC+02
10
Je-04
10440
001
Tuning width 0.5
-3++01
-1e-04
30420
02
01
-01
Figure 5: Recovery of C. Rows indicate varying level of tuning curve width y for (u). Recovery becomes more challeng- ing for smaller width since it requires a higher and higher number of bases L to approximate the true tuning bump.
Figure 6: Hyperparameter search, CLDS marginal log-likelihood on held-out vali- dation dataset set for E {0.2, 0.6, 1.0] and E {0.2, 0.5, 1.0, 1.5, 2.0}. Maximum at- tained (R, o. {0.6, 0.5).
For the peaks& & spanning regularly the interval of C for the synthetic experiment are defined as
[-T, m) and widths y,, the tuning curves the rows
(u) = {
=16 cos
u- 7
u
if E (E: - y y
(31)
else
21

---

We plot the its recovery with our inference procedure in Fig. 5 up to the invertible transform non-identifiability inherent to LDS models.
In Table T we report the results of our inference method on the first synthetic ring-attractor experiment (..2) for varying emission noise scale R of J, keeping all other experimental setup parameters exactly the same as in the original experiment and fixed. We see that even as the signal-to-noise ratio degrades and we reach a co-smoothing value of 0.21, which is lower than our fit on the monkey data we still have decent A recovery and near-perfect Rt recovery.
of $3.
True noise log scale log on
-2
-1
0
1
Recovered R log scale
-1.97
-0.98
0.02
1.02
A-) recovery eπror
0.01
0.02
0.11
0.32
Co-smoothing R²
0.99
0.94
0.68
0.21
Table 1: Impact of observational noise on quality fit. The "recovered R scale" is the square-root of the matrix 2-norm of R. The (-) recovery error the average in 2-norm error (between recovered and true) in eigenvalues over regular grid of 50 angular conditions A. The co-smoothing R² is the average validation set R² on top-5- vor-i-variance neurons, like in the main text.
B.3 Neural circuit ring attractor experiment
Let y(o,t) t) be the firing-rate "bump" on angles & E [o, 2n) with periodic boundary conditions, and let v(t) = ((t) be the angular velocity derived from the head-direction input. We consider its dynamics governed by the stochastic PDE
dy(ó, 1) = [-v((, + ((w ()) + k,,(() 0su(d, t)] di + adB())
with
f(.) pointwise ReLU nonlinearity.
The recurrent drive a circular convolution
(w (b) L 2x
w(s 4) ()),
with a rotationally symmetric kernel (in radians)
W(A)= Je exp (-)/A/))))) Jj exp (A)/)))))
for A E (- T] (wrapped), Jeh, {J4, ) parameters of excitation and inhibition.
The term Ky v(t) døy shifts the bump at speed proportional to v(t), of parameter Ke.
cylindrical Wiener process on the ring with covariance ((B(o, B(d,t)) t')] = ô(Φ o') min(t, E') ("white-noise" dB/ dB)dt).
0 the time constant.
Numerical implementation We discretize time to Al and head-direction to ti E (1,...,NY,
N = 32, regularly spanning [0, 2Tr). Let y(d,, nAt) [yn]-- Consider W circulant matrix implementing w, and D a centered difference discretization of the derivative. Our update is
yr = = + ((Wy.) + Ku 0, Dyt y++1 y, yr + At + +o V n 1(0,1
as the Euler-Mant Marryyam discretization of the SPDE.
After search over kernel hyperparameters, we found the best-performing model (Fig. 7 to encode ring of fixed points that closely resembles the synthetic HD fixed points (3.2)). and with eigenvalues closer the ones observed from the CLDS fits on the mouse ADn recordings ($3.4). The results make confident that the CLDS can capture this nonlinear structure even under model mismatch.
22

---

Food
4.1
Eqg vat. Ei vec
-10
13
2.
0.5
Figure 7: Inferred CLDS fixed points (a), local dynamics (b) and eigenvalues of A (c) function of head-direction 0 for the neural circuit ring attractor experiment. Hyperparameter values are (a, x) {1.0, 0.4), obtained held-out LL over grid search.
B.4 Data and Pre-processing
Mice Head-Direction We analyze neural recordings of antero-dorsal thalamic nucleus (ADn) from Ref. (32] of the mouse HD system in mice foraging in an open environment (Fig. 33). The data was accessed through the pynapple e (PYthon Neural Analysis Package,
MIT
License) package, with the data itself stored in NWB format on OSF at htthttsp////o 10/ jb2gd
We considered neural activity from the "wake" period, binned in 50ms fime-bins, then processed to firing rates by smoothing over a window of bins, and separated into 10s trials.
Macaque center-out reaching We analyzed neural recordings of dorsal premotor cortex (PMd) in macaques performing center-out reaching task from Ref. 33 the associated article is under CC BY 4.0 license.
The experiments supported different reach radii, and we selected only the middle reach radius at 8cm. We were left with only the angular direction the reach condition, over N = 16 possible reach angles. We aligned all trials around the go-cue, selecting 200 ms before the go-cue and 300ms after. We binned the data into 5ms bins, and performed Gaussian kernel smoothing with a standard deviation 0.5 over bins.
B.5 Hyper-parameters
Throughout all experiments, we' ve set L = 5 to balance expressivity and number of parameters. While it is possible to treat instead set I to large enough value that results in negligible error in the GP kernel approximation. As L o0, we pay a larger computational price fitting the model, but the statistical properties of the model are essentially unchanged because high-frequeng basis functions have nearly zero amplitude. This idea is well established GP regression literature [56].
hyperparameter tuned improve performance, we propose to
To select the other hyper-parametent of GP prior length-scale and scale a, we evaluated the various models held-out validation set. We plot in Fig. our search over {x, on the macaque center-out reaching data.
Finally, dimensionality of D 2 for the head-direction experiments is chosen throughout for interpretability, as the latent space thought to encode head-direction. For the macaque center-out reaching data, computed co-smoothing over D E {3,5, 10, 15) (Table 2 and selected D 5.
D
15
CLDS
0.229
0.232
0.207
0.168
LDS
0.211
0.193
0.196
0.170
Table 2: Co-smoothing R reconstruction of single held-out neurons from the test-set, over varying latent dimensionality D, averaged over two random seeds.
B.6 Compute Resources
We list below the compute resources used per experiment:
23

---

1. Synthetic HDD experiment: Results were computed on personal machine (Apple MacBook Pro, M2 Pro chip. 32G RAM). Inference converges in few seconds.
2. Mouse HD experiment: Results were computed a personal machine (Apple MacBook Pro, M2 Pro chip, 32G RAM). Inference converges in a few minutes.
3. Monkey-reaching experiment: Results were computed on an external cluster equipped with NVIDIA A100 GPUs Fitting single model typically used 60GB of GPU memory. All CLDS models were run strictly on CPU-only nodes (Intel Xeon Silver 4309Y Processor, 2.80 GHz, 8 cores/16 16 threads), using approximately 2 GB of RAM per run. Wall-clock
times were between ten minutes and one hour.
24

---

NeurIPS Paper Checklist
1. Claims
Question: Do the main claims made in the abstract and introduction accurately reflect the paper's contributions and scope?
Answer: [Yes]
Justification: The paper's main claims and contributions are provided list in the intro- duction, all with pointers to the relevant sections of the manuscript.
Guidelines:
The answer NA means that the abstract and introduction do not include the claims made in the paper.
The abstract and/or introduction should clearly state the claims made, including the contributions made the paper and important assumptions and limitations. A No or NA answer to this question will not be perceived well by the reviewers.
The claims made should match theoretical and experimental results, and reflect how much the results can be expected to generalize to other settings.
It is fine to include aspirational goals as motivation long is clear that these goals are not attained by the paper.
2. Limitations
Question:
Does the paper discuss the limitations of the work performed by the authors?
Answer: [Yes]
A discussion on the limitations can be found in the second paragraph of the conclusion.
Justification:
Guidelines:
The answer NA means that the paper has limitation while the answer No means that the paper has limitations, but those are not discussed in the paper.
The authors are encouraged to create a separate "Limitations" section in their paper.
The paper should point out any strong assumptions and how robust the results are to violations of these assumptions (c.g.. independence assumptions, noiseless settings, model sell-speciiiitiinn asymptotic approximations only holding locally). ). The authors should reflect on how these assumptions might be violated in practice and what the implications would be.
The authors should reflect on the scope of the claims made, e.g. if the approach was only tested on few datasets or with few runs. In general, empirical results often depend on implicit assumptions, which should be articulated.
The authors should reflect on the factors that influence the performance of the approach. For example, facial recognition algorithm may perform poorly when image resolution low or images are taken in low lighting. Or speech-to-tex, system might not be used reliably to provide closed captions for online lectures because it fails handle technical jargon.
The authors should discuss the computational efficiency of the proposed algorithms and how they scale with dataset size.
If applicable, the authors should discuss possible limitations of their approach to address problems of privacy and faimesss.
While the authors might fear that complete honesty about limitations might be used by reviewers as grounds for rejection, a worse outcome might be that reviewers discover limitations that aren't 't acknowledged in the paper. The authors should use their best judgment and recognize that individual actions in favor of transparency play an impor- tant role in developing norms that preserve the integrity of the community, Reviewers will be specifically instructed to not penalize honesty concerning limitations.
3. Theory assumptions and proofs
Question: For each theoretical result, does the paper provide the full set of assumptions and complete (and correct) proof?
25

---

Answer: [Yes]
Justification:
The derivations for the theoretical results of 22. can be found in Appendix $A.1] the claims on linearization are supported by Appendix $A.2 and necessary further derivations (than the main text) on model equivalences can be found in Appendix A.3
Guidelines:
The answer NA means that the paper does not include theoretical results, All the theorems, formulas, and proofs in the paper should be numbered and cross- referenced.
All assumptions should be clearly stated or referenced the statement of any theorems.
The proofs can either appear in the main paper or the supplemental material, but they appear in the supplemental material, the authors are encouraged to provide short proof sketch to provide intuition.
Inversely, any informal proof provided in the core of the paper should be complemented by formal proofs provided in appendix supplemental material.
Theorems and Lemmas that the proof relies upon should be properly referenced.
4. Experimental result reproducibility
Does the paper fully disclose all the information needed to reproduce the main ex- perimental results the paper to the extent that affects the main claims and/or conclusions of the paper (regardless of whether the code and data are provided or not)?
Question:
Answer: [Yes]
The contribution is primarily a new model, which is described in the text and example derivations for inference are fully carried out in Appendix SA.I Accompanying code is available at https: //gi thub. com/ hwurostatssssaae aatasssaaaaaa clds Finally, Appendix providees complimentary information to the text for experiments, which contains most information needed to reproduce the experimental results.
Justification:
Guidelines:
The answer NA means that the paper does not include experiments.
If the paper includes experiments, a No answer to this question will not be perceived well by the reviewers: Making the paper reproducible is important, regardless of whether the code and data are provided or not.
If the contribution is dataset and/or model, the authors should describe the steps taken to make their results reproducible or verifiable.
Depending on the contribution, reproducibility can be accomplished in various ways. For example, if the contribution is novel architecture, describing the architecture fully might suffice, or if the contribution is a specific model and empirical evaluation, it may be necessary either make it possible for others to replicate the model with the same dataset, or provide access to the model. In general. releasing code and data is often one good way accomplish this, but reproducibility can also be provided via detailed instructions for how to replicate the results, access to a hosted model (e.g.. in the case of large language model), releasing of a model checkpoint, or other means that are appropriate to the research performed.
While NeurIPS does not require releasing code, the conference does require all submis- sions to provide some reasonable avenue for reproducibility, which may depend on the nature of the contribution. For example
(a) If the contribution is primarily a new algorithm, the paper should make clear how to reproduce that algorithm. (b) the contribution is primarily a new model architecture, the paper should describe the architecture clearly and fully.
(c) If the contribution is new model (c.g., a large language model), then there should either be way to access this model for reproducing the results or a way to reproduce the model (e.g., with an open-source dataset instructions for how to construct the dataset).
(d) We recognize that reproducibility may tricky in some cases, in which case authors welcome describe the particular way they provide for reproducibility.
26

---

In the case of closed-source models, it may be that access to the model is limited in some way (e.g., registered users), but it should be possible for other researchers to have some path to reproducing or verifying the results.
5. Open access to data and code
Question: Does the paper provide open access to the data and code, with sufficient instruc- tions faithfully reproduce the main experimental results, described in supplemental material?
Answer: [Yes]
Justification: Code for the model, synthetic dataset generation and notebooks for analysis are publicly available at htttps:: htttps://githun gitbbb............ con/neurostatslabb clds Experimental data is publicly available, with pre-processing detailed in Appendix $B.4
Guidelines:
The answer NA means that paper does not include experiments requiring code.
Please see the NeurIPS code and data submission guidelines (https://nnis.cm cc/
Polccy for more details.
While encourage the release of code and data, understand that this might not be possible, so "No" an acceptable answer. Papers cannot be rejected simply for not including code, unless this is central to the contribution (e.g., for new open-source benchmark).
The instructions should contain the exact command and environment needed to run to reproduce the results. See the NeurlPS code and data submission guidelines (https
7 //nips https/cccpcbllcguuuceeeecc httpa///giiddee/eeeubbbiiiiiiiii sonddesubmissoonl Policy for more details.
The authors should provide instructions on data access and preparation, including how to access the raw data, preprocessed data, intermediate data, and generated data, etc.
The authors should provide scripts to reproduce all experimental results for the new proposed method and baselines. only subset of experiments are reproducible, they should state which ones are omitted from the script and why.
At submission time, to preserve anonymity, the authors should release anonymized versions (if applicable).
Providing as much information as possible in supplemental material (appended to the paper) recommended, but including URLs to data and code permitted.
6. Experimental setting/details
Question: Does the paper specify all the training and test details (e.g., data splits, hyper- parameters, how they were chosen, type of optimizer, etc.) necessary to understand the results?
Answer: [Yes]
Justification: Yes, see Appendix $B
Guidelines:
The answer NA means that the paper does not include experiments.
The experimental setting should be presented in the core the paper to a level of detail that is necessary to appreciate the results and make sense of them.
The full details can be provided either with the code, in appendix, as supplemental material.
7. Experiment statistical significance
Question: Does the paper report error bars suitably and correctly defined or other appropriate information about the statistical significance of the experiments?
Answer: [Yes]
Justification: We provide error bars on the co-smoothing results, the only metric used in model comparison.
Guidelines:
The answer NA means that the paper does not include experiments.
27

---

The authors should answer "Yes" if the results are accompanied by error bars, confi- dence intervals, or statistical significance tests, least for the experiments that support the main claims of the paper.
The factors of variability that the bars are capturing should be clearly stated (for example, train/test split, initialization, random drawing of some parameter, or overall run with given experimental conditions).
The method for calculating the error bars should be explained (closed form formula, call to a library function, bootstrap, etc.)
The assumptions made should be given (e.g.. Normally distributed errors).
It should be clear whether the error bar is the standard deviation or the standard error of the mean.
is OK to report 1-sigma error bars, but one should state it. The authors should preferably report a 2-sigma error bar than state that they have 96% CI, if the hypothesis of Normality of errors is not verified.
For asymmetric distributions, the authors should be careful not to show in tables or figures symmetric error bars that would yield results that are out range (c.g. negative error rates).
If error bars are reported in tables or plots, The authors should explain the text how they were calculated and reference the corresponding figures or tables in the text.
8. Experiments compute resources
For each experiment, does the paper provide sufficient information the puter resources (type of compute workers, memory, time of execution) needed to reproduce the experiments?
Question:
Answer: [Yes]
Justification: See Appendix SB.6
Guidelines:
The answer NA means that the paper does not include experiments.
The paper should indicate the type of compute workers CPU or GPU, internal cluster, or cloud provider, including relevant memory and storage.
The paper should provide the amount of compute required for each of the individual experimental runs as well estimate the total compute.
The paper should disclose whether the full research project required more compute than the experiments reported in the paper (e.g., preliminary failed experiments that didn't make into the paper).
9. Code of ethics
Question: Does the research conducted in the paper conform, in every respect, with the NeurlPS Code of Ethics
Answer: [Yes]
Justification:
The research conducted in the paper conforms with the NeurIPS Code of Ethics.
Guidelines:
The answer NA means that the authors have reviewed the NeurIPS Code of Ethics.
If the authors answer No, they should explain the special circumstances that require deviation from the Code of Ethics.
The authors should make sure to preserve anonymity (e.g., there a special consid- eration due to laws or regulations in their jurisdiction).
10. Broader impacts
Question: Does the paper discuss both potential positive societal impacts and negative societal impacts of the work performed?
Answer: [NA]
28

---

The modeling presented here aims to provide a methodology to enhance our understanding neural computation. Analysis of electrophysiological data can have long- term implications for the treatment and understanding of medical treatment and neurological disorders across different species. However, these considerations are far removed for the preliminary analyses and theoretical modeling presented, and we foresee no immediate societal consequences this work.
Justification:
Guidelines:
The answer NA means that there no societal impact of the work performed.
If the authors answer NA or No, they should explain why their work has no societal impact why the paper does not address societal impact.
Examples of negative societal impacts include potential malicious unintended uses (e.g., disinformation, generating fake profiles, surveillance), fairness considerations (e.g., deployment of technologies that could make decisions that unfairly impact specific groups), privacy considerations, and security considerations.
The conference expects that many papers will be foundational research and not tied to particular applications, let alone deployments. However, if there is a direct path to any negative applications, the authors should point it out. For example, it is legitimate to point out that an improvement in the quality of generative models could be used to generate deepfakes for disinformation. On the other hand, it is not needed to point out that a generic algorithm for optimizing neural networks could enable people to train models that generate Deepfakes faster.
The authors should consider possible harms that could arise when the technology is being used intended and functioning correctly, harms that could arise when the technology is being used as intended but gives incorrect results, and harms following from (intentional unintentional) misuse of the technology.
there are negative societal impacts, the authors could also discuss possible mitigation strategies (e.g., gated release of models, providing defenses in addition to attacks, mechanisms for monitoring misuse, mechanisms to monitor how system learns from feedback over time, improving the efficiency and accessibility of ML).
11. Safeguards
Does the paper describe safeguards that have been put in place for responsible release of data models that have high risk for misuse (e.g., pretrained language models, image generators, or scraped datasets)?
Question:
Answer: [NA]
Justification:
The paper poses no such risks.
Guidelines:
The answer NA means that the paper poses no such risks.
Released models that have high risk for misuse or dual-use should be released with necessary safeguards to allow for controlled use of the model, for example by requiring that users adhere to usage guidelines or restrictions to access the model or implementing safety filters.
Datasets that have been scraped from the Internet could pose safety risks. The authors should describe how they avoided releasing unsafe images.
We recognize that providing effective safeguards is challenging, and many papers do not require this, but encourage authors to take this into account and make best faith effort.
12. Licenses for existing assets
Are the creators or original owners of assets (e.g., code, data, models), used in the paper, properly credited and are the license and terms of use explicitly mentioned and properly respected?
Question:
Answer: [Yes]
Justification: Citations are present whenever relevant, and Appendix detaiis the specific assets and licenses related to the models and data considered.
29

---

Guidelines:
The answer NA means that the paper does not use existing assets.
The authors should cite the original paper that produced the code package or dataset.
The authors should state which version of the asset is used and, if possible, include URL.
The name of the license (e.g., CC-BY 4.0) should be included for each asset. For scraped data from a particular source (e.g., website), the copyright and terms of service of that source should be provided.
If assets released, the license, copyright information, and terms of in the package should be provided. For popular datasets, papppssipphhaaoocccappppp.app.pt thcode. com/ datasets has curated licenses for some datasets. Their licensing guide can help determine the license of dataset.
For existing datasets that are re-packaged, both the original license and the license of the derived asset (if it has changed) should be provided.
If this information is not available online, the authors are encouraged to reach out to the asset's s creators.
13. New assets
Question: Are new assets introduced in the paper well documented and is the documentation provided alongside the assets?
Answer: [Yes]
Justification:
The manuscript presents the model under a CC BY 4.0 license.
Guidelines:
The answer NA means that the paper does not release new assets.
Researchers should communicate the details of the dataset/code/model as part of their submissions via structured templates. This includes details about training, license, limitations, etc.
The paper should discuss whether and how consent obtained from people whose asset is used.
At submission time, remember to anonymize your assets (if applicable). You can either create an anonymized URL or include an anonymized zip file.
14. Crowdsourcing and research with human subjects
Question: For crowdsourcing experiments and research with human subjects, does the paper include the full text of instructions given to participants and screenshots, if applicable, as well as details about compensation (if any)?
Answer: [NA]
Justification:
The paper does not involve crowdsoureing nor research with human subjects.
Guidelines:
The answer NA means that the paper does not involve crowdsourcing nor research with human subjects.
Including this information the supplemental material is fine, but if the main contribu- tion of the paper involves human subjects, then as much detail as possible should be included the main paper.
According to the NeurIPS Code of Ethics, workers involved in data collection, curation, other labor should be paid at least the minimum wage in the country of the data collector.
15. Institutional review board (IRB) approvals or equivalent for research with human subjects
Question: Does the paper describe potential risks incurred by study participants, whether such risks were disclosed to the subjects, and whether Institutional Review Board (IRB) approvals (or an equivalent approval/review based on the requirements of your country or institution) were obtained?
Answer: [NA]
30

---

Justification:
The paper does not involve crowdsoureing nor research with human subjects.
Guidelines:
The answer NA means that the paper does not involve crowdsoureing nor research with
human subjects.
Depending the country in which research is conducted, IRB approval (or equivalent) may be required for any human subjects research. If you obtained IRB approval, you should clearly state this in the paper.
recognize that the procedures for this may vary significantly between institutions and locations, and we expect authors to adhere to the NeurlPS Code of Ethics and the guidelines for their institution.
For initial submissions, do not include any information that would break anonymity (if applicable), such the institution conducting the review.
16. Declaration of LLM usage
Question: Does the paper describe the usage of LLMs if it is an important, original, or non-standard component of the core methods in this research? Note that if the LLM is used only for writing, editing, or formatting purposes and does not impact the core methodology, scientific rigorousness, originality the research, declaration is not required.
Answer: [NA]
Justification: The core method development in this research does not involve LLMs any important, original, or non-standard components.
Guidelines:
The answer NA means that the method development in this research does not involve LLMs as any important, original, or non-standard components.
Please refer to our LLM policy (https: //tppr/////uriss.... cwaeeiippsccc/ooote for what should or should be described.
31


Reviews:
[{"review_id": "dqO3qbX4w9", "reviewer": "Reviewer_pLGj", "summary": "The authors propose a novel approach that they term “Conditionally Linear Dynamical system” (CLDS) to model the neural population dynamics underlying activity of large groups of simultaneously recorded neurons. Their method complements a growing set of related methods to fit and interpret neural population dynamics. With CLDS, the authors aim to strike a balance between explanatory power and interpretability. Neural activity is modeled with a linear dynamical system (good interpretability) whose parameters are conditioned on task and/or behavioral variables that are defined by the experimenter, which makes the resulting dynamics overall non-linear (good explanatory power). The authors show analytically how CLDS relates to and differs from several previously proposed methods to fit neural population dynamics. They apply their method to three example datasets. One example consists of simulated responses based on a model consistent with CLDS (a test of model/parameter recovery), the remaining two examples are fits to open experimental datasets from mice and monkeys.", "questions": "(1) For a given choice of conditional variables u_k, the definition of the models includes several meta parameters like the dimensionality of the latent space and the parameters of the GP kernel. The authors mention that some of these parameters are optimized with cross-validation, while others appear to be fixed ad-hoc (e.g. the parameter L). Fig. 6 only provides limited insights into this cross-validation procedure for a single dataset. The latent dimensionality D employed in the examples if quite low, 2D (?) in the mice dataset and 5D in the primate dataset. These values seem low compared to past estimates of the dimensionality of neural activity, which in similar primate datasets was often found to be in the range of 10-15 dimensions. One worry is that such a discrepancy may imply that CLDS is more data-hungry than the authors suggest, at least for the employed parameterization of the Gaussian Processes in those examples. More generally, the demonstration of a robust, systematic approach for determining all meta-parameters would seem to be required to make this method broadly usable in practice. \n\n(2) The linear models that are a component of the proposed methods are known to be ill-posed, in that “dynamics” and “input” components (A and b in Eq 1a) can be non-identifiable, i.e. different combinations of inputs and dynamics can explain the observed activity equally well. This issue of non-identifiability is separate to the one between the latent and observation stages discussed by the authors. What is the implication of this non-identifiability for the application and interpretation of CLDS fits? I would guess that the addition of the conditional dependency of the linear dynamics on u_k could create an additional “level” of non-identifiability, in that it is not clear what variables u_k should be included in a model, and how exactly they should be parameterized (see also my comments below). Does this second type of non-identifiability occur in the model? If so, do the authors have a strategy on how to deal with it?\n\n(3) Fig. 2 show the main validation of the method on simulated data. These results are encouraging, but the scope of this validation is quite limited. The authors consider one specific type of 2D latent dynamics that was produced with the CLDS model and ask if the model can retrieve the underlying parameters. However, it is not clear how well this example matches the experimental realities. For example, what happens to the fits when the observation noise is increased? It is currently not clear how if the simulated observation noise matches those required to explain neural data. Also, what happens if there is a mismatch between the structure of the ground truth model (inputs, u_k) and the structure of the fitted model? Can the model reliably retrieve other types of dynamics, including higher dimensional dynamics, rotational dynamics, or non-normal dynamics? \n\n(4) The authors say that they expect CLDS to provide a good estimate of the “composite” non-linear dynamical system when u_k and x_t “tightly” co-determine each other. But how would one decide if the later condition applies? In a setting that is more akin to a switching LDS, with switches determined by some latent variable that is not obviously linked to a specific location in state-space, CLDS would presumably fail to account for the resulting dynamics. Is there some diagnostic feature of the data or the model that could be used to identify such failures? \n\n(5) I was confused by the non-linear dynamics in Fig. 2b. The ground-truth model appears to be a CLDS model, so the dynamics in Fig 2b does not seem to “appear” anywhere in the definition of the ground truth model. In fact, it would be interesting to see how well CLSD can fit data that was generated from a ground-truth model that exactly matches that in Fig. 2b. \n\n(6) The fits of CLDS to the example experimental datasets are also promising, but again the model validation is quite limited. Fig. 3b shows that a few neurons’ dependency on theta is captured, and that for many neurons the mean activity is captured. But what else can be gained from that plot? Likewise, Fig. 3c shows that for a few neurons the tuning curve for theta is well captured. But what about all the other neurons? Dynamics as in Fig 4 are typically shown as condition-averaged trajectories (see for example the original LFADS paper). How well does the model capture those, e.g. compared to LFADS? Also, how exactly was the 3d-space in Fig. 4b defined?\n\n(7) One peculiar choice in the fits of Fig. 4 is the definition of the condition variable z_t, which switches from 0 to 1 100ms after reach onset. Notably, this variable has no immediate parallel in the task, but rather seems to have been chosen by the authors to match the time at which firing rates across the population consistently increase. I find this choice problematic, because it involves a degree of “hand-engineering” the solution found by CLDS, in that it assumes that if a change in dynamics is happening during a trial, it must happen at that (human-chosen) time. I would have found this example more convincing if the authors had chosen time within trial as one of the conditional variables and had then shown that CLDS recovers a change in dynamics (and potentially inputs) sometime after the go cue.", "limitations": "See questions above.", "rating": 4, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 4, "significance": 3, "originality": 3, "strengths_and_weaknesses": "The presented methods (including analytical derivations of expectation-maximization steps to fit the model in several settings) are described in detail and could be an interesting, original addition to current methods. But while the proposed method seems promising, the results appear quite preliminary. The three example applications are limited in scope and are somewhat insufficient to fully assess the value of the method. As I explain in more detail below, several important controls are missing and in general the limits of the models and fitting procedures are not sufficiently explored and discussed.", "was_revised": true, "final_justification": "The authors' answers have clarified several of the analyses and findings presented on the simulated and experimental datasets included in the original submission. The authors also seem to have made some effort to more rigorously test the limits of their method and how it may fail in practice, although I find it somewhat difficult to fully assess these novel contributions without seeing the revised manuscript (which is not possible based on current Neurips rules). I believe such tests will be critical to the ultimate success of this method. Nonetheless, based on the clarifications and improvements shown in the responses I have increased my rating from 3 to 4.", "rebuttal": "We are very grateful to the reviewer for their detailed and insightful comments. We address your questions point by point below, and point to our reply to Reviewer hqD4 (due to space constraints) which details a \"Mechanistic model of ring-attractor dynamics\" inspired by Question #5 in your review.\n\n## Responses to questions\n\n**(1) Hyperparameters:** Tuning hyperparameters is a common challenge across many models, with cross-validation being the *de facto* approach. Although CLDS models are a relatively simple class of models (e.g., in comparison to any deep network based approach), a comprehensive grid search of all potential hyperparameter settings is still computationally infeasible. Here, we stick to a relatively simple choice of an approximate GP kernel with a single hyperparameter governing smoothness (i.e. the lengthscale). We tune this by cross-validated co-smoothing. \n\nWhile it is possible to treat $L$ as a hyperparameter tuned to improve performance, we propose to instead set $L$ to a large enough value that results in negligible error in the GP kernel approximation. As $L \\to \\infty$, we pay a larger computational price in fitting the model (because we have more basis functions), but the statistical properties of the model are essentially unchanged because high-frequency basis functions have nearly zero amplitude. This idea is well established in GP regression literature; see for example (Greengard et al. 2025; \"Equispaced Fourier representations for efficient Gaussian process regression from a billion data points.\"). Hopefully, these clarifications make it clear that $L$ was not selected in an “ad hoc” manner, but set to a “large enough” value. We will explain this in more detail in the revision, as well as numerical error metrics. \n\nRegarding dimensionality, we are unaware of conclusive evidence that 10-15 dimensions is optimal to model nonhuman primate motor cortex. Could the reviewer give us a more concrete citation? In our hands, we find that we get substantially diminishing returns after fitting a ~5D model whether the model is a classic LDS or a CLDS. We report in **Table 1** below co-smoothing results as we take $D \\in \\\\{3, 5, 10, 15\\\\}$.  In our revision, we will include a more comprehensive sweep over this hyperparameter as a supplementary figure.\n\n| **D** | 3 | 5 | 10 | 15 |\n| ---| ---| --- | --- | --- |\n|**CLDS** | 0.229 | 0.232 | 0.207 | 0.168|\n|**LDS** | 0.211 | 0.193 | 0.196 | 0.17 |\n\n**Table 1**: Co-smoothing reconstruction of single held-out neurons from the test-set, over varying latent dimensionality $D$, averaged over two random seeds. \n\n**(2) Identifiability:** Our revision will clarify the set of equivalent CLDS model solutions, known from the literature on LDS models (Glover & Willems, 1974, “Parametrizations of linear dynamical systems: Canonical forms and identifiability”). We note that this form of non-identifiability is relatively mild compared with many alternative approaches. Anything nonlinear dynamical systems model or deep neural network approach (e.g. LFADS) will have much more extreme forms of non-identifiability.\n\nAlthough the latent space of CLDS is only identifiable up to these transforms, other quantities are uniquely identified. First, eigenvalues of the dynamics matrices are also identifiable (and thus the dynamical regimes), which we make sure to convey and asses recovery on the two head-direction tasks. Second, for a fixed set of inputs $u_1, \\dots , u_T$ and a fixed estimate of CLDS model parameters the conditional distribution of observations $y_1, \\dots, y_T$ is an $NT$-dimensional Gaussian with identifiable mean and covariance. The parameters of this Gaussian are efficiently computable via Kalman smoothing. Furthermore, these parameters have scientific value. For example, the covariance characterizes across-neuron and across-time noise correlations. See, for example, Panzeri et al. (2022, \"The structures and functions of correlations in neural population codes.\").\n\n**(3) Nonlinear experiment:** Thank you for this insightful perspective, it highlights an opportunity for improving the soundness of our approach. Our motivation behind Figure 2 was to show how a canonical nonlinear neural circuit (a ring attractor) could be instantiated by a conditionally linear model. We show (not surprisingly) that we can recover the parameters of this model when fit to noisy simulated observations. Now, as prompted by you and Reviewer hqD4, we sought to challenge the inference and **devised a new task**, specifically chosen to reflect well-established high-dimensional and non-normal (chain-like) dynamics of interest in the neuroscience community. The ground-truth model is not a CLDS, providing a model-mismatch experiment. We found that the CLDS was able to reliably uncover the ring-attractor structure of the true model, further attesting to its ability to recover the nonlinear structure we've seen in the mouse head-direction data.\n\nFinally, like any statistical model, CLDS models will fail to provide good fits in regimes where the data is too noisy or scarce. Data requirements will grow with the complexity and dimensionality of the underlying system. Unlike nonlinear/deep network alternatives, linear systems are amenable to a theoretical analysis (although this remains to be an open area of research). See, for example, Hardt et al. (2018) \"Gradient descent learns linear dynamical systems\", which shows that no local minima exist under certain conditions. Some sample complexity results are also available. See, Tsiamis et al. (2023) “Statistical learning theory for control: A finite-sample perspective”.\n\n**(4) u-x diagnostic metric:** Thank you for this important question; it is challenging but essential to get a good grasp of this relationship, given the importance of inputs in our modeling approach. Already present in the manuscript, we’ve derived in Appendix A.2. bounds to our approximation error for both the CLDS approximation to a “true” nonlinear system (eq. 23) and the composite dynamics to a “true” autonomous system (eq. 42). Both of these bounds revolve around the second moments of the conditions distributions $p(x | u)$ and $p(u | x)$, namely $\\xi :=\\mathbb{E}_{u}[\\mathrm{Tr}[\\mathrm{Cov}[x | u_t]]$. We are pleased to **provide numerical estimates of this quantity** in **Table 2** below, computed empirically from posterior estimates in a manner akin to the composite dynamics detailed in the text. A control is obtained by performing the same procedure on conditions shuffled over trials/batches. \n\n| | Fit | Control |\n| ---| ---| --- |\n|Synthetic HD | 5.7 | 13.3|\n| Mouse HD | 8.4 | 30.7|\n| Macaque reaching | 6.4 | 7.6 |\n\n**Table 2**: empirical estimates of $\\xi$. The more co-dependent these variables are, the smaller $\\xi$ will be, and the better our CLDS approximation is to a ground-truth nonlinear system. \n\nTheir relative value and difference to the control help provide this diagnostic feature. These quantities will be included in the updated manuscript, and the methods along with the open-source code.  Finally, regarding the SLDS example, we refer to points (1) and (7) on general model selection.\n\n**(5) CLDS decompose nonlinear dynamics** It is worth noting that Fig. 2b represents the ground truth composite dynamics, not the ground-truth model itself. Nonetheless, and directly to your point, we have developed above an additional example where we fit a CLDS to a mechanistic ring-attractor model (see response to hqD4). We hope that our description and the pending figure will satisfy your questions.\n\nPart of our motivation here was to show that a CLDS model can, in fact, be hand-constructed to produce simulated ring attractor dynamics. We believe this approach of constructing nonlinear dynamics from linear dynamics per condition is a conceptual point worth making to the neural modeling community.\n\n**(6)** Regarding the head direction system, we will include a supplemental figure showing all neural tuning curves (all of which are well-fit by the model). Regarding condition-averaged trajectories, thank you for pointing this out, our revision will include them—they look similar for the LFADS and CLDS models. \n\nImportantly, both the tuning curves and condition-averaged trajectories are relatively simple to fit. The hard problem is predicting the statistics of single-trial dynamics and an underlying flow field that is capable of forecasting or predicting heldout activity. This is quantified by the co-smoothing analysis.\n\n**(7) Choice of inputs:** Firing rates of motor cortical neurons rapidly rise at movement onset, see Kaufman et al. (2016, “The Largest Response Component in the Motor Cortex Reflects Movement Timing but Not Movement Type”). In principle, we could have set the step function to coincide with the movement onset of the hand, a standard way to align spike times in reaching tasks, but unfortunately, the data we’re working with didn’t have this information.\n\nThe CLDS model does involve a degree of “hand-engineering.” As we explain in the paper, we view this as both a strength and a limitation of the model. It is a limitation because we might miss out on important latent variables by taking a more supervised approach. But, on the other hand, purely unsupervised methods are difficult to fit. Our view is that there is a lack of tools that are easy-to-use and allow users to incorporate domain knowledge into their model. CLDS models fill this gap in the literature. \n\nIt is true that the inputs to the model can be engineered in a variety of ways, and therefore fair to ask, “How should one choose $u$”? We can’t give a general answer to this question— the answer will depend on the specifics of each experiment. If a practitioner identifies several possible choices for a particular experiment, we suggest that cross-validation could be used to select between these alternatives.\n\nThank you again for your time; we look forward to your final assessment."}, {"review_id": "jH320I00RZ", "reviewer": "Reviewer_JFre", "summary": "This paper introduces a framework named the CLDS: Conditionally Linear Dynamical Systems, a novel probabilistic modeling framework for the neural population activity along with the task-related or the behavioral-related covariates. CLDS first have the interpretability and tractability of classical linear dynamical systems. And also, it has the flexibility of modern non-linear powerful models by allowing system parameters (e.g., transition and emission matrices) to vary as a function of observed covariates via the pre-defined Gaussian Process (GP) priors. Like most works in this probabilistic modeling genre, the CLDS derives an efficient MAP inference method using approximate GP basis expansions and EM optimization. As for the experiment part, this work has done extensive results on both synthetic and real neural datasets of mouse and macaque. The results manifest that CLDS achieves superior performance aganist the baseline methods, especially in data-limited situations.", "questions": "1. What's the most difference and novelty of your probabilistic modeling framework compared to existing works? If there is a table that compare between these methods, that would be apparent to the audience. Also, there is a recent paper [1] whose method that can be compared with the proposed CLDS.\n2. How does u_t enables non-linearity and how much could it introduce?\n3. Could this framework extend to non-Gaussian likelihoods?\n\n[1] Multi-Region Markovian Gaussian Process: An Efficient Method to Discover Directional Communications Across Multiple Brain Regions. Li, et, al. ICML 2024.", "limitations": "Please refer to my 'Questions' section.", "rating": 5, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 3, "significance": 2, "originality": 3, "strengths_and_weaknesses": "I think that the strengths of this work lie in its new probabilistic modeling framework with certain extent of locally interpretation. The idea of adding some conditional covairates to the traditional LDS framework is interesting, which lacks some expressive power. With some new derivations, the work also presents the framework with closed-form EM updates. As for the empirical experimental results, the model demonstrates good performance on both synthetic and multiple real-world neural datasets. On the other hand, I think the concerns of this method lies in what's the background neuroscience intuition behind this proposed elegant framework, what part of the method is brain-inspired. How could you prove the soundness of CLDS besides the experimental results.", "was_revised": true, "final_justification": "See my full comments below. I maintain my opinion that the manuscript is suitable for publication and keep my score at 5.", "rebuttal": "We are grateful to the reviewer for their time and constructive comments. We are pleased to see your overall positive assessment of our work! Below are the answers to your questions.\n\n### Weaknesses\n\nThe CLDS model is a statistical description of condition-dependent high-dimensional time series with latent structure, primarily motivated by and dedicated to the analysis of neural dynamics. We would not consider the CLDS to be “brain-inspired” in the sense of a biophysically realistic circuit model. \n\nRegarding the soundness of the approach besides experimental results, we would point to the fact that the CLDS model is a natural extension of a lot of existing, influential models. In particular, we would point to past work using linear network models of neural circuits (see citations [4-8] in the paper) and more recent work that develops nonlinear dynamical systems models (e.g. citations [9-14] in the paper). Our central argument is that CLDS models provide some of the nice benefits of both.\n\nIn practice, neuroscientists are often interested in finding the relationship between experimental task variables or observed behavioral covariates with neural data. Common models (e.g. LFADS) are fit either in an entirely unsupervised manner or incorporate this task and behavioral covariates in a highly complex, nonlinear manner. CLDS models strike a nice balance between interpretability and expressivity.\n\n### Questions\n\nAnswering your questions in order\n1. There is a very broad variety of statistical models for neural data analysis that are targeted at different specific circumstances. While it is not easy to fit all of these methods into a table, we think the most important advance of our paper is to show that simple and interpretable models (i.e. linear) can be tweaked to model complex (i.e. nonlinear) neural circuit dynamics.\n2. $\\mathbf{A}(\\boldsymbol{u}_t)$, $\\mathbf{b}(\\boldsymbol{u}_t)$, $\\mathbf{C}(\\boldsymbol{u}_t)$, and $\\mathbf{d}(\\boldsymbol{u}_t)$ are all nonlinear functions of $\\boldsymbol{u}_t$ and belong to the nonparametric family of functions determined by a Gaussian Process (GP). To your question, with the coefficients framed as GPs, they allow for a nonlinear dependency on the conditions, making the system linear in dynamics and *nonlinear* in conditions. The amount of nonlinearity introduced depends on the kernel and its parameters. For example, if we use RBF kernels, the length scale determines the probability of how much high-frequency content (or wiggliness) is in the function. This provides an appropriate mechanism for controlling the level of smoothness either using some prior knowledge or in a data-driven manner.\n3. Certainly! We comment on how this can be done on lines 163-170 in the submitted manuscript.\n\nFinally, thank you for providing a pointer to this missed reference; we will make sure to incorporate it. \n\nThank you again for your time; we look forward to your final assessment."}, {"review_id": "CRzN973rIR", "reviewer": "Reviewer_PoZQ", "summary": "In this paper, the authors introduce a class of latent dynamics models for modeling neural activity. In particular, they propose time-varying linear dynamical systems, where the dynamics and readout parameters vary smoothly as a function of external inputs at each point in time. They demonstrate how these models can be efficiently trained, and show that they are performant on both synthetic and real data tasks.", "questions": "- Can the authors comment on the tradeoff between latent dimension $D$  and approximate GP prior truncation length $L$ in terms of expressivity and risk of overfitting?\n- One interpretation of the parameterization $\\textbf{b}(\\textbf{u}_t)$ is that the external inputs undergo a nonlinear featurization before being linearly read into the state updates (setting aside the dynamics matrix also depending on $\\textbf{u}_t$ for now). How well would other models do with this nonlinear input featurization step? I would imagine that even a standard LDS model could improve greatly on nonlinear tasks with this modification.", "limitations": "yes", "rating": 5, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 4, "clarity": 4, "significance": 3, "originality": 3, "strengths_and_weaknesses": "Strengths:\n- The authors convincingly show that both training and inference of CLDS models can be efficiently done via classical techniques such as Kalman smoothing and EM. \n- There are many avenues towards interpretibility of CLDS models (e.g. linearizing around a particular input, or marginalizing out the inputs to directly yield possibly nonlinear composite dynamics)\n\nWeaknesses:\n- While the authors do a great job comparing CLDS to models traditionally used in neuroscience (such as (r)SLDS), the class of models introduced here---linear RNNs with parameters being time-varying solely through input dependence---is not at all new in the ML literature, and is a popular formulation of modern selective state space models (e.g. Mamba). The authors should comment on this, and perhaps discuss how the Bayesian perspective (which is indeed distinct from that setting), or the distinct parameterization strategy for achieving input dependence (truncated approximate GP prior as opposed to an MLP), yield unique advantages.", "was_revised": true, "final_justification": "I am keeping my original positive score. There were no major issues raised.", "rebuttal": "We are very grateful to the reviewer for their time and insightful questions. We are pleased to see your overall positive assessment of our work! \n\n### Weaknesses\n\nThe point on modern state-space models is a great one. We had mentioned input-switched affine networks in the first submission, but missed Mamba and its related literature. Our revision will include a longer discussion emphasizing the following:\n\n- **Presence of noise**: SSM models developed in the ML community do not model noise in the latent space. This is thought to be important when analyzing neural data due to limited dataset sizes and noisy observations. Latent noise complicates the inference process and common strategies for inference include EM and marginal optimization for tractable models and variational inference or MCMC for intractable models.\n- **Structure of the state-space**: SSM models developed in the ML community are often applied to natural language processing tasks and therefore take in discrete tokenized inputs. The CLDS model we propose takes inputs over a *continuous* space and furthermore incorporates a smoothness prior over this space, which is an important inductive bias. Thus, in the head direction example, our model supports an exact ring attractor. Furthermore, our work shows how input-selective models can be used to model continuous attractors (e.g. ring attractors for heading direction). Given the fundamental importance of attractor dynamics to modern theories of neural circuit dynamics, we think this is an important conceptual advance which cannot be found in existing ML literature.\n- **Models of nonlinearity**: Finally, we can comment on the GP v.s. MLP modeling. GPs importantly carry a prior over functions, which is of great help in inferring the non-parametric forms of the coefficients in the low-data regimes often seen in neuroscience (as mentioned above). Furthermore, this same prior makes for the potential to obtain full posterior estimates over these coefficients modeled with GPs—we discussed in the paper how this could be done in Lines 127-128. \n\n### Questions\n\nAddressing your questions one-by-one:\n- **Dimensionality and GP trade-off**: An important aspect of this tradeoff is that the CLDS is *linear* in dynamics ($D$-dimensional) and is nonlinear over conditions (with degree of nonlinearity enforced by $L$). Hence, while both trade off in final expressivity, they do so in complementary ways that model predictions (in time v.s. over different conditions) can help further parse out. Second, on the topic of overfitting, we typically treat $L$ as fixed to a high-enough value to allow sufficient expressivity by the GP prior, as opposed to an optimized hyperparameter—see also our related comments to Reviewer pLGj on this topic. The dimension $D$ is typically variable and searched over, but selected through cross-validation and restricted to be low-dimensional, both avoiding overfitting. \n\n- **Nonlinear biases**: We agree with your intuition here. It is worth noting that the specific model mentioned (standard LDS + nonlinear $\\mathbf{b}(\\boldsymbol{u}_t)$) is a subcase of the more general CLDS with both $\\mathbf{A}(\\cdot)$ and $\\mathbf{b}(\\cdot)$ nonlinear: while the GP prior implicity favors against that, one *could* learn a degenerate set of weights emphasizing the constant function for $\\mathbf{A}(\\cdot)$, making it constant over $\\boldsymbol{u}$. Still on the topic of biases, one can bypass the dynamics entirely by using $\\mathbf{d}(\\boldsymbol{u}_t)$ to essentially fit a nonlinear tuning curve or PSTH to each neuron (disregarding how single-trial dynamics vary around these trial-average firing rates). We think it is feasible and interesting to see how all of these “ablations” of the CLDS model perform on various experimental datasets, in particular the nonhuman primate reaching dataset.\n\nThank you again for your time; we look forward to your final assessment."}, {"review_id": "2LeqW3dEyy", "reviewer": "Reviewer_hqD4", "summary": "The authors introduce conditionally linear dynamical systems (CLDS), where system parameters are modulated nonlinearly as a function of observed covariates $u_t$ (e.g., via a kernel mapping). This allows the latent dynamics to flexibly vary across contexts or conditions, while retaining the interpretability and tractability of linear state-space models. The model is fit using an approximate GP prior over parameters, enabling closed-form EM updates under Gaussian noise assumptions.", "questions": "*Comparisons/Related Work*:\nThe proposed approach shares a lot of conceptual similarities linear/structured SSMs, particularly to Mamba, where the recurrence kernel is conditioned on the input sequence, yielding a similarly dynamically reparameterized linear system. While the technical motivations differ (e.g., sequence modeling vs. spike dynamics), it could be beneficial to include this in the related work discussion, also with respect to training algorithms, benefits of linearity etc. Another recent line of work (Almost-Linear RNNs, Brenner et al., NeurIPS 2024) bears close conceptual similarity to recurrent switching linear DS, and could be mentioned.", "limitations": "Clear limitations section.", "rating": 5, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 3, "significance": 3, "originality": 3, "strengths_and_weaknesses": "*Strengths*:\nOverall, I like this paper. It’s clearly written, easy to follow, and communicates its ideas well. The motivation for combining input-dependent parameter modulation with linear dynamics is compelling, especially in neuroscience where interpretability matters. The figures are informative, and the structure of the paper makes it straightforward to trace the reasoning from model definition to inference to experiments. The theoretical results in Appendix A.2, especially the bounds on the approximation error for composite dynamics (A.2.2), are a valuable addition. Overall, it’s a useful contribution that sits in an important design space between black-box flexibility and mechanistic transparency.\n\n*Weaknesses*:\n\n*Overall Weakness*:\n While the paper presents a clean and interpretable framework, its empirical scope feels slightly narrow, and the evaluation does not seem to probe the limits or broader applicability of the method. All experiments are tailored to idealized conditions that align well with the model’s assumptions (Gaussian observations, well-structured modulators, modest data scale, and low latent dimensionality). Many of the paper’s selling points (e.g., “easy inference,” “simple extensions,” interpretability) are stated as general advantages, but it feels like they are not properly stress-tested or fully explored in more challenging settings. This leaves open the question of whether the performance gains and interpretability benefits are worth the constraints introduced by the modeling assumptions.\n\nThat said, if the authors are able to convincingly address these concerns, e.g. through additional ablations, clarifications of model limitations, and a deeper discussion of scalability and generalizability, I would be happy to support acceptance.\n\n*Spike Data*:\nThe authors emphasize that one of the key advantages of CLDS models is ease of inference, particularly under a Gaussian noise model where Kalman smoothing and closed-form EM enable efficient learning (§2.3). They also note that the extension to non-Gaussian observations (e.g., Poisson spike counts) is straightforward in principle, since the posterior remains log-concave and tractable. However, all real data applications, including the macaque spike recordings, still use the Gaussian observation model. While this is mentioned in the limitations, this feels like a missed opportunity, particularly since the whole paper is tailored to \"modeling neural activity\": if extending to count-based models is practical and aligns better with the data-generating process, it would have been useful to demonstrate this empirically. At minimum, a concrete discussion of what would be required to implement this extension (and what challenges, if any, prevented it here) would help clarify how straightforward this path really is. As it stands, the claim of “easy extensibility” feels somewhat contradicted by its absence in the main experiments, where it would fit perfectly.\n\n*EM Scaling*:\nAs a related point (given EM is highlighted as a major strength of the CLDS framework) while EM is indeed attractive for small- to medium-scale models, it is well known to scale badly with increasing latent dimensionality or large datasets (which would become relevant in many modern neural recording settings). The authors mention that “in principle” the model could be extended to use variational inference or MCMC methods, but do not test these alternatives or provide a discussion of where EM might become a bottleneck. An analysis of scalability (e.g., in terms of latent dimensionality, trial count, or observation length) would strengthen the paper’s claims about practical applicability and robustness. Even a partial exploration of variational alternatives would also help contextualize how “easy to fit” this model remains.\n\n*Kernel function*:\nThe use of an approximate GP prior over system parameters is presented as the paper’s main technical contribution. However, the empirical evaluation does not really explore how the structure or complexity of the kernel, or the form of the input u, affects performance. In the experiments, the modulators are closely linked to the observed neural dynamics (e.g., categorical or smooth latent variables), making it difficult to assess the generality or necessity of the kernel. I could image that there are forms of inputs that would benefit from different kernels. Since the GP prior induces a specific form of smoothness in the system parameters, it could strengthen the paper to include ablations comparing different modulator structures and kernel variants (e.g., alternative basis functions or GP approximations), especially given the centrality of this mechanism to the paper's novelty. For instance, what prior would you recommend if smoothness assumptions over system parameters are violated (e.g. task switches or abrupt transitions)? How do find an optimal $\\kappa$ depending on smoothness?\n\n*Input-driven DS*:\nThe paper briefly mentions an important conceptual parallel to time-varying or input-driven LDSs, where external inputs u_t​ enter the system via a factor-loading matrix B. This structure is known to induce time-dependent fixed points or more complex attractors (e.g., ring attractors) if the input signal is sufficiently expressive. While not explicitly framed this way, the proposed model can be viewed as implicitly implementing a snapshot attractor. Given the close mathematical similarity to classical input-driven LDSs, it would be helpful to more explicitly articulate what is novel about the kernel-based parameterization (particularly e.g. compared to Input Switched Affine Network that are cited by the authors). At least to my understanding, similar expressive dynamics could potentially be captured by using known u_t​ in combination with a flexible encoder. A systematic empirical comparison to such baselines could help clarify this distinction and further strengthen the case for the proposed approach.", "was_revised": true, "final_justification": "An overall interesting and meaningful contribution. I already had a positive impression of the paper. The authors mostly adressed my concerns my concerns with novel experiments and clarifications. One small remaining open point I raised (a different noise model) is not an integral part of this paper's contribution, and does not preclude me from voting for acceptance.", "rebuttal": "We are very grateful to the reviewer for their detailed and insightful comments, and overall positive review. You and other reviewers found merit in the conceptual advance, found the paper well written, and suggested better outlining and testing the scopes of applicability. This provides valuable feedback towards improving our contribution. \n\nThank you for pointing out the scope of the experiments conducted so far in the paper, which we hope to highlight model recovery and inference on real-world experiments. Your feedback on more challenging experiments targeting the model assumptions is echoed by reviewers and is a valuable concern. To address this overall weakness, we introduce below a new experiment that readily fits within the experiments considered so far while targeting the idealized conditions of low dimensionality and model alignment. Furthermore, we point you to our answer to Question 4 of Reviewer pLGj, which details **new numerical calculations** of the empirical covariance in $x$ given $u$. As per Appendix A.2., this provides a diagnostic feature and proxy for our approximation error to any true nonlinear dynamics in the system. \n\n## Mechanistic model of ring-attractor dynamics\n\nInspired by reviewer feedback, we provide an alternate synthetic ring attractor experiment based on the neuroscience literature—see review by Hulse and Jayaraman (2020, “Mechanisms Underlying the Neural Computation of Head Direction”). The intention with this additional task is to test the CLDS and its inference in a synthetic data modality, where we know the underlying computation, but which, importantly, is not generated from a CLDS model, to test the model in mispecified settings. \n\n**Model description**: We write a model of continuous ring attractor dynamics with bumps of activity integrating angular velocity (Zhang, 1996, “Representation of spatial orientation by the intrinsic dynamics of the head-direction cell ensemble: a theory”); a model and computational account of spatial orientation shared across many species. In this model, attractor dynamics are implemented through head-direction preferential units arranged in a “ring”, with short-range excitation and long-range inhibition forming a bump of activity for a specific unit in a specific head-direction, with a velocity integrator causing a shift in the location of the bump as the head-direction changes. Concretely, the dynamics of the new generative model follow\n\n$$ \\frac{1}{\\tau} \\dot u(\\phi, t) = - u + f(w * u) + v \\frac{\\partial}{\\partial \\phi} u(\\phi, t).$$\n\nHere, $u(\\phi, t)$ represents the network activity of cells with preferred head-direction $\\phi$ at time $t$, $f$ is a (ReLU) nonlinearity, $w$ is a “Mexican-hat” convolution filter, and $v_t = \\dot \\theta_t$ is the instantaneous difference in head-direction $\\theta_t$. For implementation purposes, $u$ is discretized to $u_t \\in \\mathbb{R}^N$, $t \\in \\{1, \\dots, T\\}$, with each unit  $u_i$ having preferred directions $\\phi_i$ regularly spanning the interval $[0, 2\\pi)$, and $w * u$ is implemented with a circulant matrix $Wu_t$. We further include additive Gaussian random noise to make the process stochastic. We use this activity $u$ as our observations, of dimension $N=32$. The head direction $\\theta_t$ still acts as our conditions, sampled the same way as in our previous synthetic HD model.\n\n**Results**: We know the underlying model captures high-dimensional (N=32 here) ring attractor dynamics. When fitting a low-dimensional CLDS model, we expect the latent dynamics to capture the ring-attractor structure and the nonlinearity of the dynamics, with a linear emission model bringing it back to the high-dimensional activity. This is precisely what we found when fitting the CLDS—all necessary figures to relay the results will be included in the final revision. After a search over kernel hyperparameters, we found the best-performing model to encode a ring of fixed points that closely resembles the synthetic HD fixed points, and with eigenvalues closer to the ones observed from the CLDS fits on the mouse ADn recordings. The results make us confident that the CLDS can capture this nonlinear structure even under model mismatch.  \n\n## Response to individual weaknesses\n\n**(1) Spike data**: We agree that the extension to additional noise models (e.g. Poisson) is of interest, and we intend to follow up soon with a full implementation of this. There are no technical hurdles in principle—we just have not prioritized this particular direction. In our experience, simple linear Gaussian methods perform well and are popular among practitioners in neuroscience. For example, Gaussian Process Factor Analysis (GPFA) uses a Gaussian noise model.\n\nWe agree with the reviewer that this extension would enhance the paper. However, we think the core conceptual advances of the paper do not rely on the noise model. For example, we think most neuroscientists would regard the head direction system as an inherently nonlinear system. Here, we show that one can model it using linear methods (assuming that head direction is simultaneously measured with neural activity).\n\n**(2) EM Scaling**: For now, we are focused on obtaining point estimates of model parameters &mdash; i.e. the functions $A(u), b(u), C(u), d(u)$. We mentioned MCMC and variational inference as possible extensions for obtaining an approximate posteriors over these parameters/functions. In contrast, EM only provides point estimates (the maximum of the posterior). In other words, this is an apples-to-oranges comparison. Relative to EM, we would expect that MCMC to be much more expensive (since it involves sampling) and variational inference to be moderately more expensive (since it involves an optimization problem with additional parameters).\n\nWe think that EM is scalable for problems of interest in neuroscience. Extremely high-dimensional latent spaces would be susceptible to overfitting anyways. Furthermore, it is not clear that alternative methods are more scalable. In particular, any method that computes the marginal log likelihood of $y_1, \\dots, y_T$ integrating over $x_1, \\dots, x_T$ will perform the Kalman filtering step anyways, which accounts for a large chunk of the computation time in EM. In particular, the M-step is very cheap. In principle, it is possible to do gradient ascent on the marginal log likelihood, but then one must tune learning rates and in general take much smaller steps in parameter space.\n\nIn practice, we believe the biggest challenge for scaling will be datasets of very long duration (large $T$). Here, it is possible to leverage parallel computation on a GPU to perform the Kalman filtering, reducing the scaling behavior from $T$ to $\\log T$. See Sarkka & Garcia-Fernandez (2020, “Temporal parallelization of Bayesian smoothers”)\n\n**(3) Kernel function**: Draws from a GP with an RBF kernel are infinitely differentiable. A Matern kernel would be a weaker prior (leading to less smooth functions). This would introduce additional hyperparameters and lead to a complex model selection process. For example, the kernel hyperparameters are known to be unidentifiable in a simple regression setting. See Zhang (2004, “Inconsistent estimation and asymptotically equal interpolations in model-based geostatistics”).\n\nA CLDS model can account for abrupt task switches or transitions if they are associated with a measurable behavioral event. For example, we show this in our analysis of the nonhuman primate reaching task by selecting $u_1, \\dots, u_T$ to be a step function centered at the time of movement onset.\n\nIf abrupt task switches or transitions are unobserved (i.e. latent variables), then one could consider a switching linear dynamical system. We discuss the complementary challenges and advantages of this alternative model in the paper. In particular, inference in this model is less straightforward because Kalman filtering and smoothing is not possible. Speaking from personal experience, we found that open-source implementations of these models were quite difficult to tune in practice. \n\n**(4) Input-driven DS**: The most relevant prior work we are aware of is by Foerster et al. (2017, “Input Switched Affine Networks: An RNN Architecture Designed for Interpretability”) and more recently the SSM models like Mamba (Gu & Dao, 2024, “Mamba: Linear-Time Sequence Modeling with Selective State Spaces”). Neither of these papers explicitly mention attractor dynamics or ring attractors specifically.\n\nThe main difference is that these models are applied to natural language processing tasks and therefore take in discrete tokenized inputs. The CLDS model we propose takes inputs over a *continuous* space and furthermore incorporates a smoothness prior over this space. Thus, in the head direction example, our model supports an exact ring attractor.\nWe think that an input-switched affine network could in theory learn to approximate a ring attractor if the head direction variable were binned onto a discrete number of values. However, (1) it does not incorporate noise in the latent space, a crucial component of SSMs and probabilistic modeling, and (2) the model would have a hard time learning this from a random initialization because it has no smoothness prior and therefore no inductive bias towards learning similar representations for nearby heading directions.\n\nThank you again for your time; we look forward to your final assessment."}]

Author rebuttal:

Meta-review: In this paper, the authors introduce a nice novel idea of recasting nonlinear neuronal dynamics in terms of conditionally linear dynamics, given behavioral and task variables. This yields highly interpretable and tractable models for analyzing dynamics in neural systems, without sacrificing the nonlinear nature of neuronal dynamics. 

With the initial submission there were considerable concerns, brought up by two referees, about the rather limited empirical evaluation, lacking important controls and raising issues about data requirements and model scaling. Other questions concerned the determination of the model’s hyperparameters, model identifiability, and some model diagnostics. These were mostly addressed by the authors during the rebuttal through additional numerical checks, detailed formal arguments, and by adding another evaluation on a ring attractor model. Ref. hqD4 also brought up a very valid point in my mind about the model’s similarity to other input-driven DS/ SSM formulations, which casts a bit of doubt on the novelty of the authors’ approach and should be discussed in the paper’s final version.

Three of the four referees clearly voted for acceptance, while one referee initially remained a bit doubtful, but stated that most of their concerns have been properly addressed and leaned toward acceptance in the end as well. I go with the overall positive assessment and recommend acceptance. I expect the authors to include the clarifications provided in the rebuttal and discussion, and in particular the new evaluation setting, in their final revision!
Decision: Accept (poster)