Inferring stochastic low-rank recurrent neural networks from neural data
Authors: Matthijs Pals
A Erdem Sağtekin
Felix C Pei
Manuel Gloeckler
Jakob H. Macke
Venue: NeurIPS 2024 poster (2024)
Keywords: Low-rank RNNs
dynamical systems
variational inference
sequential monte carlo
neural data
Primary area: neuroscience_and_cognitive_science
TL;DR: We fit low-rank RNNs to neural data using variational SMC, and obtain models that are both generative and have tractable low-dimensional dynamics.

Abstract: A central aim in computational neuroscience is to relate the activity of large populations of neurons to an underlying dynamical system. Models of these neural dynamics should ideally be both interpretable and fit the observed data well. Low-rank recurrent neural networks (RNNs) exhibit such interpretability by having tractable dynamics. However, it is unclear how to best fit low-rank RNNs to data consisting of noisy observations of an underlying stochastic system. Here, we propose to fit stochastic low-rank RNNs with variational sequential Monte Carlo methods. We validate our method on several datasets consisting of both continuous and spiking neural data, where we obtain lower dimensional latent dynamics than current state of the art methods. Additionally, for low-rank models with piecewise linear nonlinearities, we show how to efficiently identify all fixed points in polynomial rather than exponential cost in the number of units, making analysis of the inferred dynamics tractable for large RNNs. Our method both elucidates the dynamical systems underlying experimental recordings and provides a generative model whose trajectories match observed variability.

Inferring stochastic low-rank recurrent neural networks from neural data
Matthijs Pals A Erdem Sagtekin! 122,3 Felix Peil.
Manuel Gloeckler Jakob H Macke
Machine Learning in Science, Excellence Cluster Machine Learning, University of Tübingen, Germany Tübingen Al Center, Tübingen, Germany
GGraduate Training Centre Neuroscience, University of Tübingen, Germany *Department Empirical Inference, Max Planck Institute for Intelligent Systems, Tübingen, Germany
Abstract
A central aim in computational neuroscience is to relate the activity of large populations of neurons to an underlying dynamical system. Models these neural dynamics should ideally be both interpretable and fit the observed data well. Low- rank recurrent neural networks (RNNs) exhibit such interpretability by having tractable dynamics. However, it is unclear how to best fit low-rank RNNs data consisting of noisy observations of an underlying stochastic system. Here, we propose fit stochastic low-rank RNNs with variational sequential Monte Carlo methods. We validate our method on several datasets consisting of both continuous and spiking neural data, where we obtain lower dimensional latent dynamics than current state of the art methods. Additionally, for low-rank models with piecewise-linea. nonlinearities, we show how efficiently identify all fixed points in polynomial rather than exponential cost in the number of units, making analysis of the inferred dynamics tractable for large RNNs. Our method both elucidates the dynamical systems underlying experimental recordings and provides a generative model whose trajectories match observed variability.
Introduction
common goal of many scientific fields is extract the dynamical systems underlying noisy experimental observations. In particular, in neuroscience, much work is devoted to understanding the coordinated firing of neurons as being implemented through underlying dynamical systems (T-5).
Recurrent neural networks (RNNs) constitute common model-class of neural dynamics [6
which can be eegrse-engineered result, several recent research directions have centered on interpretable or analytically tractable RNN architectures. In particular, RNNs with low-rank structure [66 22 admit a direct mapping between high-dimensional population activity and an underlying low-dimensional dynamical system. RNNs
form hypotheses about neural computations [14 [5]. As
with piecewise- activations 9 23 26 are tractable, as they have fixed points and cycles that can be accessed analytically.
To serve as useful models of brain activity, it is important that models also capture the observed brain activity, including ttal-to-trial variability. Many methods that fit RNNs to data are restricted to RNNs with deterministic transitions 6 10 It unlikely that, in general, all variability in the data can be explained by variability in the RNNs initial state. Thus, adopting stochastic transitions is imperative. While probabilistic sequence models are used effectively in neuroscience (77) they have so far largely consisted of state space models without an obvious mechanistic interpretation 28 32
(firstname.cr secoaammmmmmmccmmmmmmmmmeeeee. tammiittubinngen.e
38th Conference on Neural Information Processing Systems (NcurIPS 2024).

---

low-rank RNN
generate
condition
access
(optional)
dynamics
Figure 1: Our goal is to obtain generative models from which we can sample realistic neural data while having tractable underlying dynamical system. We achieve this by fitting stochastic low-rank RNNs with variational sequential Monte Carlo.
Here, we demonstrate that we can fit large stochastic RNNs to noisy high-dimensional data. First, we show that, by combining variational sequential Monte Carlo methods [33 with low-rank RNNs, we can efficiently fit stochastic RNNs with many units by learning the underlying low-dimensional dynamical system. The resulting RNNs are generative models neural data that can be used to sample trajectories of arbitrary length, and also allow for conditional generation with (both time- varying and stationary) inputs. Second, we show that, for low-rank networks with piecewise- -inear activation functions, the resulting dynamics can be efficiently analyzed: In particular, show how all fix points can be found with a polynomial cost in the number of units dramatically more efficient than the exponential cost in the general case.
We first validate our method using several teacher-student setups and show that we recover both the ground truth dynamics and stochasticity. We then fit our model to several real-world datasets, spanning both spiking and continuous data, where we obtain a generative model which needs lower dimensional latent dynamics than current state of the art methods. We also demonstrate how in our low-rank RNNs fixed points can be efficiently - - potentally at lower cost than approximate methods 25 while additionally coming with the guarantee that all fixed points are found.
2 Theory and methods
2.1
Low-rank RNNs
2.1.1 Access to the low-dimensional dynamics underlying large networks
Our goal to infer recurrent neural network models of the form
dx dl
T x (Tx(()
(I)
with neuron activity x(t) E RV time-constant R recurrent weights J e N element-wiss nonlinearity Φ, an R dimensional white noise process (tt) and I'x NxR R In particular, we are interested in the case where the weight matrix has rank R S N, i.e., it can be written with M, N E N (18) Assuming that x(0) lies in the subspace spanned by the columns of M and Fx = MT, with 1, R R we can rewrite Eq. as an equivalent It dimensional system,
M
dz
= z(t) + + ((M((t) + IZE(t))
dt
where we can switch between Eq. and Eq. by means of linear projection, z(t) (MMM)--MX((0 and x(1) = Mz((). Note that we can directly extend these equations to in- clude input, representing, e.g., experimental stimuli context. Even if these stimuli are time-varying, x will be constrained to the span of the input weights and M. By including input to the RNN, we can use the fit models for conditional generation (see Supplement C.2)

---

2.1.2 Low-rank RNNs as state space models
We consider nonlinear latent dynamical systems with observations y::
T
T
p(ZI.T, N1.)) - -p(z)) II
p(ze be-1) TPLr ze);
2
f(z z(-1) = =N((F(z-1)) Mig) p(z1) = Nu 22n)s
ply: z(-1) G(zz),
where the transition distribution is parameterised by discretising low-rank RNN with timestep A., we have mean F(z,) az, + Ñ (MMM), with a = ب and Ñ - N and covariance Σg (see Supplement C.I The specific form of the observation function G. depends on the data- fata-modality, e.g., here we use a Poisson distribution for count observations. This formulation allows one to keep the one-to-one correspondence between RNN units (or a subset of those) and recorded data neurons, (as was desired in previous work, e.g.. 10 12 36). For example, assuming Gaussian observation noise, we simply that X; Mz, and define G MMM,, Eyy)
Once we learn p( VI.T), we can use the obtained RNN as generative model sample trajec- tories, and reverse engineer the underlying dynamics to gain insight in the data generation process. Given the sequential structure of the RNN, we can do model learning by using variational sequential Monte Carlo (also called Particle Filtering) methods 33:35 35
2.2
Model learning with variational sequential Monte Carlo
2.2.1 Sequential Monte Carlo
Sequential Monte Carlo (SMC) can be used to approximate sequences of distributions, such those generated by our RNN, with set of K trajectories of latents Z1:T (commonly called particles) [37]. As we do not have direct access to the posterior p(z; y1:t), we instead sample from a proposal distribution and adjust for the discrepancy between the proposal and target posterior distribution using importance weights. Thus, a crucial choice when doing SMC is picking the right proposal distribution from which we sample latents conditioned on the previous latent Z;-1 and observed data y 1:7, or a subset of those. Given initial samples z7 1:K ~ r and corresponding importance weights k (as defined below) SMC progresses by repeatedly executing the following steps:
resample
Discrete(a.-1 wh-1),
z ~ r(z) Yts =
propose
k ply,, zi (2z-T)) r(zt ye, Ye,ze11
reweight
with w Σ wt Here, the resampling step avoids most of the weights from concentrating on very few partieles. Using SMC, we obtain, at time 1, a filtering approximation to the posterior,
giie (Z Y1:t)
Σ क 8 (z..))
(3)
The unnormalised weights give an unbiased estimate to the marginal likelihood,
T
K
M(y =
(4)
K
1=1
We now detail how we pick the proposal distribution r. For linear Gaussian observations G (Wz,, M,), we set r(z, yt, 2.-1) p(z, yt, 24-1), as this is available in closed form and is optimal (in the sense that minimises the variance of the importance weights
rzz Yr. 1) = A'(I KW) + Ky, Ag)
(5)
3

---

with K the Kalman Gain: K = A2 A and A2 + W W)- (or equivalently K Eg WT (WE, WT Eyy ) and A₂ (I kW) For non-linear observations, we can not invert the observation process in closed form, so we instead jointly optimize a parameterized "encoding" distribution e(z; ye-r':e) (as in a variational autoencoder (38) ). In particular,
assume to be a multivariate normal with diagonal covariance, which we parameterize by causal convolutional neural network, such that each latent is conditioned on the l' latest observations (although sometimes non-causal encoders can be ad advantageous, see Supplement B.5). We then use the following proposal:
r(ze
Zt-1 Vt-era) X e(Zz yr-ra)p(zz (z-1),
(6)
where we now also assume p(z;
Zt-1) has diagonal covariance matrix.
2.2.2 Relationship to Generalised Teacher Forcing
In our approach, the mean of the proposal distribution at time L is a linear combination between the RNN predicted state F(z-1) and data-inferred state ž, A recent study obtained state-of- the art results for reconstructing dynamical systems by fitting deterministic RNNs with a method called Generalised Teacher Forcing (GTF), which similarly linearly interpolates between data- inferred and an RNN predicted state at every time-step [ the model propagates forward in time as Z, = (1 a) F + a²2 Hess et al. B showed that by choosing the appropriate a, one can completely avoid exploding gradients, while still allowing backpropagation through time, and thus obtaining long-term stable solutions [39]. The optimal a can be picked based on the maximum Lyaponuv exponent of the system (a measure of how fast trajectories diverge in a chaotic system).
By including the RNN in the proposal distribution, gation through time through the sampled trajectories. The linear combination is given by a = (E + - where Σz, is the predicted variance of the encoding network. Thus, Iistead of interpolating based on an estimate of how chaotic the system is, our approach combines RNN and data inferred states adaptively (every time step, if Eq. i used) based on how relatively noisy the transition distribution is with respect to the data-inferred states at time L, analogous to, e.g., the gain of Kalman filter. In the formulation of GTF of Hess et al. [ invertable observation model is required. By learning an encoder that predicts a distribution over latents, our method naturally extends to models with non-invertable (e.g., Poisson) observations.
similarly GTF allow backpropa-
in Eq. and similarly in Eq. 6 a a (E₂ -
2.2.3 Variational objective
We can fit our RNNs to data by using SMC to specify a variational objective (33 In variational inference, we specify a family of parameterized distributions Q, and optimize those parameters such that a divergence (usually the KL divergence) between the variational distribution q(21r)) E Q and the true posterior p(Z1.T YI:r) minimized. We do this by maximising a lower bound (ELBO) to the log likelihood p(yr.)).
). In particular, we can use Eq. 4 to specify the ELBO objective B3
L= R Bome (x .atk (yirr)
[lg Mogg(():)))
with Gsme the sampling distribution:
4smc (z a y1:T) = I r(z) y1)
IT=2 r(z) 2 Discreteaat -
During each training iteration, we run SMC, using the closed form optimal proposal (Eq. 5) if
observations are linear Gaussian, otherwise the proposal includes parameterised encoder (Eq. 6)
We can then use the resulting unnormalised importance weights (Eq. 4 to estimate the ELBO, which
we maximise with backpropagation (through time). As suggested in previous studies B3 35 40)
use biased gradients during optimization by dropping high-variance terms arising from the resampling.
2.3 Finding fixed points in piecewise- -ineer low-rank RNNs
After having learned our model, we can gain insight into the mechanisms underlying the data generation process by reverse engineering the learned dynamics [15), e.g., by calculating their fixed points. Here, show that the fixed points can be found analytically and efficiently for low-rank networks with piecewise-linear activation functions. This class of activation functions ((x.) =
(d) Σd bid max(x; h," ,0) includes, e.g., the standard ReLU (d(x.) = max(x, h,, 0))

---

the "clipped" variant (o(x,) = max(x; + h;, 0) - max(x, 0)) with real-world data here.
8 which we used in all experiments
Naively, the cost of finding all fixed points piecewise-linear networks scales exponentially with the number of units in the networks: we would have to solve (1++))* systems of N equations p 24 If
networks are low rank, it is straightforward to show that we can reduce this cost to solving (D + 1 >
systems of R equations (See Supplement [A.I). In addition, however, we show that the computational cost can be greatly reduced further: One can find all fixed points in a cost that is polynomial instead of exponential in the number of units:
Proposition 1. Assume Eq. [ with J of rank R and piecewise- activations c. For fixed rank R and fixed number of basis functions D, we can find all fixed points in the absence of noise, that is all x for which dx de 0, by solving at most (NN)) linear systems of R equations. Proof. See Supplement A
Sketch. Assuming D = 1, activations & max x; -:);; N units will partition
the full phase space into 2 regions in which the dynamics are linear (2 units, 4 regions in Fig. 2. We can thus, in principle, solve for all fixed points by solving all corresponding linear systems of equations
x,=h
( [4 If dynamics are confined the R-dimensional subspace spanned by the columns of M, only subset of the linear regions (3 Fig. 2) can be reached. Each unit partitions the space spanned by the columns of M with a hyperplane (pink points in Fig. (). The amount of linear regions in M, becomes equivalent 'how many regions can we create in R-dimensional space with N hyperplanes?" Using Zaslavsky's 's theorem [41], we can show that this most Σ (+) O(NB) (for fixed Rt.
Figure Proof sketch.
3 Empirical Results
3.1 RNNs recover ground truth dynamics in student-teacher setups
We validated our method using several student-teacher setups (Fig. 3 additional statistics in Fig. $4). We first trained a 'teacher' RNN, with the weight matrix constrained to rank 2, to oscillate. We then simulated multiple trajectories with high level of stochasticity in the latent dynamics (Fig.Ba. top left) and additional additive Gaussian observation noise (Fig. [3.. top right) on the observed neuron activity (y: NNx,,, o.). with Mz). A second 'student" RNN was then fit to the data drawn from the teacher, and both recovered the true latent dynamical system, as well as the right level of
stochasticity (Fig. 3a bottom; Fig. 3d).
We also verified that we can obtain covariance matrices E.g that are numerically close to the ground truth, for teacher networks with various levels of noise (Fig. $5). When using the bootstrap proposal (i.c., sampling from the prior; r = p(z, =1)))), or too few particles, the right level of stochasticity is not obtained, indicating that the use of multiple particles and proposal that conditions observed data is indeed beneficial.
Given that neurons emit action potentials, which are commonly approximated as discrete events, we repeated the initial teacher-student experiment with Poisson observations generated according to y; Pois(softplus( w;x; b;)). The student RNN again recovers the oscillatory latent dynamics. Note that because of the affine transformation in the observation model, the inferred dynamics can be scaled and translated with respect to the teacher model. To verify that samples from our inferred model follow the same distribution as samples from the teacher model, we computed several statistics, which all show close match (Fig. 3e: Fig. $4..
In our final teacher-student setups, we verified the ability to recover dynamics when there are known stimuli or contexts. In particular, we trained rank-2 RNN on task where, at each trial, it receives a transient pulse input corresponding to a particular angle (given as sin(f) cos(f)), and is asked to provide output matching the input after stimulus offset. The teacher RNN learns to perform the task by using an approximate ring attractor which the student RNN accurately infers (Fig. 3D). Here, we inferred all fixed points by making use of Preposition 1. To demonstrate that our method also works when inputs are strongly time-varying, we included an additional setup where the teacher network was asked to report the sign of the mean of noisy stimulus (Fig. S6).
5

---

b
latents
observed rates
latents
observed spikes
latents
angle
0.00m
generated by teacher RNN
0,25m
1.75m
N
N
saddle
stable
-3
33
-3
3
5
generated by student RNN
N'
N'
-9
-3
3
0
100
14
0
100
-3
3
2,
timesteps
2
timesteps
2,
d
autocorr.
SD rates
mean rates
mean ISIs
example rate dist.
student teacher
22
student teacher
M
W
studemt
M
S
100 timelag
0.1 teacher
0.3
0
20
teacher
Figure 3: RNNs recover dynamics in teacher-student setups. a) Example ground truth latent trajectory and phase plane of low-rank RNN trained to oscillate (top left) and noisy observations of neuron activity (top right; 6/20 shown). A second low-rank RNN trained the activity of the first recovers ground truth dynamics. b) Same set-up, but with Poisson observations. c) The teacher network was trained task where it has to provide output corresponding an input cue. The student network, when given the same input during fitting, recovers the approximate ring attractor. d) Mean (+1SD) autocorrelation of the latents of the models from panel a, show the oscillation frequency is captured, as well as the decorrelation due to recurrent noise. The scale of the observed rates also agrees between student and teacher. e) Mean rates and ISI between student and teacher units of panel b match. f) Example rate distribution of one unit of the teacher and student RNN (of panel c), after onset of the 8 different stimuli.
different angles depending
3.2 Stochasticity allows recovering low-dimensional latents underlying EEG data
After validating our model on toy example, we went on several challenging world
EEG
generated
datasets. We first used an EEG dataset [2 43
channel name
with 64 channels containing one minute of con- tinuous data sampled at 160 Hz (Fig, 4). This
dataset was recently used in a study where gen- eralized teacher forcing (GTF) was used to fit deterministic RNNs with low-rank structure [8]. The GTF method obtains state-of-the-an results on several dynamical systems reconstruction tasks. It outperformed SINDy 44 neural differ- ential equations 455). Long- Eeee--pprees-eeeeeeee [66, and other methods, while using smaller latent dynamical system.
Fc1
5
time (s)
Figure 4: Example ground truth EEG 42 43 and (unconditionalll) generated traces by our model. Shown are 5/64 EBG channels.
Here we show that using a stochastic RNN with SMC instead of deterministic RNN with GTF, we can decrease the latent dimensionality even further, from 16 to just latents, while matching the original reconstruction accuracy (Table I.. We hypothesize this is because the data can be

---

well explained by stochastic transitions with simple underlying dynamics opposed to complex deterministic chaos.
Table 1: Lower dimensional latent dynamics than SOTA at same sample quality. We report median = median absolute deviation over 20 independent training runs, 'dim' refers the dimensionality of the model's s underlying dynamics and |0] denotes the total number of trainable parameters. Values for GTF taken from Hess et al. g
Dataset
Method
Dasp
Du
dim
[0]
GTF adaptive GTF g
2.1 + 0.2
0.11 0.01
16
17952
EEG (64d)
2.2
0.13 +.0+
16
17952
SMC (ours)
0.1
0.11
3920
We evaluated samples from our RNN with two measures which were used in previous work [8]. one KL divergence-based measure between the states (Dstsp), and measure over time, based on
the power spectra of generated and inferred dynamics (Dμ; see Supplement D.3.3). Unlike Hess et al. (8) who applied smoothing, we optimized also fit stochastic full-rank RNNs with variational SMC, however these models tend to have worse performance on this task, while also being less interpretable (Fig. $7).
models directly the raw EEG data. We
3.3 Interpretable latent dynamics underlying spikes recorded from rat hippocampus
a
b
spikes from rat hippocampus
generated spikes
mean rates (hz)
mean ISls (s)
neuron
gen train
10
10
1.0
1.0
lestigen testtrain
0.1
2
3
1.0 10
0.1
1.0
10
time (s)
test
test
pairwise corr.
d
latents and LFP
e
psd
f
coherence w.r.t. LFP
0.07
2,
gen train
22
0.1
0.2
posterior gen. Istents
Zg
0.01
0.1
LFP
Mww
-0.03
-0.03
0.07
2 time (s)
3
8 15 frequency (hz)
15 frequency (hz)
8
test
Figure 5: RNNs reproduce the stationary distribution of spiking data. a) We fit a rank-3 RNN to spike data recorded from rat hippocampus [47 48 (left), and generate new samples from the RNN (right). b) Single neuron statistics. Mean rates and means of interspike interval (ISI) distributions of long trajectory of data generated by the RNN (gen) match those of held-out of data (test). As reference we additionally computed the same statistics between the train and test set. c) Population level statisties. We plot the pairwise correlations between all neurons for generated data against the pairwise correlations in the test data. d) The corresponding latents generated by running the RNN look visually similar to the local field potential (LFP). e) The peak in the power spectrum matches between latents and LFP. f) The posterior latents show coherence with the LFP. As reference, we compute the coherence between the LFP and the latents generated by the RNN.
We next investigated how well our model can capture the distribution of non non-continuous time series. In particular, we used publicly available electrophysiological recordings from the hippocampus of rats running to drops of water or pieces of food [47 48). We binned the spiking data into 10ms bins and fit rank-3 RNN to ~50 of data. Samples generated by running the fit RNN autonomously closely matched the statistics of the recordings (Fig. 5 c). Previous investigations into this dataset have examined the relationship between spikes and theta (5-10 Hz) oscillations in the local field potential ( (177), and found that units were locked to the LFP rhythm, with the relative phase depending on the
7

---

subregions from which the units were recorded. The latents generated by the RNN are visually similar the average local field potential (Fig. ( and match power speetrum (Fig. Dp. While the model was solely trained on the spikes, the posterior latents (Eq. have clear phase relationship with the LFP, as evidenced by a high coherence between the posterior latents and L.P. contrast, and as expected, latents from running the RNN are not correlated with the LFP (Fig- 5D. The correspondence between generated latents and LFP was absent when we use related method for fitting RNNs (with deterministic transitions) to neural data (LFADS 77: Fig. $8 Using the bootstrap proposal also led to lower-quality samples (Fig. S
Units in rat hippocampus have been shown to code for position, e.g., through place cells [49], which tend fire if the animal specific location. To further investigate how well we can model recordings from the hippocampus, we fit rank-4 RNN to ad- ditional set of recordings of rats running on linear track 520 (Fig. STOp.
location
M
Zhou and Wei [S spikes recorded while the rat is moving, which we bin into 25 ms bins. The RNN again accurately reconstructs the distribution spikes and again has latent oscillations. Here the frequency at which power peaks is slightly higher than that of the LFP, potentially re- lated to phase precession (154)). While solely trained spikes, the posterior latents also allowed us to predict the position of the rats with reasonable accuracy (R² = 0.79 + 0.05 mean = SD, N 4 RNNs; Fig. 6 We also fit rank- 12 RNNs to around 15 minutes recording (again with 25 ms bins), which includes long intermediate periods where the rat is stationary. Here our generative model learns have higher theta power during running bouts, in line with the
we first focus only on the
predicted
30
time (s)
Figure 6: Posterior latents of our model (fit solely spikes) can be used to predict rat posi- tion.
data (Fig. $11
3.4 Extracting stimulus-conditioned dynamics in monkey reaching task
pre-movement
movement
decoded reaches
mean firing rate dissimilarity
true generated
posterior
y posdidoo
condition
distance
NS'
0.15
14
0.00
Zz
Z,
mean ISI dissimilarity true generated
generated
1.0
y imositen
condition
0.5
1,
2,
Z6
x position
condition
condition
Figure 7: Inferred and generated dynamics from the model fit to macaque spiking activity during reaching task. a) Latent states inferred from the macaque spiking data prior to movement initiation ("pre-movement') and during movement execution ('movement') ), colored by the intended reach target. b) Reach trajectories decoded from model-inferred neural activity. e) Dissimilarity matrices computed across the seven conditions (i.e., the seven colors in b) for per-neuron mean firing rate and ISI. We generate neural activity from the model by providing the same conditioning stimuli as in the real data. Then, for each statistic, compute and show the correlation distance between conditions in the real data (left) and model-generated data (right). e) Same as a, b, but with latent activity and behavioral predictions generated from the model with conditioning inputs including directions not seen in the real data (e.g., lime green). For clarity, we show only a subset of conditions in the decoded reaches.

---

We further investigated how well we can recover stimulus-condittoonss dynamics. We applied our method to spiking activity recorded from the motor and premotor cortices of macaque performing a delayed reaching task. This type of data has been popular for investigating neural dynamics underlying the control of movement 2 3 and evaluating neuroscientific latent variable models (755 (56. We first validated the ability of our method to obtain a sensible posterior by evaluating it on the Neural Latents Benchmark [56) (Supplement B.5 Table $2
We then went on to a set-up where we explicitly conditioned our model on external context. For simplicity, we constrained our experiment trials with straight reach trajectories in the data. We fit rank-5 model to these data while conditioning the RNN dynamics on the target position by providing the target position as input. Our model was able to infer single-trial latent dynamics and neuron firing rates that predict reach velocity with high accuracy at lower latent dimensionalities than models without inputs (Fig. 7p. R² 0.90 for this model, see Table $ for additional statistics).
We examined the posterior latents inferred by the model and found that our model recovers structured and interpretable latent dynamics. Before movement onset, latent states corresponded to the intended reach targets, which were near the edges of a rectangular screen (Fig. 7a. left), in line with [55]. During the movement period, the latents followed parallel curved trajectories that preserve target information (Fig. Za. right) and can be decoded to predict monkey reach behavior (Fig. 7b).
We then generated neural data from the RNN conditioned on stimulus input. Again, the distribution of spikes is well-captured (Fig. $T2). We additionally evaluated whether the model faithfully captures differences in spiking statistics across the seven reach directions, finding reasonable correspondence in dissimilarities between conditions in the generated and the real data (Fig-Ze) Finally, we simulated trained RNN with conditioning inputs, including reach directions not present in the data, and found that the structured latent space recovered by the model enables realistic generalization to unseen reach conditions (Fig. 7 e, lime green condition).
3.5 Searching for fixed points
In Proposition 1, we derived a bound on the number of systems of equations one has solve in order find all fixed points in piecewise-lineat low-rank RNNs. Re- cently, approximate algorithm for finding fixed points in piecewise- networks was proposed [25]. Here, we perform an exploration into how this compares to our an- alytic method by searching for fixed points of the RNN in Fig. 33 (top). For the same number of matrix inverses com- puted by our analytic method, the approximate method generally does not find all 17 fixed points (Fig. (). We note, however, that (unlike ours) the convergence of the approximate method depends on the dynamics of the RNN, and as result, there are theoretical scenarios where the approximate method can be shown to be faster. Yet we empirically also found scenarios where the approximate methods failed to converge within the time-frame our experiments (Fig. $13).
# fixed points found
approximate combined analytic
10
1000
4000 7000
# inverses
Comparison of our ana- lytic method (star) and the approximate method proposed in Eisenmann et al. 25 (blue) for finding the fixed points of the teacher RNN Fig. 33. We can also use Proposition to constrain the search space of the approximate method (or- ange). Error bars denote the minimum and maximum amount of fixed points found over 20 independent runs of the algorithm.
Figure 8:
Our analytic method relies on the insight that only a subset of all linear subregions formed by the piecewise-linee -inear acti- vations can be reached in low-rank networks. For networks with moderate rank, the cost of searching through all of the subregions might still be too high. We can, however, hugely reduce the search space of the approximate method
(25) (from (D + 1) to > D (M) ). at upfront cost (Supplement B.7 orange line in Fig.

---

4 Discussion
Here we proposed to fit low-rank RNNs to neural data using variational sequential Monte Carlo. The resulting RNNs are generative models with tractable underlying dynamics, from which we can sample long, stable trajectories of realistic data. We validated our method on several teacher-student setups and demonstrated the effectiveness of our method on multiple challenging real-world examples, where we generally needed latent dynamical system with very few dimensions to accurately model the data. Besides our empirical results, we obtained a theoretical bound on the cost of finding fixed points for RNNs with piecewise-lineal activation functions when they are also low-rank.
Adding stochastic transitions to low-rank RNNs can potentially hugely reduce the rank required to accurately model observed data, as demonstrated here with network fit to EEG data where we could reduce the dimensionality from 16 to just 3. While many methods that fit RNNs to neural data (e.g., 6 12) assume deterministic transitions, there is a rich literature concentrating on probabilistic sequence models in neuroscience (e.g., (2)) In particular, recent work termed FINDR [I] uses variational inference (but not SMC), to similarly find very low-dimensional dynamical systems underlying neural data. These stochastic dynamical systems were parameterized using neural differential equations [55. While Eq. [ can be seen as a neural differential equation with one hidden layer, our particular formulation allows us to find fixed-points effectively and map back to regular, mechanistically interpretable RNN (Eq. T after fitting, which enables additional investigations into neural population dynamics [8 20 22
We here - similar to FINDR (and (57]) did not use the adjoint method as is typical in the neural differential equation literature, but rather simple Buler-Maruyama discretisation scheme and standard backpropagation through time. However, one could investigate how we can integrate our approach with variational approaches that use adjoint methods when fitting latent neural SDEs 58 59 as well as with filtering approaches for continuous time systems [60]. This could be especially relevant for irregularly sampled time-series.
The reason we can do the mapping between a low-rank RNN (Eq. 1 and a latent dynamical system (Eq.2) 2 crucially relies on our assumption that samples from the recurrent noise process are correlated, such that they lie within the column-space of M. Valente et al. [61] showed that for linear low- rank RNNs arbitrary covariances in the full N dimensional space can be used, when increasing the dimensionality of the latent dynamics to twice the rank R (to the column space of both M and N), this however does not generalise to our non-lineat setting. We do expect correlated recurrent noise to be appropriate for modeling stochasticity arising from unobserved inputs or from partial observations [6]] additionally, correlated noise constituted a pragmatic choice that allows building an stochastic model that can allow for trial-by-trial variability while maintaining the tractability of low-rank deterministic RNNs.
Still, future work can investigate training networks with more relaxed assumptions on the recurrent noise models, including extensions non-Gaussian noise-processes. The latter could be of particular interest if more biologically plausible (i.e., spiking) neurons were used in the recurrence [3 62]
Our results also open up further avenues to explore questions in neuroscience. The relation between LFP and spike (phase) in the hippocampus has been of great interest [47] [54] 663, (64. While we performed some preliminary investigation into the relation between the inferred latents and the local field potential, further studies could perform a systematic investigation into their relation, for instance, by using multi-modal setup [13], to investigate multi-region temporal relationships and interactions [0]]
Taken together, by inferring low-rank RNNs with variational SMC, we obtained generative models of neural data whose trajectories match observed variability, and whose underlying latent dynamics are tractable.
Code availability
Code to reproduce our results is available at httt::://gittuu...cm/mmallllllll
10

---

Acknowledgments
This work was supported by the German Research Foundation (DFG) through Germany's Excel- Ience Strategy (EXC-Number 2064/1, PN 390727645), SFB 1089 (PN 227953431) and SPP2041 (PN 34721065), the German Federal Ministry of Education and Research (Tübingen Al Center, FKZ: OIIS 1803;; DeepHumanVissoon FKZ: 03 100199B)) and the European Union (ERC, Deep- ComechTome, 101089288), the "Certification and Foundations of Safe Machine Learning Systems in Healthcare" project funded by the Carl Zeiss Foundation. MP and MG members of the In- ternational Max Planck Research School for Intelligent Systems (IMPRS-IS). We thank Cornelius Schröder for feedback on the manuscript, and all members of Mackelab for discussions throughout the project.
References
[1] M. M Churchland, B. M Yu, M. Sahani, and K. V Shenoy. Techniques for extracting single-trial activity patterns from large-scale neural recordings. Current Opinion Neurobiology, 17(5):609-618, 2007.
[2] K. Shenoy, M Sahani, and Churchland. Cortical control arm movements: A dynamical systems perspective. Annual Review of Neuroscience, 36(1):337 (00993339999900 2013.
[3] Gallego, M Perich, L Miller, and S Solla. Neural manifolds for the control of movement. Neuron, 94: 978-984. 2017.
[4] S Vyas, M. D Golub, D Sussillo, and K. V Shenoy. Computation through neural population dynamics. Annual Review of Neuroscience, 2341)224 2020.
[5] D. Barack and W Krakauer. Two views on the cognitive brain. Nature Reviews Neuroscience, 22(6): 359-371, 2021.
[6] D Sussillo and L. F Abbolt. Generating coherent patterns of activity from chaotic neural networks. Neuron, 63:544-557 2009.
[7] C Pandarinath, D. J OO'Shea, J Collins, R Jozefowicz, S. D Stavisky, J. C Kao, F. M Trautmann, M. T Kaufman, S. Ryu, L. R Hochberg, M Henderson, K. V Shenoy, I., F Abbott, and Sussillo. Inferring single-trial neural population dynamics using sequential auto-encoders. Nature Methods, 15(10)8805 815, 2018.
[8] F Hess, Z Monfared, M Brenner, and D Durstewitz. Generalized teacher forcing for learning chaotic dynamics. In Proceedings of the 40th International Conference Muchine Learning, ICML'23, 2023.
[9] D Durstewitz. state space approach for piecewise-liness recurrent neural networks for identifying computational dynamics from neural measurements. PLOS Computational Biology. 13(6):1-33, 2017.
[10] Perich, Arlt, S Soares. M. E Young. P Mosher, J Minxha, Carter, Rutishauser, P. Rudebeck, C. D Harvey, and K Rajan. Inferring brain-wide interactions using data-constrained recurrent neural network models. bioRxiv:2020.1 12. 18, 8,342338, 2021.
[11] A Valente, J. W Pillow, and S Ostojic. Extracting computational mechanisms from neural data using low-rank mns.. In Advances in Neural Information Processing Systems, volume 35, 2022.
[12] F Dine, Shai, M Schnitzer, and H Tanaka. CORNN: Convex optimization of recurrent neural networks for rapid inference neural dynamics. In Thirty-sevnnth Conference on Neural Information Processing Systems, 2023.
[13] M Brenner, F Hess, G Koppe, and D Durstewitz. Integrating multimodal data for joint generative modeling of complex dynamics. In forst International Conference Machine Learning, 2024.
[14] O Barak. Recurrent neural networks versatile tools of neuroscience research. Current Opinion in Neurobiology. 46:1-6. 2017. Computational Neuroscience.
[15] D Sussillo and Barak. Opening the Black Box: Low-Dimensional Dynamics in High- High-Dimensional Recurrent Neural Networks. Neural Computation, 25(3):626669, 2013.
[16] H. S Seung. How the brain keeps the eyes still. Proceedings of the National Academy of Sciences, 93(23): 13339 13344, 1996.
11

---

[17] C E'liasmith and C. HI Anderson. Neural Engineering Computational Neuroscience Series); Computational, Representation, and Dynamics in Neurobiological Systems. MIT Press, Cambridge, MA, USA, 2002.
[18] F Mastrogiuseppe and S Ostojic. Linking connectivity, dynamics, and computations in lon-rank recurrent neural networks. Neuron, 99(3):609 623.e29, 2018.
[19] F Schuessler, Dubreuil, Mastrogiuseppe, Ostojic, and Barak. Dynamics of random recurrent networks with correlated low-rank structure. Physical Review Research, 2(1):013111, 2020.
[20] M Beiran, A Dubreuil, A Valente, Mastrogiuseppe, and Ostojic. Shaping Dynamics With Multiple Populations in Low-Rank Recurrent Networks. Neural Computation, 33(6):1572-1611 2021.
[21] A Dubreuil, A Valente, M Beiran, Mastrogiussppee and S Ostojic. The role of population structure in computations through neural dynamics. Nature Neuroscience, 25(6):783-794 2022.
[22] M Pals. J. H Macke, and O Barak. Trained recurrent neural networks develop phase-locked limit cycles in working memory task. PLOS Computational Biology, 20(2):1 22, 2024.
[23] C Curto, J Geneson, and K Morrison. Fixed Points of Competitive Threshold- Networks. Neural Computation, 31(5):94-155, 2019.
[24] M Brenner, F Hess, J. M Mikhaeil, I.. Bereska, Z. Monfared, P.-C Kuo, and D Durstewitz. Tractable dendritic RNNs for reconstructing noulinear dynamical systems. Proceedings of the 39th International Conference on Machine Learning, volume 162 of Proceedings of Machine Learning Research, 2022.
[25] L Eisenmann, Z Monfared, N Göring, and D Durstewitz. Bifurcations and loss jumps in mn training. In Advances in Neural Information Processing Systems, volume 36, 2023.
[26] K Morrison, A Degeratu, V Itskov, and C Curto. Diversity of emergent dynamics in competitive threshold- linear networks. SIAM Journal on Applied Dynamical Systems, 23(1 23(1):855 (:855)8884 2024.
[27] J. P Cunningham and B. M Yu. Dimensionality reduction for large-scale neural recordings. Nature Neuroxcience, 17(1 (7(11))15001159, 2014.
[28] B Petreska. B. M Yu, P Cunningham, G Santhanam, S Ryu, K. V Shenoy, and M Sahani. Dynamical segmentation of single trials from population neural data. Advances Neural Information Processing Systems, volume 24, 2011.
[29] J. H Macke, L Buesing, P Cunningham, B. M Yu, K. V Shenoy, and M Sahani. Empirical models of spiking in neural populations. In Advances in Neural Information Processing Systems, volume 24, 2011.
[30] S I Linderman, M Johnson, A Miller, R Adams, D Blei, and I. Paninski. Bayesian learning and inference in recurrent switching linear dynamical systems. In Proceedings of the 20th International Conference Artificial Intelligence and Statistics, volume 54 of Proceedings of Machine Learning Revearch, pages 914-922. 2017
[31] T. D Kim. Luo, Can, Krishnamurthy, J. Pillow, and C. D Brody. Flow-field inference from neural data using deep recurrent networks. bioRxiv:2023.1 11.14.6671.66 2023.
[32] Y Zhao, Nassar, Jordan, M Bugallo, and I Park. Streaming variational monte carlo. IEEE Transactions on Pattern Analysis & Machine Intelligence, (1000)):110-1100 2023.
[33] T. Le, M Igl. T Rainforth, T Jin, and Wood. Auto-encoding sequential monte carlo. In International Conference on Learning Representations, 2018.
[34] C. J Maddison, J Lawson, G Tucker, N Heess, M Norouzi. A Mnih, A Doucet, and Y Teh. Filtering variational objectives. In Advances in Neural Information Processing Systems, volume 30, 2017.
[35] C Naesseth, S Linderman, R Ranganath, and D Blei. Variational sequential monte carlo. In Proceed- ings of the Twenty-Firs First International Conference on Artificial Intelligence and Statistics, volume 84 of Proceedings of Machine Learning Research, 2018.
[36] C Sourmpis, C Petersen, W Gerstner, and Bellec. Trial matching: capturing variability with data- constrained spiking neural networks. In Advances in Neural Information Processing Systems, volume 36, 2023.
[37] A Doucet and A. M Johansen. A tutorial on particle filtering and smoothing: Fifteen years later. The Oxford Handbook of Nonlinear Filtering, pages 656-704, 2011.
12

---

[38] D. P Kingma and M Welling. Auto-Encoding Variational Bayes, In 2nd International Conference on Learning Reprexentations. ICLR, Conference Track Proceedings, 2014.
[39] K Doya. Bifureations of recurrent neural networks in gradient descent learning. IEEE Transactions on Neural Networks, 1993.
[40] Zenn and R Bamler. Resampling gradients vanish in differentiable sequential monte carlo samplers. In The First Tiny Papers Track at ICLR 2023, Tiny Papers @ ICLR 2023, 2023.
[41] T Zaslavsky. Facing up to arrangements: face-count formulas for partitions space by hyperplanes. Memoirs of American Mathematical Society, 154:1-95, 1975.
[42] G Schalk, D McFarland, T Hinterberger, N Birbaumer, and J Wolpaw. Bci2000: general-purpose brain-computer interface (bci) system. IEEE Transactions on Biomedical Engineering. 1043, 2004.
[43] G Moody, R Mark, and A Goldberger. Physionet: a research resource for studies of complex physiologic and biomedical signals. Computers in cardiology, 27:179-88 82, 2000.
[44] S. L Brunton, J. L Proctor, and J. N Kutz. Discovering governing equations from data by sparse iden- tification of nonlincar dynamical systems. Proceedings of the National Academy of Science, 113(15); 3932-3937, 2016.
[45] R. T. Q Chen, Y Rubanova, J Bettencourt, and D. K Duvenaud. Neural ordinary differential equations. In Advances in Neural Information Processing Systems, volume 31, 2018.
[46] T. Rusch, Mishra. N. Erichson, and Mahoney. Long expressive memory for sequence modeling. In International Conference on Learning Representations, 2022.
[47] K Mizuseki, A Sirota, E Pastalkova, and G Buzsáki. Theta oscillations provide temporal windows for local circuit computation in the autoriinal-hipoocampee loop. Neuron, 64(2):267-280, 2009.
[48] K Mizuseki, Sirota, Pastalkova, and Buzsáki. Multi-unit recordings from the rat hippocampus made during open field foraging. Database: CRCNS, 2009.
[49] O'Keefe. Place units in the hippocampus of the freely moving rat. Experimental Neurology, 51(1): 78-109, 1976.
[50] A. D Grosmark and G Buzsáki. Diversity in neural firing dynamics supports both rigid and learned hippocampal sequences. Science, 351 (6280)))): 1400) 1440-1443, 2016.
[51] Z Chen, A. D Grosmark, H Penagos, and M. A Wilson. Uncovering representations of associaadd hippocampal ensemble spike activity. Scientific Reports, 6, 2016.
[52] I.. J Grosmark, A.D. and Buzsáki. Recordings from hippocampal area cal, pre, during and post novel spatial learning. Database: CRCNS, 2016.
[53] D Zhou and X.-X Wei. Learning identifiable and interpretable latent models of high-dimensional neural activity using pi-vac. In H Larochelle, M Ranzato, R Hadsell, M Balcan, and H Lin, editors, Advances in Neural Information Processing Systems, volume 33, pages 7234-7247. Curran Associates, Inc., 2020.
[54] O' Keefe and M. L. Recce. Phase relationship between hippocampal place units and the eeg theta rhythm. Hippocampus,, 333:3317-330, 1993.
[55] G Santhanam, B. M Yu, V Gilja, S. Ryu, A Afshar, M Sahani, and K. V Shenoy. Factor-analysis methods for pireeacperformmacc neural prostheses. Journal of Neurophysiology. 102(2)11315-1330, 2009.
[56] F Pei, Yc, D. M Zoltowski, A Wu, R. H Chowdhury, H Sohn, J. E O'Doherty, K. V Shenoy, M. Kaufman, M Churchland, M Jazayeri, L. E Miller, J Pillow, I. M Park, E. L Dyer, and C Pandarinath. Neural latents benchmark '21: Evaluating latent variable models of neural population activity. In Advances in Neural Information Processing Systems (NeurIPS), Track Datasets and Benchmarks, 2021.
[57] C Versteeg, A. R Sedler, McCart, and C Pandarinath. Expressive dynamics models with nonlinear injective readouts enable reliable recovery of latent features from neural activity. In Proceedings of the 2nd NewrIPS Workshop on Symmetry and Geometry in Neural Representations, volume 228 of Proceedings of Machine Learning Research, pages 255-278. PMLR, 16 Dec 2024.
13

---

[58] X Li. T.-K. L Wong, R. T. Q Chen, and D Duvenaud. Scalable gradients for stochastic differential equations. In Proceedings of the Twenty Third International Conference Artificial Intelligence and Statistics, volume 108 of Proceedings of Machine Learning Research, pages 3870-3882. PMLR, 26-28 Aug 2020.
[59] R Deng, M. A Brubaker, Mori, and A Lehrmann. Continuous latent process flows. In Advances in Neural Information Processing Systems, volume 34, pages 5162-5173, 2021.
[60] T Sottinen and Särkkä. Application Girsanov theorem to particle filtering discretely observed continuous time non non-linear systems. Bayesian Analysis, 3(3):555 584, 2008.
[61] A Valente, Ostojic, and J. Pillow. Probing the Relationship Between Latent Linear Dynamical Systems and Low-Rank Recurrent Neural Network Models. Neural Compulation, 34(9):1871-33 2022.
[62] Cimeša, L Ciric, and S Ostojic. Geometry of population activity spiking networks with low-rank structure. PLOS Computational Biology, 19(8):1334, 2023.
[63] G Buzsáki. Rhythms of the Brain. Oxford University Press, edition, 2006.
[64] S Liche, Niedick, M Pals, T. Reher, J Faber, J Bostroem, C. E Elger, J. H Macke, and F Mormann. Phase of firing does not reflect temporal order sequence memory of humans and recurrent neural networks. wtt.22222..9.2555097777 2022.
[65] Schläfli. Theorie der vielfachen Kontinuität. Birkhäuser Basel, Basel, 1901.
[66] R. C. Buck. Partition of space. The American Mathematical Monthly, 50(9):541-544, 1943.
[67] R Stanley, An introduction to hyperplane arrangements. Geometric Combinatorics, 13:389 496, 2007.
[68] M. R Keshtkaran, A. R Sedler, R. H Chowdhury, R Tandon, D Basrai, S. L Nguyen, H Sohn, M Jazayeri, L. E Miller, and C Pandarinath. A large-scale neural network training framework for generalized estimation of singlc-trial population dynamics. Nature Methods, 19(12):1572- -5777, 2022.
[69] B. M Yu, J. Cunningham, G Santhanam, S Ryu, K. V Shenoy, and M Sahani. Gaussian-proees process factor analysis for low-dimensional single-trial analysis of neural population activity. In D Koller, D) Schuurmans, Bengio, and I. Botton, editors, Advances Neural Information Processing Systems, volume 21, 2008.
[70] Ye and C Pandarinath. Representation learning neural population activity with neural data transformers. Neurons, Behavior, Data analysis, and Theory, 5(3), 2021.
[71] D Lawson, Raventós, Warrington, and S Linderman. Sixo: Smoothing inference with twisted objectives. In Advances in Neural Information Processing Systems, volume 35, pages 38844-38858 2022.
[72] A Paszke, S Gross, Massa, A Lerer, Bradbury, G Chanan, T Killeen, Lin, N Gimelshein, Antiga, A Desmaison, A Kopf, E Yang. 7. Vetioo, Raison, A Tejani, Chilamkurthy, B Steiner, I. Fang, Bai, and S Chintala. Pytorch: An imperative style, nigh-performance deep learning library. In Advances in Neural Information Processing Systems, volume 32, 2019.
[73] A Orvieto, S. L Smith, A Gu, A Fermando, C Gulcehre, R Pascanu, and S De. Resurrecting recurrent neural networks for long sequences. In Proceedings of the 40th International Conference Machine Learning. ICMIL'23, 2023.
[74] L Liu, H Jiang, P He. W Chen, X Liu, J Gao, and J Han. On the variance of the adaptive learning rate and beyond. In International Conference Learning Representations, 2020.
14

---

A Supplemental material
A.1 Proof of preposition
A.1.1 Problem definition
We are interested in finding all fixed points of the following equation:
dx
x() Jo(x()))
dt
with x(t) RR, element-wise nonlinearity 6(x,) = bdd max(x; h (dd) and low-rank matrix J MN with M, N E RN and Rt S N. Since only scales the speed of the dynamics, we will, for convenience and without loss of generality, assume J
A.1.2 Preliminaries: Fixed points in Piecewise-linear RNNS
First, we briefly repeat results from Durstewitz [9] Assume D = 1, o(x,) max(x, - h2). To find all fixed points of Eq. start by redefining by introducing diagonal indicator matrix:
dl
d2
Dg
(9)
dx
with d, = 1, x; > h;
0, otherwise
Then our RNN equation, for given and corresponding Dg reads:
dx
dt x(t) JDnx(() - JDnh.
Each of the 2 configurations of Da corresponds to region in which the dynamics are linear. Thus, for each configuration, we can solve:
- JD Dep_DDDDJDD.. - JDph,
x = (JD, - JDah.
Next, we check whether the obtained x* is consistent with the assumed Dg (Eq, 9). If so, we found a fixed point of the RNN. We have to check, as the solution to the system of linear equations can lie outside of the linear regions specified by Dg. Note that if for some Dg the matrix JDn 1 is not invertible, then there is no single fixed point, but we still can find structure of interest (e.g., direction with eigenvalue corresponds marginal stability, i.e., a line attractor).
A.2 Preliminaries: Fixed points in Piecewise-linear low-rank RNNs
First, assume x(0) is in the subspace spanned by the columns of M. With the low-rank assumption, we can rewrite Eq. g for all [0, oo), by projecting it M [18 20 21:
dz = 2(()) N. (M2(( h)
(10)
di
with x(t) Mz(t).
Now assume x() contains some part x=(()) not in the subspace spanned by M, i.c., we have x() Mz(0) (0). The dynamics of (t) are simply given by de = (1) which will decay to its stable point at 0 irrespective of z(t), and can thus not contribute additional fixed points.
15

---

0
1
D
M
x,=h,
0
]
D
D
x,=h,
Supplementary Figure 1: Proof sketch including Dn's. The phase-space of RNN with N (here 2) units with activation max(0, -h,)
h,) is partitioned into 2 (here 4) regions in which the dynamics are configuration of Dg. If dynamics
linear, each corresponding subspace spanned by the columms of M, only a subset (here 3) can be reached. Each unit intersects the space spanned by the columns of M with hyperplane (the pink points in the Figure). The amount of linear regions in M, thus becomes equivalent to "how many regions can we create in R-dimensional space with N hyperplanes?"
confined to the R- R-dimensional
Naively, using the same strategy as before to obtain all fixed points z, we would need solve 2 linear systems of It equations (again for all configurations of Dn):
x
(11)
A.3 Preliminaries: Hyperplane arrangements
In the subsequent section, will turn the question of how many equations need to solve find all possible fixed points. Recall that it is possible to calculate the fixed points analytically because piecewise-linee linear nonlinearities partition space into subregions in which dynamics are linear. Each of the linear regions corresponds to configuration of Dg. For networks with low-rant connectivity, we have to consider only a small subset of those, only a small subset of all configurations Dn correspond to x's within the column space of M (See Fig. S To find out exactly how many regions lie within the column space, we will need to answer the question:
how many regions can we divide R-dimensional space, with N hyperplanes? To answer this question in general, we will need theorem from the field of hyperplane arrangements
65 67 Here we give a brief introduction.
Introduction to hyperplane arrangements: A finite arrangements of hyperplanes is a set of {a....... in some vector space
N affine subspaces an R dimensional subspace defined by linear equation a; := {v V I Jv some m E h E R. Note that any linear system of equations Mv h with M E g x equivalents defines arrangement of N hyperplanes R dimensional space. In Fig. $Da.c, we show arrangements of 3 hyperplanes in IR2 In this case, hyperplane is a line, but there are infinitely many possibilities on how we can arrange these lines in -io-dimensional space. We are interested in
RR Recall a hyperplane is h for
N(AA == number of regions A partitions R.
where regions correspond to the connected components of IR A. this simple case, we can visually verify that the arrangements in Fig. $2h partitions the space into regions, whereas the arrangement in Fig. She partitions the space into only 6 regions. Clearly, the number regions A partitions space in strongly related to the number of unique intersections of lines. We have fewer regions in Fig. s simply because all lines intersect at the same point. If we can wiggle the hyperplanes little, and not change the number of regions (as can do in Fig. $2h. but not Fig. $2 we call the hyperplanes in general position (see Theorem [ for a formal definition).
To count the amount of regions for any arrangement of hyperplanes, we can leverage an algebraic construction called the intersections posel L(A). This the set of all nonempty intersections of hyperplanes in A and includes V. Elements of this set are generally referred to as flats. The flats
16

---

a
Q2 0 a3
01
d1
02
a3
R
1
2 d1 0gg nas
02
a3
01
a1
02
03
-1
IR2
Supplementary Figure 2: a) An arrangement of hyperplanes a1, d2 and d3 in general position. b) the associated intersection poset the arrangement a. e) An alternative arrangement with its associated intersection poset. d). Blue numbers indicate the value of the Möbius function.
x y in the intersection poset. We visualized example intersection posets of the previous examples (Fig. $2b, d). Here we organized the flats by dimensionality (such a visualization called Hesse Diagram). Importantly for any real arrangement A, A. (LA) solely depends on ((A) (Corollary 2.1, (67])
ordered by reverse inclusion S
To calculate (AA) from I.(A), we need one last construction, namely the Möbius function, recur- sively defined by
((X (- (- ((X, MX,
if X if C X.
(12)
8',
The numerical values for the example are shown in Fig. $
Theorem (Zaslavsky's Theorem;
67)).. Given vector space V RR and an arrangement of N hyperplanes A = fai...., aN) } on V. then the number of regions A partitions V in (denoted N (A),, can be expressed as follows
N(A)
Σ SEL(A)
I(R 8) ()(-))dim())
furthermore, holds that
M(A) S Σ N
(13)
with equality if and only if A is in general position i.e., A must satisfy
(i) {ai...., ap} C A and p I R + dim(Ni (_ a;) = N
(ii) ap} S A and
R i= = 0
One can verify this fact for the given example shown in Fig. S We refer to Stanley 67 for an in-depth formal introduction to this topic. Fundamentally. it is based on the following recursion that number of regions for any arrangement satisfies
N(AU (oN+1)) =N(AA) (AANN1)
17

---

where A := {aN +1 Najla, A, aN+1 Na: # 00 0,4111 & a (Lemma 2.1, Stanley [67 Note that AN++ itself an R dimensional vector space, and each intersection 4N+1 N a, an R dimensional hyperplane within aN+1 (e.g., the intersection of two planes is a line within the planes). Hence 1 itself an arrangement of N hyperplanes, but in an R dimensional subspace. In fact, the intersection poset exhaustively enumerates the elements of all possible A4t and the Möbius function can be shown to satisfy the above recursion.
h; 0) i.e D = 1, each neuron would partition space by a single hyperplane x; h, or equivalently the R dimensional subspace by the hyperplane M,z h;- Hence, the hyperplane arrangement is determined by the matrix M and offset h. As these quantities are learned during training of the RNN, this arrangement is often simply numerically unlikely that two hyperplanes are exactly parallel or intersect in exactly the same "point". This does, however, change in the general case 1, for which we derive a tighter bound in the section below.
If we choose (xx)) = max
x; m
general position because is
Arrangements of parallel families For the general case o(x;) -d== by max(x; h h (..0)) x; each neuron will partition space with D hyperplanes (,0)x;- b (d) 2 x; h (d) before; cquivalently each neuron partitions the dimensional subspace with D hyperplanes M, h (
.0)
h;
Notably, all the D hyperplanes here will share the same of M, and thus they parallel. Clearly, any such arrangement cannot be in general arrangement by definition.
The resulting arrangement will have very specific structure. Let's define
A; == {a.1...., a,D]
a family of D parallel hyperplanes. Any pair of hyperplanes au; Oim A, parallel. A low- rank RNN with N neurons and a general piecewise-linee activation function will thus lead to an arrangement consisting of N families of parallel hyperplanes.
We can use this specific structure to obtain a tighter bound.
Lemma 1. Let A A U... Ax-1 be an arrangement of N families of D parallel lines, then it satisfies the following recursion
D EA
N(AUAN)= =
Furthermore, denote by ( R. D) the maximum number of regions attainable by any arrangement of N families of D parallel hyperplanes in R dimensional space then
- N(N,R, (N(N-,,R) + DD N-N 1, R 1)
Proof. To add AN to A. we have add D new parallel hyperplanes. We can do so by iteratively applying Lemma 2.1 [67]. We obtain
N(AU (ANI, (aNI,..., ..n.)) = N(AU (aNI, (N(D-1))) + -N((UU )) aND) (
aNd
=((A)+ ΣN
AU
[uuUt HaNi)
Now note that A = {aNz O dim amm E A,aN; O dim + 0, UNj & am hence by definition only hyperplanes that intersect with aNj are included in this set. As aNz is parallel to any other ax; for all * j. all aNj O aNi cannot be in the set. Hence for any d, we have that
d-1 [aane]
GNd
( i=!
((AuN))
which proves the first equation.
Recall that we define (NN, R, D) as the maximum number of regions attainable by any arrangement. Notice that A by construction is an arrangement of N families of parallel hyperplanes in dimension. Thus by definition (AA) NNN - 1, R, D).
18

---

As noted before, the intersection set of two hyperplanes in dimension R is itself hyperplane of dimension 1. Furthermore the intersection sets of D parallel hyperplanes with ayd. remain parallel. Hence A is an arrangement of at most N - 1 families of D) parallel hyperplanes in R - 1 dimensions. Thus (A S M(N - 1, R 1, D) leaving us with
M MANU............
- <N(N - D, D -N(N - D-N(N - D,
As this holds for any arrangement, also holds for the arrangement that has (NN,, R, D) regions (i.e., which maximizes the number of regions) and, therefore, proves the second equation.
Lemma 2. Let A be arrangement of N families of D parallel hyperplanes. Then, it holds that
((A) S 1 D (7)
F= 0
with equality if each family is in a general position, i.e., that every subarrangement {a15 "NNN } for all 1 < J4 D general position.
Proof. We will first construct an intersection poset L(A) on the level of families A, in general position. After all, the intersection properties between these families is the same as between their elements, e.g., if an intersects 4₁1 then also all lines in A, intersect all lines in Aj.
The resulting intersection poset T.(AA be clustered into the corresponding families. We visualize the construction in Fig. $3
At each rank (level from bottom to top), we can choose exactly (N) families of hyperplanes that intersect (exactly the case if we just have N hyperplanes in general position). To obtain a flat of dimension R r we have to choose r out the N hyperplane families without replacement.
If, c.g., two families of parallel hyperplanes A,, A, intersect, then any element ag will intersect with any element 0g for all S k < D leading to at most D flats within each family (there can be less other families might intersect in the same "point"). In general, each cluster of intersections of families will contain at most D flats.
By construction of L(A) and Theorem I the lemma follows directly.
To show that this construction indeed an upper bound for all arrangements, we can Lemma I There, we established a recursion. which any such upper bound must satisfy. Hence, assume (,N = -H Dr (M). Notice that using Pascal's identity, we can rewrite
N
W(N,R, Σ D D
r
R
((N-) )
1)
N -
r- 1
+
P=0
Σ D"
R
(N-1)
N
+ Σ D"
R
- 1
r-
r=0
N -1
N -
N -1
Dr
Σ
D"
r=0
=
R-1
N - 1)
- WVNN-,,,D))
Σ Dr+1
-
- D.
19

---

N Ax {N akdy akdy,d,, ................................
Id.,..., de - D)
A Nassisl D
A n Av {au MaNylini D)
A 1 {an/i <= <= D}
An Haxili <= D}
R
Supplementary Figure 3: Construction of the intersection poset L(A) for an arrangement of N families A, of D parallel hyperplanes in "general position".
A.4 Proof of proposition
Using the previously derived techniques, we will prove here the main proposition. Furthermore, in Algorithm| L pseudo-code is given to compute all fixed points in practice.
Proposition 1. Assume the RNN of Eq. [ with J of rank R and piecewise- linear activations: ((x;) = - b(4) max(x; 1 h(()) 0). For fixed rank R and fixed number of basis functions D,, we can find all fixed points in the absence of noise, that is all x for which dx = 0, by solving at most dt
R) linear systems of R equations (for fixed R).
Proof. By definition, each neuron partitions R- in D+1 linear regions with D hyperplanes described by x = h (d) for the i'th neuron. Using that in the columnspace of M, we have x = Mz, follows that each neuron partitions the R dimensional subspace spanned by columns of M, with D hyperplanes described by Σ M2 = h (4) Notice that these hyperplanes are parallel, as they all share the same coefficients M, but have different offset h (()). Using Lemma we know that there can only be Σ/= Dr ( such regions.
How do we find those regions? Let's 's first consider the case of D = 1, and assume that the hyperplanes in general position. We can find the corresponding configurations of Da as follows. We first obtain the set of all intersections of R hyperplanes. For this we try solve (N) systems of R equations. Let MR E IR Rx be the matrix obtained by choosing different rows I,..., R of M E NR (i.e., picking It neurons), then we may find the corresponding intersection of R hyperplanes by solving the following linear system of R equations
On MR hR
and
Xn
M
which will always have unique solution if all hyperplanes are in general position, as then all M have rank R. Each Xn has 2 possible bordering linear regions. We can find the corresponding matrices of each of those subsections as follows. First d, = T(xn < 0)
Da diag([4 (d), .., [dl....,dd for all i <= N. By construction 1,..., from it dR can become either zero one, depending in which region we and up. Hence, the 2 regions correspond to one which either combination of neurons I,... R is active (meaning that it is above the threshold). We thus just have check all combinations d1,,. d...... dn E {0, 1)
Xn will be exactly at the threshold, by moving away
Using this. we will find at most D ER (A) unique configurations (as this is the maximal number of regions possible for D = 1). To find all the fixed points we hence have to solve Eq. T1 for each configuration. We thus end up with solving (N) systems of R linear equations to find all regions, and another E (*) E R) (for fixed R) systems of R linear equations to find all fixed points.
R
Let us now consider the case for D > 1. Note that an RNN with N units and D basis functions per unit, can be expanded an RNN with ND units with activation ((x;) max(x; ; h,, 0)
Theorem 1). Any fixed point can then still be analytically computed using Eq. 11 We expand
the network but keep track of all Σ R D Dr(N) possible intersections. It still holds that from each
20

---

intersection, we can reach 2 regions. In total, we will now find at most Σr= (r regions (Lemma 2 To find all the fixed points, we hence have to solve (+) Dr + +E (+) systems of R linear equations, which for constant D and R has a cost ((N))
Finally, let's consider the case when hyperplanes are not in general position (which is unlikely to happen when doing numerical optimization). If there are intersections of more than hyperplanes, proceed as before, but in case the intersection of It hyperplanes intersects additional hyperplanes, set the diagonal elements of Dp corresponding to these additional hyperplanes arbitrarily to I (as intersections including the additional hyperplanes are considered separately). On the other hand, in case some hyperplanes are only part of intersections of less than R hyperplanes (because they became parallel), proceed follows. Instead of considering only intersections R. hyperplanes, we now also consider all possible intersections r hyperplanes, with 1 r r S R. For this, we solve no more than Σ,' (?) systems of equations. Let M, E Rx*x be the matrix obtained by choosing r different linearly independent rows 1,... of M RNNRR, then we may find point on the corresponding intersection of hyperplanes (note that the intersection itself can now also be hyperplane) by to solving the following linear system of r equations
currently considering
nn M
and
X-
Mzn
with being the pseudoinverse. We here now end up with solving no more than Σ R (*) systems of linear equations find all regions, which has equal cost in N as the previous cases.
We here provide pseudocode. For simplicity, we restrict ourselved to the case of = and assume that the arrangement specified by M and h is in general position. This be generalized to the general setting as presented in the proof.
Algorithm 1: Improved exhaustive search for all fixed points
Data: E Nx M NxR E RV
Result: set set of all fixpoints, D_set the set of all relevant Dg configurations.
D_set }:
ses = U
idr [1,..... N]:
Il Find feasible configurations idr comb all (n) combinations of indices idæ;
for (i ... (i,....,ig) in ded__comb do
MR (M,
hg hh((i....iiz)))) ix)];
// M is invertible as the arrangement is in general position
Zn solve(MR. hg);
X- Mzni
d_inil Xn h; for (Ul, (Ul,..., UR) {0, 1) do
d = init (i),-- sete,,,,.
a diag(d);
D set = D set U (D));
end
end
II Find fixed points, for the at most ER (*) configurations for Da in Det do
solve(N TR,NNDDDM--
x* = Mz':
I/ check if fixed point is consistent with assumed Do
if diag(x" > h) = Dg then
set set }z
end
end
21

---

B
Additional figures & tables
B.1 Additional statisties for Teacher-Student setups
b
pairwise correlations
student
0.2
0.0
-1
-1
-1
0 teacher
0.00
0.25
-1
Supplementary Figure 4: a-c) Pairwise correlations between units of the modes for panel a-c) of respectively. Note that e is computed over all conditions.
Fig. 3
k=64
student
22
k-1
Boctstrap (k-64)
.05
Isacher
Supplementary Figure 5: Our method allows recovering the true latent noise in student-teacher setups. We repeated the experiment of Fig. 3a for teacher networks with three levels latent noise, with diagonal covariances matrices Σg II. For each teacher trained 5 student networks with varying number of particles (k), as well as using the bootstrap proposal (i.e., sampling from the prior). The standard deviations of the latent noise process only matches between the student and teacher, we use enough particles during training. The bootstrap proposal (with k=64) is not as reliable as the optimal proposal.
inputs
teacher
student
.5
5
-5
-55
0
50
100
10
0
10
10
0
10
timesteps
2
Supplementary Figure 6: To demonstrate that our method works with time-varying input, rank-1 teacher RNN with 60 units trained to report the sign of time-varying stimulus (left). We then generated 400 trials of data with observation noise covariance Ex .01I and latent noise covariance Ex .00251. We trained student the observed activity of the teacher for 400 epochs. The matching latent dynamics of the student and teacher lie in the column space of the recurrent and input weights (right; coordinates and respectively; see Supplement
22

---

B.2 EEG: Inferring full-rank RNNS
Hank (Nmb112)
Dtep
D22
Rank (N=30)
Full Rank (N-128)
2
Supplementary Figure 7: We fit full-rank RNNs to EEG data, by parameterising the mean of the transition distribution as F(z,) = az, + (1 aaJJaa (Zz). We trained full-rank RNNs with units (a roughly similar amount of parameters as our rank RNNs with 512 units), as well as full-rank RNNs with 128 units (over 10 times more parameters). The KL. divergence based measure (Dstsp). between generated samples and data is worse for the full-rank RNNs, while also being less interpretable.
factors
psd
generator
factors
psd
2,
wMM
LFP
sampled Inferred
Zg
Zg
mmm
.01
.01
2g
wmm
Zg
15 frequency (hz)
15 frequency (hz)
time (s)
time (s)
time (s)
Supplementary Figure Latent Factor Analysis Dynamical Systems (LFADS) is current method that can be used to fit RNNs neural data, and while is generally used for inference, also be used for generation. We here explored sampling from LFADS, both using the autonomous version, and when using the controller (i.c., stochastic inputs). We fitted LFADS the HPC-2 dataset (using AutoLFADS for model selection [68]). a) In the case of an autonomous LEADS model, one samples initial condition from the prior and then simulates deterministic RNN forward. As long sequences not variability can be explained by variability in the initial condition, the latents end up not representing any variability that resembles the underlying system (cf. Fig. 5). b) In the case of a full LFADS model with the controller, one can sample both an initial condition and time-varying inputs from the controller's auto-regressive prior. Here the full model seemed rely overly on the controller's s data-inferred inputs for inference, which deviated quite strongly from the samples from the controller's rugo-regressive prior. As a consequence, the generated latents do not seem to represent variability that is meaningful.
B.3 HPC-2, additional results
LFP
spike stats
9
kv64
7
D rT.5
k-1
bootstrap (lmB4))
trainitest
3
0
Maae
Musi
pW corr.
Supplementary Figure 9: The Hellinger distance (D#) between the power spectrum of latents and LFP of HPC-2 is lower when using the bootstrap proposal or too few particles (left), and simulated data is slightly worse when using the bootstrap proposal (right)
23

---

B.4 HPC-11, additional results
b
data
generated spikes
mean rates (hz)
CV ISIs
10
gen. train
nouro
lest/gen testitrain
10 time (s)
10
test
test
e
pairwise corr.
latents
latents and LFP
psd
0.2
2,
LFP
wwmm
gen train
Z2
Z
0.1
Zz
Zg
Zg
0.001
Z4
Z4
-0.1
-0.1
test 0.2
10 time
20
0
0.2
10.0 20.0
time (s)
frequency (hz)
Supplementary Figure 10: a) We fit a rank RNN spikes recorded from hippocampus 150-52)
and generate new samples from the RNN (right), taking only the part of the recording where the rat running. b) Single neuron statistics. The mean rates and coefficient of variations of interspike interval (ISI) distributions of long trajectory of data generated by the RNN (gen) match those of held-out set of data (test). As a reference we additionally computed the same statistics between the train and test set. c) Population level statistics. The pairwise correlations between neurons for generated data and the test data. d) The corresponding latents generated by the RNN consists of 10Hz top of slower oscillations. e) Latents with further zooming in (on time),
(fast theta) oscillations shown together with the LFP signal. f) The power spectrum of latents sampled from the RNN, which show power at slightly higher frequency than that of the LFP [54
b
mean rates (hz)
ISIs
pairwise corr.
latents
0.2
gen train
testlgen
teatttrai
-0.1
0
4
-0.1
0.2
10 time (s)
20
test
test
test
d
predicted location
generated location
psd
predicted
model running model static LFP rat running
0.1
rat
LFP rat static
0.001
200
400
200 time (s)
400
0 2
20.0 frequency (hz)
10.0
time (s)
Supplementary Figure 11: We additionally fit rank- 12 RNN to whole recording (2067 seconds resampled to 40 Hz) which includes long bouts where the rat is stationary. a) We generate new sample from the RNN and obtain matching spike statistics. b) Generated latents by our model. e) Inferred posterior latents can again be used to predict the location of the rat on a held-out set. d) Using the same decoder generated latents, of stationarity and running. e) The mean power of the latents at theta frequency during running bouts higher then during stationarity bouts. The increased theta power during running is also there in the LFP data again with latent dynamics obtained from spiking data oscillating at a slightly higher frequency than the LFP.
obtain model that also predicts alternating bouts
24

---

mean rates
SD rates
mean ISI
SD ISI
pairwise corr.
angle 135" gen/train
lest/frain
0.1
10
0.5
0.00
0.25
angle = 90 gen/train
0.1
10
0
50
0.0
0.5
0.00
0.25
angle = 45" gen/train
0.1
0
0.0
0.5
0.00
0.25
angle = 0 genitraie
0.1
1
10
0
50
0.0
0.5
0.00
0.25
angle = 45 genttrain
0.1
10
0
50
0.0
0.5
0.00
0.25
angle= 135
gen/train
0.1 1 10
0
50
0.0
0.5
0.00
0.25
angle = 180'
gen/train
0.1
10
0
50
0.0
0.5
0.00
0.25 test
test
test
test
test
Supplementary Figure 12: Spiking statistics of model-generated (teal) and train data (brick red) compared against test data.
25

---

B.5 Neural Latents Benchmark evaluation
We applied our method to the MC_ MC_Maze dataset of the Neural Latents Benchmark (NLB) [56] at 20 millisecond bin size (Table 2) by using our method to obtain expected Poisson rates, given a filtering posterior over latents. The benchmark evaluates methods smoothing bits-pere bits-per-spike) assesses the quality of firing rate predictions for a set of held-out neurons that are unobserved in the test data, evaluated with the Poisson likelihood of the true spiking activity given the rate predictions. 'vel R2' evaluates how well the model's inferred firing rates can predict the subject's s hand velocity. 'PSTH R2' evaluates how well peri-stimulus time histograms (PSTHs) computed from model-inferred rates match empirical PSTHs from the data. 'fp-bps' evaluates predictions on heldout timesteps (which we predict by running the RNN forward from the last data-inferred step). We found that our method outperforms classical methods (GPFA [69] and SLDS [30] while certain state-of-the-art deep learning (LFADS 77 688 Neural Data Transformer (70) are slightly better than our method on the 'co-bps" metric, but our method matches them in the 'vel R2' metric (in case we include smoothing information in the proposal). We do note that NLB metrics center around evaluating the quality of smooth rates inferred from spikes, which is not the central focus of our method, which is generation, i.e., sampling noisy trajectories that reproduce variability in the data. We here found that the quality of inference increased when using an non-causal CNN encoder as part of the proposal distribution, and additional gains might be obtained by also changing the target distribution to a smoothing (instead of a filtering) one [71].
number of metrics: 'co-bps' (co-
While our method also has comparatively lower dimensionality than the other deep learning ap- proaches, a latent dimensionality of 36 is still considerably higher than all networks considered in the Main text. We reason that we need high number of latents, because the full MC_ MC_Maze dataset has a large number of conditions (108), spanning multiple conzgugoiiuuttionn,, which may be difficult to fully model with autonomous oow-dimessional latent dynamics.
Table 2: Performance of our method on the MC SC_Maze dataset of the Neural Latents Benchmark, 'dim' refers to the dimensionality of the model's underlying dynamics (where possible).
method
dim
co-bps f
vel R2
PSTH R2 1
fp-bps f
Spike smoothing
137
0.2076
0.6111
-0.0005
GPFA
52
0.2463
0.6613
0.5574
SLDS
38
0.2117
0.7944
0.4709
-0.1513
LFADS
100
0.3554
0.8906
0,6002
0.2454
NDT
274
0.3597
0.8897
0.6172
0.2442
Ours
36
0.3225
0.8479
0.5927
0.2184
Ours (non-causal)
36
0.3407
0.8902
0.5963
0.2417
B.6 Stimulas-conditioninng in monkey reaching task
For the experiment with stimulus-conditione dynamics in the monkey reaching task, tested the performance of models with and without the conditioning inputs. We found that the conditioning inputs allow the networks perform better velocity decoding at lower dimensionalities.
Table 3: Performance benefits of conditioning for monkey reaching task.
conditioning
dim
vel R2 f
5
0.7897 + 0.0687
w/o conditioning
0.8944 + 0.0039
8
0.9085 + 0.0048
16
0.9196 + 0.0041
with
5
0.8589 + 0.0493
conditioning
6
0.9018 + 0.0114
26

---

Following the analysis in Fig. 7 we also further visualized the model's match to the spiking statistics, including mean and standard deviation (SD) of spiking rate, and mean, SD, and coefficient of variation (CV) of inter-spike intervals. We observed a good match to the mean and SD of the spiking rate across all conditions. Match to ISI statisties is also quite reasonable given the noise observed between estimates of the statistics from train and test.
B.7 Comparison to an approximate method for finding fixed points
# fixed paints found
17-
approximate combined analytic
5000 25000 50000
# inverses
Supplementary Figure 13: Repetition of the experiment of Fig. 8 but now with a rank-2 RNN with 128 units. Again, we show the number of fixed points found as a function of the number of matrix inverses computed, with errorbars denoting the minimum and maximum amount of fixed points found over 20 independent runs of the algorithm.
Recently an approximate method for finding fixed points in piecewise-lineat RNNs was proposed 25] The method proceeds by randomly selecting a linear region (a configuration of Do. see Supplement 1.2 and calculating the corresponding fixed-point. If it is indeed a 'true' fixed point of the RNN (it is consistent with the assumed Dn), we store If the fixed point was inconsistent with the assumed Dg. we iteratively initialize Dt according to the 'virtual' fixed point found and calculate the new fixed point corresponding to this Dg, until we either reach 'true' fixed point or reach a certain amount of iterations. Then, we reinitialize at a randomly selected new configuration of Dn and repeat the procedure.
Under some conditions, the approximate method can be shown converge in linear time (|M ||a]]| < 1) g. where will be faster than our exact method - however in general the convergence of the approximate method strongly depends the dynamics of the networks. In partieular, there are reasonable settings where the approximate method fails to find all fixed points, such as of a rank-2, 128 unit RNN with 17 fixed points (trained similarly the teacher RNN of Fig. Be: Fig. ST3). While an in-depth study of the approximate method is out of scope. we hypothesize that the failure to converge is because when initializing with randomly selected Dos out of (D + 1)^ possible configurations, the approximate method tends to converges to the same set of Dos.
Our method is completely independent of the dynamics of the system and has fixed cost, after which one is guaranteed that all fixed points are found. However, we do note that there can be scenarios where our exact method is still too costly. In this scenario, we propose
the approximate method, with one adjustment we first pre-compute the subset of E D (N) configurations that can contain fixed points, and then initialize the approximate method using randomly selected Dq from this subset. Empirically, this leads to better convergence in least some scenarios (Fig. Fig. $(3)
For the approximate method, we used code from https://ggthbb.ccm////sswtt////t 2023 which was released with the GNU General Public License.
C
Additional details of low-rank RNNs
C.1
Discretisation
Given
dz dt
= -z()) + N (T((M(()) + +Z2(())
27

---

Using the Euler-Mant Maruyyaaa method with timestep As:
A
N (M + VA.
2++1 = (1
7 2z +
Izce:
and with ((, I), define a = -- - N N. and Eg Trr distribution used our experiments. Note the slight 'overloading' of here, as the discrete time indice / of e.g., Ze corresponds to the continuous time z((1 1)
x. we obtain the transition
C.2 Conditional generation
Given input weights H E NNN,, and stimulus R,,, we define our model
dx x(t) + Jo(x(t)) Hs(t) + Ex.
dt
Using the same assumptions as before, x can be described by R N, variables
dz dt
z(t) N ((Mz(t) + Hz(()) + Ex.
ds
8(L) + s((),
di
with x = Mz + H5, and
(IM, HJ M,
HJ'x
We can write the distribution generated after discretization as:
T
pi21:r, N1:1, $1:11 11:11 p(zz)(((5)))) ((5=) II
p( St-1)p(Ze 5t-1, 2
II
plye Sc, 2c),
1=2
[=1
p(z: (;-1, =
Ex), p(2)) = NN Exn);
p(5; 1; 8(as; +
a a--aa)s 1). p(5)) = 5(0).
where the mean of the transition distribution is F(S,, z.) = az, + (To(M + Hā:).
For constant input s, š will converge to s, and we can ignore the additional N, variables, assuming
x(0) Mz(0)+ + Similarly if varies on a time scale slower than т, S is good approximation
21 Here, for experiments where the input is constant context signal (Fig. 7 we substitute s for š and consider the dimensional system described by z (which now has additional conditioning on s):
p(ZITT, Vier S1:T) p(21) II p(z, St-1, 2t
II ply: St, 2t), t=1
1=2
p(2 St-1, 1-1) = NNF((sF(ss 2 11); En), p(z)) = N N Eni):
where F(s z.) az, OMMz Hs,).
C.3
Linear transformations of the latent space and orthogonalisation
Given
X++1 = ax; + M o(x:) + 6x
2+11 = 0Z; NT ((M,.) + 6x
with (O, Ez), (0, ) We can do any linear transformation of the latent dynamics z: à Az, long A has rank R, without changing the neuron activity x. To see this, define M MA -1 N AÑ, and cg ~ (O, AE,AT), T), giving us:
X =1=ax, ax; + MN o (x:) + +
2 = az, + N (M2,) + +
which will leave x unchanged, while our latents are expressed in a new basis. We typically got more interpretable visualization of the latents by orthonormalising the columns of M. Thus we applied for all visualisations after training A U,M, with M = U, where U are the first R left singular vectors of J MN
28

---

D Details of empirical experiments
D.1 Training details
D.1.1 Initialisation
Our models are (unless noted otherwise) initialized as follows:
Ni ~ U1 vx vx
My U
v vi!!
2
Wa N0,
R
H U
Miny Nan
h, U
WR VN'V
x 0,
+ .9.
Eg f 0111
Σx I, Uz1 f 0, where W and b the output weights and biases respectively. For Gaussian observations initialise Σy + .01I.
For experiments with Poisson observations, we jointly optimized a (usually causal) CNN encoder as part the proposal distribution. The CNN was conditioned on observations and predicted the mean and log variance of normal distribution. It consisted of common initial layers consisting of ID convolutions, with a GeLU activation function, and separate output convolution for the predicted mean and (log) variance. The CNN was initialized to the Pytorch (72) defaults, except for the bias of the log variance output layer, to which we added log(.01) term, such that the output matches the initially predicted variance of the RNN. The exact number of layers and channels are reported in the sections for each experiment.
For the teacher-student setups, we used as non-linearity x(x;) = max(x; h, 0) for both the students and the teachers, and for all experiments with real-world data, we used the 'clipped" ((x:) max(x; - maxxxx;; 0) [
D.1.2 Parameterisation
We constrain a be between 0 and 1 by instead optimising à with the following (sigmoidal) parameterisation a = exp(- exp(-exp(ā)) ex)( (73.. In experiments with the optimal proposal, we estimate the full Σz, which we constrain be symmetric positive definite, by optimizing a lower triangular matrix C such that Σg = CC where we additionally constrain the diagonal of C to be positive using C; = exp( Ca:/2) 22). For all diagonal covariances, we parameterise the diagonal elements using Σg exp(E.)). For Poisson observations, we apply a Softplus function to rectify the predicted rate.
D.1.3 Optimisation
with stochastic gradient During training we minimise the variational SMC ELBO 33
(Eq 7
descent, using the RAdam (74 optimiser in Pytorch [72] We generally use an exponentially decaying learning rate (details under each experiment).
D.2
Teacher student experiments
D.2.1 Dataset description
We created datasets by first training 'teacher' RNNs to perform task and then generating observations by simulating the trained teacher RNNs.
29

---

For Ba b we used code from 22 http://tthuub.oommmccllllele licence) to train rank-2 2 RNNS to produce oscillations, using sine-wave with periodcityy of 50 time-steps as a target and an additional I.2 regularisation on the rates. After training, we extracted the recurrent weights M. N and biases h, orthonormalized the columns of M, and created a dataset by simulating the model for 75 timesteps, with Σ, .04I. For Fig. Ba we used N = 20 units and generated observations according to = (Mzz, ) with Σy = .011. Fig. 3 we used = 40 units and generated observations according to G Pois(Softpluss) Mu. b). with b 3.
Apache
4 and
For Fig. Be. we followed similar procedure but now trained the teacher RNN on a task where it has to use imput. After an initial period of 25 time steps, stimulus was presented for 25 timesteps consisting of sinn(f)coos()) where 9 was randomly selected every trial out of 8 fixed angles. The RNN was tasked to produce output that equals the transient stimulus for the next 100 time-steps. Here we used N = 60 units, Eg .00251 and generated observations according to G (M2,, Ey), with by = .025II. The training data for the student RNN was included for each trial the corresponding stimulus.
D.2.2 Training details
The 'student' RNNs had 20, 40, 60 units, respectively and rank R = 2, matching that of the teacher linear Gaussian according to G (MM,, Σ
RNNs. For Fig. 3a. c. The observation model and we used the optimal proposal distribution. For Fig. 3 we used G = Pois(Softplun WMz, b)). with W diagonal matrix (scaling the output of each unit individually). For Fig. 3b, we used causal CNN encoder as part of the proposal distribution. It consisted of 3 layers, with kernel sizes (21, 1), and channels (64, 64, 2). We used (causal) circular padding.
For all three experiments, we used k = 64 particles, batch-sizes of 10, and decreased the learning rate exponentially from 10-3 3 to 10 For Fig. 3a we trained for 1000 epochs of 200 trials, for Fig. 3) for 1500 epochs of 400 trials and for Fig. 33 for 200 epochs of 800 trials. We used a workstation with NVIDIA GeForce RTX 3090 GPU for these runs. One model took about 3 to hours to finish training.
D.2.3 Evaluation setup
For Fig. we generated long trajectories of 10000 time-steps of data for both the student and teacher RNNs. To facilitate visual comparisons between student and teacher dynamics, we also orthonormalized the columns of the students weights M after training, and for Fig. picked signs of the columns of M such that the student and teacher match (note that after orthormalizing, the columns of M are equal to the non-zero singular vectors of the full weight matrix J, which are only unique up to sign flip). As noted before. this leaves the output of the model unchanged. The autocorrelation in Fig. 3a was computed by convolving a sequence of lag- 120 steps of data with itself (with duration 2x lag), and normalising such that lag=0 corresponds to a correlation of 1. We repeated this for 80 sequences starting at different time-points of the whole trajectory.
D.3
EEG data
D.3.1
Dataset description
We used openly accessible electroencephalogram (EEG) data from 42
43
.0.00/ ODC-BY licence). The data was recorded from human subject sitting still with eyes open (session (001R01), and sampled at 160 Hz. Like 8. we used the full minute of recording, but unlike [8) we did not smooth the data (but just standardized the data). Thus, to compare our performance to (8), who ran their evaluation using the smoothed data, we smoothed our generated samples equivalently, using Hann filter with a window length of 15-time bins, so that we can also compare our samples to the smoothed data.
D.3.2 Training details
We used N = 512 units, and rank R = 3. The observation model was linear Gaussian conditioned on the hidden state and we used the optimal proposal distribution. We trained for 1000 epochs consisting of 50 batches of size of 10, and l 10 particles. The learning rate was decreased
30

---

exponentially from 10 10-3 to 10-6, Models were trained using NVIDIA RTX 2080 TI GPUs on compute cluster, single model took between 4 and 5 hours finish training.
D.3.3 Evaluation setup
We used our RNN to generate long trajectory of T = 9760 steps of data, Yt (after discarding the first 2440 steps), which we compare to the EEG data, y,, using two evaluation measures from 8 24] (using code from
GNU General Public License):
This is an estimate of the KL divergence between the ground truth and generated states. To compute this, obtained kernel density estimates of the probability density functions (over states, not time), using a Gaussian kernel with standard deviation 0 1. We get for the EEG data:
Dstsp:
ply) = + WSte I), and for the generated data g(y) = + E Wyr,, I). We then used the
following Monte Carlo estimate of the KL divergence: Datsp R 1 Σ log P(y') using n 1000 samples y' drawn randomly from the EEG data.
g(y")
DH: This is an estimate of the difference in power spectra between the ground truth and generated states. We first computed for each data dimension the spectra y yu for the EEG and generated data, respectively. We used a Fast Fourier Transform, smoothed the estimates with a Gaussian kernel with standard deviation 20, and normalized the spectra so they to 1. We computed the mean of the Hellinger distances between the spectra: Dи
64 Ei V
vyull-
D.4
Hippocampus HC-2
D.4.1 Dataset description
We used openly accessible neurophysiological data recorded from layer CAI of the right dorsal hippocampus [47 (8) (https:/ccress.orr http:/s..oo/aattaaaaaooaa wets/hc/he-2/abbuu/
Signals were recorded the rats engaged in an open field task, chasing drops of water or pieces of food that were randomly placed. We used the session ec013.527 from rat ID ec13, which is approximately 1062 seconds long. From 37 units (neurons) we used 21 neurons that have maximal spike counts, discarding the rest of the comparatively silent neurons. We binned the spike data 10ms. We used the first 80 percent of the data for training, and the rest was saved for testing purposes.
D.4.2 Training details
We used N = 512 units, and rank 3 for the run that was used our Fig. 5 CNN encoder as part of the proposal distribution, which consisted of 3 layers with kernel sizes (150, 11, 1), with (64, 64, 3) channels. During our study, we swept over multiple ranks and found that theta oscillations consistently emerged from rank 3 onwards, after which reconstruction accuracy was relatively stable. For each rank, we used three different seeds and two different first layer sizes for the encoder, 25 or 150. The duration of randomly sampled trial (sequence length) from the whole data was 94 time steps when the first layer size was 25, and 219 when the first layer size was 150. We, however, also found that the choice of the duration did not affect the results much. We trained the model using 3000 epochs, each epoch consisting of 3000 trials with 64 batches and k k= 64 particles. The learning rate decreased exponentially from 10 10-3 to 10 (0-6 A single model took approximately 21 hours to finish training on a NVIDIA RTX 2080 TI GPU on compute cluster.
We used causal
D.4.3 Evaluation setup
We used our RNN to generate data that matches the duration of the test data, which 20810 time steps (~2088 s) (after discarding the first 1000 steps). We compare different spike statistics of generated data with test data, and for comparison purposes, we also compared the same statistics measurements between train and test data as well. We calculated the mean firing rate of each neuron, mean of ISI distributions, and pairwise correlations. We used a band-pass filter 1-40 Hz for the latents and the LFP signal before calculating the powerspectrogram (Fig 5e).
31

---

D.5
Hippocampus HC-11
D.5.1 Dataset description
We used openly accessible neurophysiological data recorded from hippocampal CAI region 52 www.e//esss.g///aaattaeeee//eeee epoch, where rat was running
1/about-hc-11), We used the subset of the dataset called the maze 1.6-meter linear track, with rewards located at each end (left and right). Throughout this task, neural activity was recorded from 120 identified pyramidal neurons. As in [33, we only used 60 neurons that had sufficient actívity and discarded rest of the units. We used code from [53] btttps:/gitthub.cott com/ som//hd96/pi-valt
to preprocess the spike data, and only use data corresponding to the running and the location data being available for the results shown in Fig-
and Fig. STO For the model shown in Fig. STT we used the last 1350s of data of the Maze epoch,
which also includes bouts where the is stationary. We used 25ms bins.
D.5.2 Training details
We used N = 512 units, and rank R= 4 for Fig. and Fig. STO and rank R = 12 for Fig. STI We used a causal CNN encoder with zero padding with layers (24, 11, 1), (64, 64, 4) channels, and layers (150, 11, 1), (128, 64, 12) channels, respectively. The models were trained for 3000 epochs, each epoch having 3000 trials with a sequence length of 94 time bins (2.35 s) and 219 bins (5.48 s), respectively. We used batch sizes of 64 and k = 64 particles. The learning rate was decreased exponentially from 10-3 to 10-6, single model took approximately 21 hours finish training on NVIDIA RTX 2080 TI GPU on a compute cluster.
D.5.3 Evaluation setup
We used our RNN to generate data that matches the duration of the test data, which is 4289 time steps (~107 s) (after discarding the first 1000 steps) for Fig. and Fig. $[0 16539 time steps (~413 s) for Fig. STT We calculated the mean firing rate of each neuron, coefficient of variations of ISI distributions, and pairwise correlations. For the R reported in the Main text, we fit ridge regression model to the posterior latent variables on the training data, after smoothing with a Hann window of size 100, and apply the regression model to latents inferred from the test data.
D.6
Monkey Reach
D.6.1
Dataset description
We used the publicly available MC_ MC_Maze dataset from the Neural Latents Benchmark (NLB) [56] htw.aaadddddcciiieeoo///ddadddid
CC-BY-4.0 licence). The data were recorded from macaque performing a delayed center-out reaching task with barriers, resulting in variety of straight and curved reaches. For simplicity, we took only the trials with no barriers and thus straight reach trajectories, resulting in 592 training trials and 197 test trials. We binned the data at 20 ms and aligned each trial from 250 ms before to 450 ms after movement onset.
To create conditioning inputs for the model, we took the x and y coordinates of the target position for each trial and scaled them to be between -1 and 1. We then provide this scaled target position as constant context input to the RNN for the duration of the trial.
D.6.2 Training details
We ran random search of 30 different models with rank E 3,4, 6 and particle number & E 16, 32, 64. All models had 512 units and used causal CNN encoder with kernel sizes (14, 4, 2) and channels (128, 64, r). We used (causal) reflect padding. We trained each model for up to 2000 epochs, terminating training early no improvement was seen for 50 epochs. Each model took around 3 to 4 hours to train on an NVIDIA RTX 2080 TI GPU on compute cluster. Seeing that rank of was sufficient for velocity decoding R² as 0.9, we took the best-performing rank-5 model for subsequent analyses.
32

---

D.6.3 Evaluation setup
For qualitative evaluation of replication of cross-condition differences, grouped the reach targets in the data into 7 conditions, one at each corner and the midpoint of each edge of the rectangular reach plane, excluding the midpoint directly at the bottom. We then generated data from the model RNN using conditioning inputs from the test trials of the real data. Then, for the test data and the model-ggeneatee data, we computed mean firing rate and inter-spike interval for each neuron for each condition. We then computed correlation distance (1 r, where is the Pearson correlation coefficient) the neuron statistics between conditions in the test data and model-generated data.
For generation of data for Fig. 7d. e, we selected target locations by choosing angles from 0 360", evenly spaced by 22.5", and determined the corresponding reach endpoint
square spanning from (-1.-1) to (1, 1). We then constructed conditioning inputs similar to the real data using these target locations and simulated the RNN with them. To decode the reaches, we used linear decoder trained from inferred firing rates reach velocity from the real data.
D.7 Neural Latents Benchmark
D.7.1 Dataset description
We again used the publicly available MC MC_Maze dataset from NLB (see Supplement D.6.1). We resampled the data to 20 ms bin size and followed the standard data preprocessing procedures for the benchmark, as described in [56].
D.7.2 Training details
We ran a random search of 30 different models with varying rank from 12 to 40 and particle number I E 16, 32, 64, All models had 512 units and used a used a CNN encoder with kernel sizes (14, 4, and channels (128, 64, 36). either with causal reflect padding or acausal zero padding. We trained each model for up to 2000 epochs, terminating training early if no improvement was seen for 50 epochs. Each model took around 10 hours train on an NVIDIA RTX 2080 TI GPU on compute cluster.
Because the primary task of the benchmark is co-smoothing, i.e., prediction of held-out neuron firing rates from held-in neurons, we provide the encoder with only the activity of held-in neurons. However, the observation likelihood component of the ELBO is computed all neurons. held-in and held-out.
After training, we selected the model with the best co-smoothing score the validation split and submitted its predictions to the benchmark for the final evaluation.
D.7.3 Evaluation setup
Automated evaluation was performed on the benchmark platform, described in [56].
We used for the prediction timestep f, the expected Poisson rate of held-out neurons, conditioned on the activity of held-in neurons at the current and previous timesteps, by making use of the filtering Posterior (Eq. 3). We averaged over 32 sets trajectories with 192 particles each.
33

---

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: Our main claims a) That we can fit stochastic (ow-dimensional) low-rank RNNs to noisy neural data with variational SMC, and b) that we can analytically obtain the fixed points of the RNNs, if they sufficiently low-rank. Claim a) is supported by experiments with real datasets (after validation on student-teacher setups), where we obtain RNNs whose generated samples match the statistics of those datasets while having an underlying dynamical system that is 3-5 dimensional. Claim b) is supported by an extensive proof, as well as 2 numerical experiments.
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 in the paper and important assumptions and limitations. A No or
NA answer to this question will not perceived well by the reviewers.
The claims made should match theoretical and experimental results, and reflect how much the results be expected to generalize other settings.
It is fine include aspirational goals as motivation as 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]
Justification: In our discussion, we discussed restrictions on the noise model of our transition distribution, which currently has to generate (correlated) samples, such that they lie in the space of M. In addition, discuss the limitations of our analytic method for finding fixed points when the model's rank is not low.
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 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 (e.g., independence assumptions, noiseless settings, model well-specification, 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 the scope of the claims made, e.g., if the approach only tested on a few datasets or with a 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-text system might not be used reliably to provide closed captions for online lectures because it fails to 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 faimess.
While the authors might fear that complete honesty about limitations might be used by reviewers as grounds for rejection, worse outcome might be that reviewers discover limitations that aren't acknowledged in the paper. The authors should use their best
34

---

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 of assumptions and a complete (and correct) proof?
Answer: [Yes]
Justification: We have extensive proof of our theoretical result fixed points in low-rank piecewise-linent RNNs in the appendix.
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 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
Question: Does the paper fully disclose all the information needed to reproduce the main ex- perimental results of the paper to the extent that affects the main claims and/or conclusions of the paper (regardless whether the code and data are provided or not)?
Answer: [Yes]
Justification: Our Supplement contains for each experiment separate section that includes the exact experiment setup, hyperparameters, and evaluation details.
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 a specific model and empirical evaluation, may be necessary to 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 to accomplish this, but reproducibility instructions for how to replicate the results, access to a hosted model (e.g., in the case of large language model), releasing of model checkpoint, other means that appropriate to the research performed.
also be provided via detailed
While NeurlPS 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 it clear how reproduce that algorithm. (b) If the contribution is primarily the architecture clearly and fully.
model architecture, the paper should describe
(c) the contribution a new model (e.g., a large language model), then there should either be way to access this model for reproducing the results or way to reproduce the model (e.g., with an open-source dataset instructions for how to construct the dataset).
35

---

(d)
We recognize that reproducibility may be tricky in some cases, in which case authors are welcome describe the particular way they provide for reproducibility. In the case of closed-source models, it may be that access to the model is limited in some way (e.g., to 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 to faithfully reproduce the main experimental results, as described in supplemental material?
Answer: [Yes]
Justification: Code to reproduce results is available at sttts://giitun com/ mackelab/smc_ rnns
Guidelines:
The answer NA means that paper does not include experiments requiring code.
Please see the NeurlPS code and data submission guidelines (https://niis.cc)
Soniides/CodeSumiisst
for more details.
While we encourage the release of code and data, we 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 NeurIPS code and data submission guidelines (https
for more details.
The authors should provide instructions on data access and preparation, including how 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. If 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 much information possible in supplemental material (appended the paper) is 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: Hyperparameters, data splits and so on are described for each experiment in the Supplement.
Guidelines:
The answer NA means that the paper does not include experiments.
The experimental setting should be presented in the core of 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, or 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: Where-ever we show uncertainty, the figure or table caption reports what is over.
Guidelines:
36

---

The answer NA means that the paper does not include experiments.
The authors should answer "Yes" if the results are accompanied by enror bars, confi- dence intervals, or statistical significance tests, at least for the experiments that support the main claims of the paper.
The factors of variability that the error bars capturing should be clearly stated (for example, train/test split, initialization, random drawing of some parameter, or overall
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 the standard error of the mean.
is OK to report I-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 of range (e.g., negative error rates).
If error bars reported in tables plots, The authors should explain in the text how they were calculated and reference the corresponding figures or tables the text.
8. Experiments Compute Resources
For each experiment, does the paper provide sufficient information the com- puter resources (type of compute workers, memory, time of execution) needed to reproduce the experiments?
Question:
Answer: [Yes]
Justification: In the experiment details of the Supplement, we report the GPU used, as well the approximate compute-time.
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 compute required for each of the individual experimental runs 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 or failed experiments that didn't make it into the paper).
9. Code Of Ethics
Question: Does the research conducted in the paper conform, in every respect, with the NeurIPS Code of Ethics https: Ethhes/attpss///nurrppee
Answer: [Yes]
Justification: There is nothing that disagrees with any of the statements in the Ethics guideline.
Guidelines:
The answer NA means that the authors have not 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 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]
37

---

Our work is fundamental, we do not directly apply our approach to any problems with societal impact.
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 (c.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 as 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.
If there are negative societal impacts, the authors could also discuss possible mitigation strategies (c.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:
We do not foresee that our models have a high risk for misuse.
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 challenging, and many papers do not require this, but we 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: We included links to all assets used, cited the requested papers, and mentioned the licenses if available.
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 dataset.
38

---

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 use in the package should be provided. For popular datasets, paperswi 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 has changed) should be provided. this information 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 the documentation provided alongside the assets?
Answer: [Yes]
Our code includes new method to find fixed points in piccewise-lineal low-rank RNNs, available /ttttp:://gittub..ooo
Guidelines:
The answer NA means that the paper does not release new assets.
Researchers should communicate the details the dataset/code/moodl as part of their submissions via structured templates. This includes details about training, license, limitations, etc.
The paper should discuss whether and how consent was obtained from people whose asset 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 details about compensation (if any)?
Answer: [NA]
Justification: Our work does not involve crowdsourcing nor research with human subjects.
Guidelines:
The answer NA means that the paper does not involve crowdsoureing nor research with human subjects.
Including this information in 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 in the main paper.
According the NeurIPS Code of Ethies, workers involved in data collection, curation, or 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 equivalent approval/review based the requirements of your country or institution) were obtained?
Answer: [NA]
Justification: Our work does involve crowdsourcing nor research with human subjects.
Guidelines:
39

---

The answer NA means that the paper does not involve crowdsoureing nor research with
human subjects.
Depending on 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.
We recognize that the procedures for this may vary significantly between institutions and locations, and expect authors adhere the NeurIPS Code of Ethics and the guidelines for their institution.
For initial submissions, do not include any information that would break anonymity (if applicable), such as the institution conducting the review.
40


Reviews:
[{"review_id": "wfjC0JxKdt", "reviewer": "Reviewer_nuem", "summary": "This work focused on inferring stochastic low-rank structure from neural data. It developed a low-rank RNNs as state space models, and using sequential monte Carlo (SMC)  to learn the model's parameters. The proposed model is efficient in finding all fixed points in a polynomial cost instead of exponential cost. It is evaluated on multiple settings, one from student-teacher setups, given ground truth teacher RNN, and infer the structure and statistics, and demonstrated on extracting latents from EEG data, and recover interpretable latent dynamics from spike recordings in hippocampus, and allow position decoding.", "questions": "1. Explain the limited performance of model as mentioned above, how does GTF perform in the same task?\n\n2. Ablation studies of stochasticity and low-rank constraints.", "limitations": "No potential negative societal impact.", "rating": 5, "confidence": 4, "soundness": 3, "presentation": 3, "contribution": 2, "strengths": "**Motivation**\nThis work focused on important question in neuroscience, given partial observation and initial states, building stochastic models is essential to model complex neural dynamics. \n\n**Method**\nThe proposed method is efficient in terms of the low-rank structure, and as well as low complexity in finding fixed points. \n\n**Evaluation**\nThe method is demonstrated on extensive neural datasets, as well as synthetic networks. Evaluations are compared with deterministic model GTF.", "weaknesses": "**Novelty**\nLow-rank RNN has been demonstrated on many applications in recent works, the novelty is limited.\n\n**Experiments**\nAdding simulation with partial observations to demonstrate the effectiveness of stochastic low-rank RNN.\n\n**Results**\nLimited performance in recovering latent structure and predicting dynamics,  how does GTF perform in this task?\n1. Fig 3 a and b.\n2. Fig 4, the generated trace is not similar to EEG signal ground truth in the qualitative evaluation.\n3. Fig 5c, pairwise correlation is lower than 0.10.\n4. Fig 7e, the mismatch between true and generated in mean ISI dissimilarity.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We appreciate the reviewer's acknowledgement of the importance of the question we focus on, and the efficiency of our proposed methods for fitting low-rank models and finding their fixed points. We here hope to clarify some one of the raised concerns.\n\n>Low-rank RNN has been demonstrated on many applications in recent works, the novelty is limited.\n\nWhile there has indeed been extensive applications of deterministic RNNs, we here provide a method for fitting *stochastic non-linear* low-rank RNNs to single-trial neural data.  We also show the importance of using stochastic RNNs to capture observed neural variability.  Rebuttal Fig. B-D in the attached pdf highlight the importance of being able to infer the right level of latent noise. Related works (that we also cite) consider stochastic *linear* low-rank RNNs [1] or fit *deterministic* low-rank RNNs to trial-averaged data [2]. \n\n\n>Adding simulation with partial observations to demonstrate the effectiveness of stochastic low-rank RNNs\n\nThe low-rank setting naturally allows for partial observations. We have added a new experiment where we only observe 10 of the 20 units (Rebuttal Fig. A). The ground truth dynamics are still accurately captured.\n>Limited performance in recovering latent structure and predicting dynamics, how does GTF perform in the same task\n>Explain the limited performance of the model as mentioned above, how does GTF perform in the same task\n\nWe disagree with the reviewer's assessment of  \"limited performance\". In our student teacher setups we accurately capture the true latent dynamics (as well as the true level of latent noise, see Rebuttal Fig. G), and our empirical results on e.g., EEG data match state-of-the art. As noted in the summary response, we did not compare GTF on data with Poisson observations as GTF requires an invertible observation model, a limitation our method also overcomes. We have, however, now run GTF (deterministic) on one of the teacher-student setups with continuous observations, where (unlike our method) it  fails to capture the stochastic dynamics of the teacher model (Rebuttal Fig. A).\n\n>Fig 5c, pairwise correlation is lower than 0.10.\n\nFig 5c shows the pairwise correlation between neurons’ activity, plotted against the pairwise correlation between neurons of our fit network. The fact that spikes of real neurons have a low correlation with each other (the values of which our model captures) can hardly be seen as a limitation of our method. We will clarify this plot in the camera ready version.\n>Fig 4, the generated trace is not similar to EEG signal ground truth in the qualitative evaluation\n\nWhile we agree with the reviewer that there is no perfect match between generated and ground truth traces, we note that quantitatively the (smoothed) traces generated by our model match state-of-the art, while needing only 3 latent dimensions. The best-performing method in [3] used 16 latents.\n>Fig 7e, the mismatch between true and generated in mean ISI dissimilarity\n\nWe note that the ISI distribution is generally relatively noisy, see for instance the shift between mean train and test ISI per condition (supplementary figure 6), which can be in the 100s of ms for some neurons. We would argue that, given this, we do capture the differences in distributions between conditions well: The median absolute deviation of off-diagonal elements in the dissimilarity matrices is below .1 between our simulated data and the test set.\n>ablation studies of low-rank constraints\n\nWe have now fit full-rank RNNs to the EEG dataset, both models with 30 units (which have a similar number of parameters to our rank 3, 512 unit RNN), and models with 128 units (which have over 10x more parameters). We note that the KL divergence between data and simulated data is higher for the full rank models, while additionally being less tractable than our low-rank RNN (Rebuttal Fig. F).\n>ablation studies of stochasticity\n\nWhile our method is fundamentally probabilistic, GTF can be used to fit deterministic low-rank RNNs (we describe the relation in the manuscript section 2.2.2). We now show that when the true underlying model is stochastic, fitting a deterministic model with GTF does not allow one to accurately capture the true underlying model (Rebuttal Fig. B). Besides this, empirically, on the EEG dataset, stochasticity allowed us to obtain a reconstruction similar to GTF with only 3 latents, as opposed to having to learn a complex deterministic 16-dimensional chaotic system.\n\n[1] Valente, Pillow & Ostojic 2022. Probing the Relationship Between Latent Linear Dynamical Systems and Low-Rank Recurrent Neural Network Models\n\n[2] Valente, Pillow & Ostojic 2022. Extracting computational mechanisms from neural data using low-rank RNNs", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "hwLKuoZoJB", "reviewer": "Reviewer_drM9", "summary": "The authors describe an elegant method to infer a low-dimensional description of stochastic neural dynamics using variational Sequential Monte Carlo and Low Rank RNNs. They apply this method to simulated data and to three experimental datasets (EEG, hippocampus, and motor cortex). They also describe an elegant analytic approach to determine the fixed points of the learned dynamical system.", "questions": "n/a", "limitations": "The limitations are well discussed.", "rating": 8, "confidence": 5, "soundness": 4, "presentation": 4, "contribution": 4, "strengths": "The paper is very well-written. The motivating problem is well-established as an important question in computational neuroscience. The discussion and summary of prior work on this topic is executed very well. In addition to showing strong empirical results in simulation and real experimental analysis, the authors provide an illuminating mathematical analysis of low-rank RNNs and a procedure to efficient numerical procedure to find fixed points.\n\nA key modeling advance is the ability to fit a model with stochastic and nonlinear dynamics. This is in contrast to most RNNs, which currently use a deterministic transition.", "weaknesses": "I enjoyed this paper and found very few weaknesses. However, I would ideally like to see the variational SMC method benchmarked against simpler baselines -- mostly I would be interested in SMC methods with simpler proposals (e.g. bootstrap particle filtering). Additionally, I would like to see more details about performance as a function of (a) number of particles, and (b) amount of observed data.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for their positive feedback on our manuscript and their appreciation for modeling RNNs with stochastic dynamics! We have provided some additional analyses based on your suggestions, which we believe strengthens the original manuscript.\n\n>mostly I would be interested in SMC method with simpler proposals (e.g. bootstrap particle filtering).\n\nWe included new experiments using the bootstrap proposal. When rerunning the student-teacher setup of Fig 3A, the student networks are not reliably able to recover the true latent noise of the teacher if the bootstrap proposal is used instead of the optimal proposal, and generally has a large variance in performance (Rebuttal Fig. G). When fitting models to the HPC-2 data, both the quality of the generated data, as well as the Hellinger distance between the power spectra of the latents with that of the local field potential is worse when the bootstrap proposal is used (Rebuttal Fig. H).\n\n>I would like to see more details about performance of (a) number of particles and (b) amount of observed data\n\n(a) We included new experiments where we vary the number of particles. Concretely, for the teacher student setup of Fig 3A, we show that with 1 particle (i.e., no resampling) the true underlying noise of the teacher network is not recovered, unlike when using 64 particles (Rebuttal Fig. G). When using 10 particles we still get close (in line with previous work on variational sequential monte carlo, where it is suggested to use a small number of particles during training for efficiency, and then a larger number for evaluation). For the HPC-2 dataset, we similarly obtain more informative latents when using multiple particles (Rebuttal Fig. H). \n\n(b) We also included a new experiment where we ran both our method and GTF on 30 seconds (instead of 60) of the EEG data (Rebuttal Fig. I).  While both methods have a slightly worse KL-divergence score ($\\mathsf{D\\_{stsp}}$), this is in line with train/test shift between the first and second 30s of the EEG data. Additionally, we included a new student-teacher experiment (Rebuttal Fig. A) where the student only sees half of the teacher's neurons and still learns to approximate the teacher’s latent dynamics.", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "o3GpDWzpYC", "reviewer": "Reviewer_bWek", "summary": "This paper proposes a low-dimensional, nonlinear dynamics model of neural data based on low-rank RNNs. In particular, the model dynamics are a discretization of the dynamics of a low-rank RNN in the space spanned by the column factors of the low-rank RNN matrix. The model is fit using variational SMC. The approach is validated in simulation and applied to three neural datasets of EEG, hippocampal, and motor cortex recordings. Additionally, for low-rank RNNs the authors introduce an exact fixed point finding procedure that runs in polynomial time.", "questions": "- The authors point out in the discussion that the proposed model can be viewed as a neuralODE, when discussing the FINDR method. This paper could be improved by more clearly and directly discussing how the model and inference approach relate to other models and fitting methods in the neuralODE literature (another neuro example is ODIN from Versteeg et al). Since this model can be viewed as a neuralSDE, should it be using techniques from the neuralSDE literature for fitting? \n\n- Generally, low-rank RNNs can have transient dynamics outside the column space of the weight matrix factors, e.g. via inputs as the paper mentions. The current approach may be limited to handling only certain types of inputs. I think it would be helpful for the authors to describe this difference and the limitations more clearly in the main text.", "limitations": "Yes", "rating": 6, "confidence": 4, "soundness": 4, "presentation": 3, "contribution": 3, "strengths": "The modeling approach and fitting methods are generally clearly described. The model is applied in multiple synthetic and real applications, where it recovers known ground truth or produces sensible results. Additionally, the proposed analytical technique for finding fixed points of low-rank RNNs appears to be generally useful for anyone working with such models.", "weaknesses": "It appears that the formulation of the model in the low-dimensional space loses some generality relative to the original low-rank RNNs. As the authors point out, for non-constant inputs the activity can lie outside the low-d subspace. While the proposed approximation in appendix C.2 may work for the settings in this paper, it is not necessarily clear how the approach would perform if the inputs to the original low-rank RNN are significantly non stationary. \n\nAll of the models fit in this paper have small latent dimensionalities (typically 2-3 latent dimensions). It is not clear how well the approach scales to larger dimensionalities than those considered here.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for acknowledging that the paper is generally clear, and we hope to address their concerns, in particular with respect to when the inputs are time-varying.\n>It appears that the formulation of the model in the low-dimensional space loses some generality relative to the original low-rank RNNs. As the authors point out, for non-constant inputs the activity can lie outside the low-d subspace. While the proposed approximation in appendix C.2 may work for the settings in this paper, it is not necessarily clear how the approach would perform if the inputs to the original low-rank RNN are significantly non stationary\n\nWe would like to clarify that our method also works in settings where the inputs are non-stationary, although unfortunately we did not describe or use this in the main paper. The reviewer is correct in that with time-varying inputs, the activity of the neurons in the network can lie outside of the column space of the recurrent weights: it will however still be constrained to be in the the space spanned by the columns of the recurrent and the columns of input weights (as we also describe in  first paragraph of supplement C.2). While for the models used in the paper, we indeed make an approximation ($\\tilde{\\mathbf{s}} \\approx \\mathbf{s}$), that allowed us to ignore the additional input dimensions, one does not need to make this approximation, and can simply consider the augmented system of $[\\mathbf{z}, \\tilde{\\mathbf{s}}]$.\n\nWe will explicitly include the equations for the distribution generated by the augmented system in the camera-ready version of the manuscript to clarify this.\n\nTo explicitly demonstrate this, we now ran new experiments with a student-teacher setup, where a rank 1 teacher network was trained to report the sign of the mean of a time-varying noisy input (Rebuttal Fig. E).  We visualize the closely matching dynamics of the student and teacher network dynamics, in the plane spanned by the one column of the left recurrent connectivity vector $\\mathbf{M}$ and the input weights vector (coordinates $\\mathbf{z}$ and $\\tilde{\\mathbf{s}}$ respectively). This demonstrates that our method naturally works, even when there are time-varying inputs.\n\n>All of the models fit in this paper have small latent dimensionalities (typically 2-3 latent dimensions). It is not clear how well the approach scales to larger dimensionalities than those considered here.\n\nWe would like to point out that our appendix already contained a model with 36 latents (Supplementary table 2). Additionally Rebuttal Fig. F contains full-rank models (rank 30 and rank 128) trained on the EEG dataset. We have since also trained new models on the whole rat HPC 11 dataset (instead of only the subset of data where the rat is moving), and we get good performance with 12-16 latents, we would be happy to include these in the camera-ready version of the paper. \n\n>The authors point out in the discussion that the proposed model can be viewed as a neuralODE, when discussing the FINDR method. This paper could be improved by more clearly and directly discussing how the model and inference approach relate to other models and fitting methods in the neuralODE literature (another neuro example is ODIN from Versteeg et al). Since this model can be viewed as a neuralSDE, should it be using techniques from the neuralSDE literature for fitting?\n\nWe thank the reviewer for this comment, and will include a more extended discussion of related works as ODIN and neuralSDEs in our camera-ready version. A prominent line of work in the neural ODE/SDE literature concerns overcoming memory requirements of backpropagating through operations in an ODE solver, by for example using the adjoint method (i.e., substituting 'optimize then discretise' for 'discretise then optimize'). We here — similar to FINDR and ODIN — do not use the adjoint method, but rather a simple Euler-Maruyama discretisation scheme and 'standard' backpropagation through time. We note that this generally performed well for our use-cases. However, one could investigate how we can integrate our approach with variational approaches that use adjoint methods when fitting latent neural SDEs [1,2] as well as with filtering approaches for continuous time systems [3]. This could be especially relevant for irregularly sampled time-series.\n>Generally, low-rank RNNs can have transient dynamics outside the column space of the weight matrix factors, e.g. via inputs as the paper mentions. The current approach may be limited to handling only certain types of inputs. I think it would be helpful for the authors to describe this difference and the limitations more clearly in the main text.\n\nPlease see our response to the first comment and Rebuttal Fig. E — our method naturally works with any kinds of inputs (irrespective of whether they are time-varying or not).\n\nReferences:\n\n[1]  Li, Wong, Chen & Duvenaud. 2020. Scalable Gradients for Stochastic Differential Equations\n\n[2] Deng, Brubaker, Mor & Lehrmann. 2021. Continuous Latent Process Flows  \n\n[3] Särkkä & Sottinen. 2008. Application of Girsanov Theorem to Particle Filtering of Discretely Observed Continuous-Time Non-Linear Systems", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "Myu5Aklz5g", "reviewer": "Reviewer_Strk", "summary": "The authors propose fitting stochastic low-rank RNNs to neural recordings using variational sequential Monte Carlo methods. Such techniques permit modeling of noisy sequences (i.e., trial-to-trial variability), identification of a low-dimensional latent dynamical system, generative sampling of neural trajectories with realistic variability, and interpretation via fixed point analysis. The technique is applied successfully to recover the ground truth dynamics in two synthetic systems, and then to model EEG recordings, hippocampal spiking data, and motor cortical spiking data.", "questions": "- What is meant exactly by \"tractable dynamics\"? Is there a distinction between \"tractable\" and \"interpretable\"?\n- Why is the proposed technique compared to Generalised Teacher Forcing in the EEG experiments, but not in the synthetic setups (Fig 3) or the spiking data (Figs 5-7)?", "limitations": "- The authors note just one assumption / limitation of their approach: an assumption of correlated Gaussian noise in the recurrent dynamics. How does the technique fair when observations from the true underlying system reflect private noise processes (e.g., measurement noise, variability in neurotransmitter release, etc)?\n- One limitation that was not discussed or addressed is that the approach does not model the effects of unobserved inputs (beyond noting that correlated noise in the dynamics may arise due to unobserved inputs). Does this imply that the technique as presented is only appropriate for settings where the true dynamics can be reasonably modeled as an autonomous system?", "rating": 6, "confidence": 2, "soundness": 4, "presentation": 3, "contribution": 3, "strengths": "- Originality: The technique presented appears original in its combination of existing ideas from low-rank RNNs, variational inference, and sequential Monte Carlo. \n- Quality: Thoughtful comparisons were made to existing approaches. The results improve upon state-of-the-art techniques (e.g., Generalised Teacher Forcing) in certain settings.\n- Clarity: The figures and tables are clearly presented and quite interpretable. Much of the writing is clear, although see Weaknesses and Questions for suggestions here.", "weaknesses": "- Section 2.2 could benefit from being made more accessible to readers who are not experts in the subdomains of variational inference and sequential Monte Carlo methods. \n- The authors should make clear how to explicitly implement the technique.\n- The approach does not outperform state-of-the-art techniques in the Neural Latents Benchmark (NLB). While the authors mention that the \"NLB metrics center around evaluating the quality of smooth rates inferred from spikes, which is not the central focus of our method. Rather, we aim to fit an RNN, from which -- by design -- we can sample noisy latent trajectories that reproduce variability in the data.\" But doesn't LFADS (and NDT?) allow generation of noisy latent trajectories (in the LFADS \"factors\") that reproduce variability in the data?", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for recognizing that our approach \"appears original in its combination of existing ideas” and that we provided “thoughtful” evaluations and demonstrations of our method. We clarified some raised issues below\n>Section 2.2 could benefit from being made more accessible to readers who are not experts in the subdomains of variational inference and sequential Monte Carlo methods.\n>The authors should make clear how to explicitly implement the technique.\n\nWe would be happy to make these points more clear in the camera-ready version! We did provide additional implementation details in the appendix plus some of the code (which we by now also cleaned and re-organised, and will link in the camera-ready version).\n\n>But doesn't LFADS (and NDT?) allow generation of noisy latent trajectories (in the LFADS \"factors\") that reproduce variability in the data?\n\nThe reviewer is correct that LFADS can also be used to generate latent trajectories (instead of inferring them), though it is not commonly used in this way. We now explored sampling from LFADS, both using the autonomous version, and when using the controller, after training on the rat HPC 2 dataset (Rebuttal Fig. C, D). In the case of an autonomous LFADS model, one samples an initial condition from the prior and then simulates a deterministic RNN forward. As on long sequences not all variability can be explained by variability in the initial condition, the latents end up not representing any variability that resembles the underlying system (Rebuttal Fig. C, cf. Manuscript Fig. 3). In the case of a full LFADS model with the controller, one can sample both an initial condition and time-varying inputs from the controller’s autoregressive prior. The full model seemed to rely overly on the controller’s data-inferred inputs, which deviated quite strongly from the samples from the controller’s autoregressive prior (Rebuttal Fig. D left). As a consequence, the latents do not seem to represent variability that is meaningful (Rebuttal Fig. D right).\n\nIn general, accurately inferring firing rates from data does not necessarily translate to generating realistic data. NLB metrics do not evaluate generative quality, and in fact many high-performing methods in NLB are not used as generative models at all (e.g., NDT). We stress that performance on NLB is not the central focus of our method, and when it comes to generation our method outperformed relevant baselines.\n\n>What is meant exactly by \"tractable dynamics\"? Is there a distinction between \"tractable\" and \"interpretable\"?”\n\nWith “tractable” we mean that properties of the dynamics, such as fixed points, can be obtained (exactly) in a computationally feasible way. The ability to compute fixed points of the dynamics efficiently can aid in interpretation and understanding the model’s underlying computations (see e.g., [1]). We will clarify this in the camera-ready version.\n\n>Why is the proposed technique compared to Generalised Teacher Forcing in the EEG experiments, but not in the synthetic setups (Fig 3) or the spiking data (Figs 5-7)?\n\nWe thank the reviewer for pointing this out. We have now run GTF on one of the teacher student setups of Fig. 3A (see Rebuttal Fig. B). The deterministic student model that is fit with GTF doesn’t capture the latent dynamics of the stochastic teacher model well, and generated samples from the deterministic student are not as similar to the teacher model as when the student is fit with our method.\nAs noted in the summary response, we did not compare GTF on data with Poisson observations, as GTF requires an invertible observation model, a limitation our method also overcomes. \n\n>How does the technique fair when observations from the true underlying system reflect private noise processes (e.g., measurement noise, variability in neurotransmitter release, etc)?\n\nMeasurement noise that is unique to every neuron can be (and is in our manuscript) straightforwardly taken into account by choosing the appropriate observation model. Private noise that is iid over neurons and feeds back into the recurrency, reflecting e.g., variability in the neurotransmitter release, is however currently not modeled. There are cases where one can show that iid noise has little influence on the underlying population dynamics – e.g., in linear models with many more units than latent dimensions, where such noise becomes effectively orthogonal to the recurrent weights. However, as noted in our discussion, an interesting future direction would be investigating how to best incorporate private noise processes that feed into the recurrence. \n\n>Does this imply that the technique as presented is only appropriate for settings where the true dynamics can be reasonably modeled as an autonomous system?\n\nThe reviewer raises a valid point. Our method does not explicitly model unobserved inputs. However, in principle, the model’s recurrent dynamics can capture unobserved inputs, as long as these inputs can be modeled as being generated by Gaussian noise, or by a subnetwork of the RNN (or in other words, as terms of a stochastic differential equation). Previous work on low-rank RNNs has focused on finding subpopulations in trained low-rank RNNs [2], which could potentially be applied here to find out which subpopulations in the trained RNN provide input to each other. Such an approach has been successfully applied to (deterministic) RNNs fit to trial-averaged data in [3]. \n\n[1] Sussillo & Barak, 2013, Opening the black box: low-dimensional dynamics in high-dimensional recurrent neural networks\n\n[2] Dubreuil, Valente, Beiran, Mastrogiuseppe & Ostojic. 2022. The role of population structure in computations through neural dynamics\n\n[3] Valente, Pillow & Ostojic 2022. Extracting computational mechanisms from neural data using low-rank RNNs", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}]

Author rebuttal:
**General response**

We thank the reviewers for their extensive comments and insightful feedback on our manuscript. Our paper introduced a method for obtaining low-dimensional descriptions of stochastic neural dynamics, tackling a “well-established (drM9)” and “important question in neuroscience” (neum). Additionally, we proposed a method for finding fixed points in piecewise-linear low-rank RNNs, using an “elegant analytic approach” (drM9) that “appears to be generally useful for anyone working with such models” (bWek).

The reviewers  appreciated the "strong empirical results" (drM9) that “improve upon state-of-the-art techniques (e.g., GTF) in certain settings" (Strk), which were ”demonstrated on extensive neural datasets, as well as synthetic networks" (neum). Moreover the reviewers provided positive feedback on the presentation, stating that we provided an "illuminating mathematical analysis of low-rank RNNs" (drM9) and that the "figures and tables are clearly presented and quite interpretable" (Strk).

Nevertheless, the reviewers highlighted important points which we have now addressed with new figures and explanations. We here summarize the main new analyses, which fully support our original claims.


**Importance of fitting stochastic models**

Both reviewer neum and Strk ask for more comparisons to generalised teacher forcing  (GTF) [1] — a related method for fitting RNNs with deterministic transitions, as we discuss in the manuscript. We repeated one of the student-teacher setups by fitting a student with GTF, which demonstrates that — if the underlying RNN is indeed stochastic — deterministic GTF is not adequate (Rebuttal Fig. 1B). Thus, one needs methods for stochastic dynamics  (such as the one we propose), both to obtain good reconstructions and to recover the true latent dynamics. This point is also reinforced by a new experiment where we fit LFADS [2] to the HPC2 dataset (Strk; Rebuttal Fig. 1C,D).
Furthermore, in the formulation of GTF in [1], the authors require an invertible observation model (see section 3.4 in [1]). We can therefore not directly compare our method to GTF in the experiments with spiking observations (Manuscript Figs 5-7). Indeed, one advantage of our method is that it overcomes the limitation of needing an invertible observation model (by using an encoding network that predicts a distribution over latents).


**Additional analysis of hyper-parameters and training setup**

We performed additional experiments where we demonstrate that: 
- In a teacher-student setup, we can recover the underlying latent dynamical system, even when we only fit to  partial observations (neum; Rebuttal Fig. 1A).
- We can recover the true latent noise in student teacher setups, when we use the optimal proposal ($ p(\mathbf{z}\_t | \mathbf{y}\_t, \mathbf{z}\_{t-1}) $ ) with 64 particles, and can get close when using 10 particles. When using a bootstrap proposal ($ p(\mathbf{z}\_t|\mathbf{z}\_{t-1})$) or only 1 particle (i.e., no resampling) this is no longer the case (drM9; rebuttal Fig 1G). The bootstrap proposal also leads to worse performance on the HPC2 dataset (drM9; rebuttal Fig 1H). This indicates that our strategy of using multiple particles and conditioning the proposal distribution on observed data, is indeed beneficial relative to alternatives.
- We fit full-rank models (32 and 128 units) to the EEG dataset, and show that the performance of full rank RNNs is worse than those of  low-rank RNNs, on top of the models being less tractable (neum; Rebuttal Fig. 1F).


**Generalization to time-varying inputs**

Reviewer bWek asks how the approach would perform if the inputs to the original low-rank RNN are non stationary. We unfortunately did not describe or use this in the main paper, but our method works with time-varying inputs (see the equations in the first paragraph of Supplement C.2 which define a low-rank RNN with arbitrary inputs). Our probabilistic setup therefore also allows fitting models with time-varying inputs. To demonstrate this, we now added a new teacher-student setup where the teacher network integrated a time-varying stimulus (Rebuttal Fig. 1E), and show that the student network successfully learned the latent dynamics of the teacher using our method. We will update the paper accordingly.

[1] Hess, Monfared, Brenner & Durstewitz. 2023. Generalized Teacher Forcing for Learning Chaotic Dynamics  
[2] Pandarinath, O’Shea, Collins, Jozefowicz, Stavisky, Kao, Trautmann, Kaufman, Ryu, Hochberg, Henderson, Shenoy, Abbott & Sussillo. 2018. Inferring single-trial neural population dynamics using sequential auto-encoders

Meta-review: This work introduces stochastic low-rank RNNs, trained by sequential Monte-Carlo in a way that is related to generalized teacher forcing. The authors further describe a nice technique for finding all fixed points of the low-rank model in polynomial rather than exponential time. In addition to explicitly accounting for noise in the process, the proposed method finds much lower-dimensional latent spaces which is a major plus. While some issues with novelty and lack of comparisons were brought up in the reviews, these were largely addressed in the rebuttal, and all referees finally supported acceptance. I fully agree and think this paper is not only technically strong and solid, but contributes substantial advances to the field. 

Just a little note from my side: The authors claim (and show for a very particular example) that GTF cannot fit stochastic dynamics, but in the cited paper it is indeed trained on systems with noise and seems to perform quite well?
Decision: Accept (poster)