SIMPL: Scalable and hassle-free optimisation of neural representations from behaviour
Authors: Tom George
Pierre Glaser
Kim Stachenfeld
Caswell Barry
Claudia Clopath
Venue: ICLR 2025 Poster (2025)
Keywords: neuroscience; place cells; grid cells; representations; neural data; hippocampus;
Primary area: applications to neuroscience & cognitive science
TL;DR: A fast and simple technique for optimising neural representation

Abstract: Neural activity in the brain is known to encode low-dimensional, time-evolving, behaviour-related variables. A long-standing goal of neural data analysis has been to identify these variables and their mapping to neural activity. A productive and canonical approach has been to simply visualise neural "tuning curves" as a function of behaviour. However, significant discrepancies between behaviour and the true latent variables -- such as an agent thinking of position Y whilst located at position X -- distort and blur the tuning curves, decreasing their interpretability. To address this, latent variable models propose to learn the latent variable from data; these are typically expensive, hard to tune, or scale poorly, complicating their adoption. Here we propose SIMPL (Scalable Iterative Maximization of Population-coded Latents), an EM-style algorithm which iteratively optimises latent variables and tuning curves. SIMPL is fast, scalable and exploits behaviour as an initial condition to further improve convergence and identifiability. It can accurately recover latent variables in spatial and non-spatial tasks. When applied to a large hippocampal dataset SIMPL converges on smaller, more numerous, and more uniformly sized place fields than those based on behaviour, suggesting the brain may encode space with greater resolution than previously thought.

Published as a conference paper at ICLR 2025
SIMPL: SCALABLE AND HASSLE-FREB OPTIMISATION OF NEURAL REPRESENTATIONS FROM BEHAVIOUR
Tom M George
Pierre Glaser
Sainsbury Wellcome Centre, UCL tom. gomrge 20@ucl .c..uk
Gatsby Computational Neuroscience Unit, UCL
Kimberly Stachenfeld
Caswell Barry
Google DeepMind & Columbia University
Dept. of Cell and Developmental Biology, UCL
Claudia Clopath
Sainsbury Wellcome Centre, UCL Imperial College London c.clopath@imerrial..............
ABSTRACT
Neural activity the brain known to encode low-dimensional, time-evolving, behaviour-related variables. A long-standing goal of neural data analysis has been to identify these variables and their mapping to neural activity. A productive and canonical approach has been to simply visualise neural "tuming curves" func- tion of behaviour. In reality, significant discrepancies between behaviour and the true latent variables-such position X-distort and blur the tuning curves, decreasing their interpretability. To address this, latent variable models propose to learn the latent variable from data; these are typically expensive, hard to tune, scale poorly, complicating their adoption. Here we propose SIMPL (Scalable Iterative Maximization Population-coded Latents), an EM-style algorithm which iteratively optimises la- tent variables and tuning curves. SIMPL. fast, scalable and exploits behaviour as an initial condition to further improve convergence and identifiability. It can accurately recover latent variables in spatial and non-spatial tasks. When applied to a large hippocampal dataset SIMPL converges smaller, more numerous, and more uniformly sized place fields than those based on behaviour, suggesting the brain may encode space with greater resolution than previously thought.
agent thinking of position Y whilst located
INTRODUCTION
Large neural populations in the brain are known to encode low-dimensional, time-evolving latent variables which are, oftentimes, closely related to behaviour (Afshar et al. 2011 Harvey al. 2012
Mante et al.) 2013) Carnevale et al.] (015). Coupled with the advent of modern neural recording techniques (Jun al.] 2017 Wilt al. (2099) focus has shifted from single-cell studies the joint
analysis of hundreds of neurons across long time windows, where the goal is to extract latents using
variety of statistical (Yu al. 2008a| Cunningham
Yu 2014: Kobak et al. 2016 Zhao & Park
2017 Williams et al. 2020 Bjerke et al. 2022)) and computational (Van der Maaten & Hinton 2008
Pandarinath al. 2018 Mackevicius et 2019) methods.
This paradigm shift is particularly pertinent in mammalian spatial memory and motor systems where celebrated discoveries have identified cells whose neural activity depends on behavioural variables such as position Keefeefe & Dostrovsky 1971] Hafting et al.) (005), heading direction (Taube ct al. 1990), speed (McNaughton et al.) 1983), distance to environmental boundaries/objects (Lever al. 2009 Høydal et al. 2019) and limb movement direction (eeeriooooorroooooooo et al. 1986) through complex and non- non-linear tuning curves. Characterising neural activity in terms behaviour remains cornerstone practice in these fields however the implicit assumption supporting
that the latent variable encoded by neural activity is and only is the behavioural variable increasingly being
called into question (Sanders al. 2015 Whittington et al. 2020 George al. 2024b).

---

Published as a conference paper at ICLR 2025
The brain is not a passive observer of the world. The same neurons which encode animal's 's current position//ehavionn state are also used to plan a future routes (Spiers & Maguire, 2006), predict up-
coming states (Muller & Kubie 1989; Mehta et al.) 1997| Stachenfeld et al. 2017)
recalllrrpplay"
past positions (Squire et al. 2010) Carr et al. 201 necesarily causing the encoded latent variables to deviate from behaviour. Nor is the brain a perfect observer; uncertainty due to limited, noisy or ambiguous sensory data can lead to similar discrepancies. Measurement inaccuracies can contribute further. These hypotheses are supported by analyses which show that it is rarely, if ever, possible to perfectly decode "behhaviur" from neural data (Glaser et al., (020) and the observation that neurons show high variability under identical behavioural conditions (Fenton & Muller 1998 Low et all] 2018). All combined, these facts hint at a richer and more complex internal neural code. When this is not accounted for tuning curves will be blurred, distorted or mischaracterised relative to their true form. For example, consider an animal situated at position X 'imagining' or 'anticipating' a remote position, Y, for which a place cell is tuned. This might trigger the cell to fire leading to the mistaken conclusion that the cell has place field at location X.
Nonetheless, the fact that behaviour is often a coose-butiimpefee proxy for the true latent motivates searching for techniques which exploit this link. Most existing methods for latent discovery don't 't ex- ploit behaviour (Gao et al.) 2016, Gondur et al.] (2023) at the cost of complexity and interpretability.
Others don't 't model temporal dynamics(Zhou & Wei 2220; Schneider al. 2023 Lawrence] 2003).
don't scale to large datasets (Wang et al. [005]] Nam |013 [u et al.] 201 7, can't 't model complex
non-linear tuning curves (Pandarinath al. 2018| Hurwitz et al. 202 Duncker et al. 2019 Lin-
derman et al. 2016 Gondur et al. 2023). or aren't t designed for spiking datasets(Lawrence, 2003
Krishnan al., 2015). Moreover, many of these methods are conceptually complex, lack usable
code implementations, or necessitate GPUs limiting their accessibility.
Contributions Here we introduce SIMPL (Scalable Iterative Maximisation of Ropulation-coder Latents), a straightforward yet effective enhancement to the current paradigm. Our approach fits tuning curves to observed behaviour and refines these by iterating two-step process. First the latent trajectory is decoded from the current tuning curves then, the tuning curves are refitted based on this decoded latent trajectory. SIMPL imposes minimal constraints on the tuning curve structure, scales well to large datasets without relying on neural networks which can be expensive to train. Theoretical analysis establishes formal connections to aspectation-maximistt
(EM, Dempster et al. 1977) for a flexible class of generative models. By exploiting behaviour as an initialisation, SIMPL converges fast and helps mitigate local minima and identifiability (Hyvärinen & Pajunen 1999 Locatello al. 2019) issues. This allows it to reliably return refined tuning curves and latents which remain close to, but improve upon, their behavioural analogues readily admitting direct comparison. All in all, SIMPL is able to identify temporally smooth latents and complex tuning curves related to behaviour, while remaining cheap and natively supporting spiking data distinguishing set of features in the field of latent variable models for neural data analysis.
We first validate SIMPL on a dataset of synthetically generated 2D grid cells. Next, we apply SIMPL to rodent electrophysiological hippocampal data (Tanni et al., (022) and show it modi- fies the latent space in an incremental but significant way: optimised tuning curves are better at explaining held-out neural data and contain sharper, more numerous place fields allowing for a rein- terpretation of previous experimental results. Finally, we apply SIMPL to somatosensory dataset for a monkey performing a centre-out reaching task (Chowdhury et al.) 2020). SIMPL, with a 4D latent space, provides good account of the data with the latent variables initialised to (and remain- ing correlated with) the monkeys hand-position and hand-velocity. With only two hyperparameters, SIMPL can be run quickly on large neural datasets without requiring GPU. It outperforms pop- ular alternative techniques based on neural networks (Schneider et al. 2023 Zhao & Park 2017) or Gaussian processes(Lawreet
2003 Wang et al., 2005) and is over 15x faster. This makes it practical alternative to existing tools particularly of interest to navigational or motor-control commu- nities where abundant data is explained well by measurable behaviours (position, hand dynamics). We provide an open-source JAX- IAX-optimised (Bradbury et al. (2018) implementation of our codd")
'One-hour recordings of 200 neurons (10" spikes) takes 1 minute to run on a CPU laptop.
"Code and a demo be found at: httpa////ggiuuubb.oo////eeeeeeee 1234 1234/8 IMP1,
2

---

Published as a conference paper at ICLR 2025
2 METHOD
Here we give a high-level description of the SIMPL algorithm. Comprehensive details and a theo- retical analysis linking SIMPL to expecation-maxiiisttio,, are provided in the Appendix.
2.1 THE MODEL
SIMPL models spike trains of the form $ == (st) (1...= where Sti represents the number of spikes emitted by neuron i between time (t (6-1) 1) di and t dt. We denote == ( ,1N) the vector of spike counts emitted by all neurons in the t-th time bin. SIMPL posits that such spike trains modulated by latent, continmously-valueed low- now-dimesional,, time-evolving variable x = (x..,=1,....T E RP through the following random process:
i=.....N
X X₁
N(x,, o1
(Latent dynamics)
(1)
Sti
X₁
Poisson
(Emission model)
dt and Xo (00 (I). This generative model enforces tunable (through the velocity hyperparameter u) amount of temporal smoothness in the trajectories. At each time step the latent variable x, determines the instantaneous firing rate of all neurons via their intensity functions f. (hereon called tuning curves, collectively denoted f), which are unknown priori, and which SIMPL will estimate. Moreover, we make the common assumption that all neurons are conditionally independent given X₁, i.e. p(s,|x, II. ((St)(x)) Finally, we assume the latent variable x is Markovian, a common assumption in the neuroscience literature. This model has been previously studied in the literature (Smith & Brown (2003) Macke et al. 201T), albeit using highly restrictive tuning curve models, something which SIMPI. avoids.
where or
(a) SIMPL
M-step spike smoothing
"M-step" fit tuning curves
tuning curves, f(()
latent trajectory, x(t)
spikes,
s(t) true dynamics
E-step MLE + Kalman
1+1
A
8,
8
"E-step" decode position
argmax, log p(s x)
(c) behaviour
(d)
true dynamics
2
din
x(O)
f(0)
X
dim
Figure 1: Schematic of SIMPI.. (a) A latent variable model (LVM) for spiking data (x)()), x(t)) is optimised by iterating a two-step procedure closely related the pnnectation-maxiniisi to an initial estimate of the latent trajectory (an "M-step"). The latent is then redecoded from these tuning cuves (an "E-step"). (b) SIMPL fits tuning curves using kernel-smoothed estimate (lop) and decodes the latent maximum likelihood estimates. (c) Measured behaviour is used to initialise
First, tuning curves are fitted
variables by the algorithm as is often closely related the true generative LVM (d).
2.2 THE SIMPL ALGORITHM
Outline We now seek an estimate of the true, unknown latent trajectory x* and tuning curves f* that led observed spike train, SIMPL does by iterating two-step procedure closely related to
3

---

Published as a conference paper at ICLR 2025
the expectation expectation-maxiiisiii maximisation (EM) algorithm: first, tuning curves fitted to initial estimate of the latent variable (the "M-step"), which are then used to decode the latent variable (the "E-step"). This procedure is then repeated using the new latent trajectory, and so until convergence.
The M-step In the M-step (or "fitting" step) of the e-th iteration SIMPL fits tnng curves to the current latent trajectory estimate x(*) using a smooth, kernel-based estimate
21 Sti (xx, X; (e))
# spikes at x # visits to
(x)
Et=1 k (e)
for some kernel In practice, we use a Gaussian kernel with small bandwidth o. Such tuning curve model is conceptually simple and free from the optimisation, misspecification or interpretabil- ity issues of most parametric models. It constitutes notable departure from alternatives which use neural network (Zhou & Wei, [2220) [Schneider et al. (203)) to model tuning curves and is particularly well suited to low-dimensional latent spaces.
The E-step In the E-step SIMPL seeks to infer (or "decode") a new estimate of the latent from the spikes and current tuning curves, x ((++1) = ((<)) [x]. Directly performing this inference from the spikes is difficult due to the non-linearity and non- non-Gaussianity of the emission model in Eq (2) Instead, SIMPL. first calculates the maximum likelihood estimate (MLE) of x, denoted X. approximation p(x, x,) = Axx,; Σ.), the variables (x, X) form a Linear Gaussian State Space Model (LGSSM) fully characterised by of (the transition noise covariance) and Σ, (the observation noise covariance). This enables efficient inference via Kalman smoothing of the MLEs in order to approximate x(++1) [p(x{s][x] (schematic in Fig. [b).
Then, by making
( = arg max log p(|x( p(e)
(4)
x (e+1) = (p( [x] R KalmanSmnooth(((()) o0 E.)
Crucially, the non-linear and Gaussinn by design, but on p(X[x), a quantity which provably asymptotically Gaussian in the many-neurons regime (theoretical argument and an explicit formula for Σ, B.I..
approximation is not made on the spiking emissions p(s)x), which is
Behavioural initialisation Spike trains often come alongside behavioural recordings x" thought to relate closely to the true latent variable x² ~ x*. SIMPL. leverages this by setting the initial decoded latent trajectory, to measured behaviour x (0) + x². We posit that behavioural initialisation will place the first iterate of SIMPL within the vicinity of the true trajectory and tuning curves, accelerating convergence and favouring the true latent and tuning curves (x",f*) over alternative isomorphic pairs (o(x*), I o Φー² whose latent space is warped by an invertible map & but which would explain the data equally well. This amounts to an inductive bias favouring tuning curves close to those calculated from behaviour. Through ablation studies we confirm these beneficial effects.
All in all, SIMPL is interpretable and closely matches common practice in neuroscience (e.g. kernel-based curve fitting. MLE-based decoding); moreover, can be formally related to a gen- eralised version of the EM-algorithm, for which theoretical guarantees may be obtained. We leave to the appendix detailed theoretical arguments justifying the validity of SIMPL as well as its con- nection to EM.
Algorithm SIMPL: An algorithm for optimizing tuning curves and latents from behaviour
S NXT
> Spike count matrix Initial latent estimate e.g. measured position of animal
2: (0) DXT
3: procedure SIMPL(s,x())
4: for e t to E do
> Loop for E iterations
5: f(e) <- (uiTuuniggCCrres((
b The "M-step" The "E-step"
6: (c)1) DecoodLLaent(((
7: end for
8: return (1++1,f(F5)
The optimised latent and tuning curves
9: end procedure
4

---

Published as a conference paper at ICLR 2025
3 RESULTS
3.1 CONTINUOUS SYNTHETIC DATA: 2D GRID CELLS
First tested SIMPL realistic navigational task by generating large artificial dataset of spikes from population of N 225 2D grid cells - a type of neuron commonly found in the medial entorhinal cortex (Hafting et al. 2005) in a square environment. All grid cells had maximum firing rate of 10 Hz and were arranged into three discrete modules, 75 cells per module, of increasing grid scale from 0.3-0.8 m (Fig. 2 A latent trajectory, x", was then generated by simulating agent moving around the environment for 1 hour under smooth continuous random motion model. Data was sampled 800,000 spikes). All data was generated using the Rat InABox package (George et al. 2024a).
rate of 10 Hz giving a total of
36, 000 time bins (~
The initial trajectory, x ("),, was generated by adding smooth Gaussian noise to the latent such that, on average, the true latent and initial condition differed by 20 em (Fig. 2. top panel). This discrepancy models the agent's 's internal position uncertainty and/or a measurement error. It sufficed to obscure almost all structure from the initial tuning curves f(0) (x) (Fig. 2 top). To assess performance track to the log-likelihood of training and test spikes (see Appendix C.55 for how we partition the dataset). We also calculate the error between the true and latent trajectory the epoch-to-epode
change the tuning curves and the negative entropy (hereon called "spatial info"" of the normalized tuning curves as a measure of how spatially informative they are (Fig. 2t).
SIMPL. was then run for epochs (total compute time 39.8 CPU-secs on a consumer grade laptop). The true latent trajectory and receptive fields were recovered almost perfectly and the log- leellikeldoodd of both train and test spikes rapidly approached the ceiling performance with negligible overfitting. As expected, SIMPL performs better on larger datasets, Fig. 2e. however performance remains good even with substantially smaller datasets (e.g. 50 cells for duration of 5 minutes). We also swept across the velocity and kernel bandwidth hyperparameters (v, o) and found SIMPL was surprisingly robust to changes in these hyperparameters within reasonable limits (see D.2).
(a) latent trajectory, (t)
estimate, true,
(b) tuning curves, (x))
e 10
pos
time [min] Atfl/Ae
(d) log- likelihood error
spatial info.
oeror [om]a0
(c)
true tuning curves, fxx())
-0.05
(mins)
train
80
test
0.23
Figure 2: Results on synthetic 2D grid cell dataset. (a) Estimated latent trajectories (epochs 0, 1 and 10). Initial conditions generated from the true latent (black) by the addition of slow Gaussian noise. Shaded zones show discrepancy between the true and estimated latent. (b) Tuning curve estimates for exemplar grid cells. (c) Ground truth tuning curves. (d) Performance metrics: Left: lieliioid of the training and test spikes (averaged per time step, dotted line shows ceiling performance on a model initialised with the true latent). Middle-left: Euclidean distance between the true and estimated latent trajectories (averaged per time step). Middle-right: Epoch-to-epoch change in the tuning curves showing they stabilise over iteration. Right: Cell spatial information. Violin plots, where shown, give distributions across all neurons. (e) sweep over number of cells and the duration of the trajectory.
Finally, despite having an implicit prior for temporally smooth latent dynamics, further synthetic analysis revealed SIMPL is still able recover discontinuous latent trajectories (for example those
5

---

Published as a conference paper at ICLR 2025
containing jumpy-like "replay" events, Appendix D.3 even discrete latents in non- dynamical task akin to a discrete two-alterativv forced choice task (2AFC, see Appendix[D.]]
3.2 HIPPOCAMPAL PLACE CELL DATA
Having confirmed the efficacy of SIMPL on synthetic data, we next tested it on real dataset of hippocampal neurons recorded from rat it foraged in large environment TTanni et al. 2022)
This dataset consists of = 226 neurons recorded over 2 hours, binned at Hz giving = 36, 000
data samples and 700,000 spikes. Many of these cells are place cells OO'Keefe Dostrovsky 1971 which, in large environments, are known to have multiple place fields (Park et al. 2011).
We initialised with the animal's position, measured by LED located between its ears, and optimised for 10 epochs. The log-likelihood of test and train spikes both increased, converging after epochs (Fig. 3b) compute time of <40 CPU-secs. We then analysed the shapes and statistics of the tuning curves: After optimisation, tuning curves were visibly sharper, Fig. [a: previously diffuse place fields contracted (e.g. the third exemplar tuning curve) or split into multiple, smaller fields (second exemplar). Occasionally, new place fields appeared (fourth exemplar) or multiple place fields merged into a single larger field (fifth exemplar). Statistically, tuning curves had significantly more individual place fields (+19%, mean 1.14->1.41 per cell, 0.0035 Mann Whitney U tests), substantially higher maximum firing rates (+45%, median 4.2 4.2->61 Hz, = 9.8 x 10-7) and were more spatially informative (p 0.038). Individual place fields became smaller (-25%, median 0.0037).
0.54 m²) and rounder (+8% median 0.63 0.63->0.68,
To ensure these observed changes weren't merely an artefact of the optimisation procedure generated a control dataset by resampling spikes from the behaviour- tuning curves, Scom ~ p(-Ix("), p(0)). Control spikes thus had very similar temporal statistics and identical tuning curves to those in the hippocampal dataset but, critically, were generated from known ground truth model exactly equal to their initialization. Thus, any changes in the control tuning curves Post-SIMPL must be artefactual. Indeed, significant changes area (Fig. pe, grey) providing strong evidence the significant changes observed the real data (e.g. the decrease in field area) were genuine, reflecting the true nature of hippocampal tuning curves.
observed besides a slight increase in field
(a) tuning curves, (()(
(b) Log-likelhhodd
0,12
uan
tstt
Control
3:5
-0.15
epoch
10
(c)
no. of place fields
field area, A [m²)
field 'roundness' max. firing rate [Hz]
spatial info.
before after control
6
1.5
15
-8
-
(d) latent trajectory, (()
before (behaviour) after (latent)
(o) x-difference
(f) area vs dist. to wall
3,5
1.5 1.5
opos [m]
A = loxtlott
E
(mm
time [min]
d [m]
Figure 3: Results on a hippocampal place cell dataset collected by Tanni et al. 2022 (a) Exemplar tuning
curves before and after optimization. Automatically identified place field boundaries shown in white. (b)
Log-likelihood of test and train spikes. Control model shown in grey. (c) Statistics analysis of place fields. Violin plots show the distributions over all place fields cells. (d) The final latent trajectory estimated from SIMPL (green) overlaid on top of the measured position the animal (used as initial conditions, yellow). (e) Behavioural discrepancy map: the average discrepancy between the latent and behaviour as function of optimised latent x (10), Overlaid is snippet of the behavioural optimised true latent trajectory. (f) Place field area as function of the distance to the nearest wall.
The optimized latent trajectory x (10) remained highly correlated with behaviour (R² = 0.86, Fig. 3d) occasionally diverging for short periods it "jumped" to and from location, as if the

---

Published as a conference paper at ICLR 2025
animal was mentally teleporting itself (an example is visualized in Fig. 3e). We calculated the difference between the optimised latent and the behaviour at each time point, A = ||
x(11) I|2.
and visualized this as a heat map overlaid onto the latent space (Fig. 3d). We found that the latent discrepancy was minimal near the edges of the environment and peaked near the centre, perhaps because sensory input is scarce in the centre of the environment due fewer visual and tactile cues.
Tanni et al. (2022) observed that the size a place field size increases with its distance to a wall. Our observation-than the latent discrepancy is highest in the centre of the environment- suggests one possible hypothesis: behavioural place fields merely appear larger in the centre of the environment because they blurred by the correspondingly larger latent discrepaney. true, this trend should weaken after optimisation, once the "true" latent has been found.
To this we plotted field size against distance-to-war (Fig. 3); op- timized fields, like behavioural fields, were small very near to the walls and grew with distance (replicating the result of Tanni et al. ), but this correspondence stopped after 0.5 beyond which the optimized place fields size grew more weakly with tottace-toowall. This supports hypothesis, suggesting substamtial fraction of the correlation between size and distance isn't a fundamental feature of the neural tuning curves but an artefactual distortion in the tuning curves, something which can be cor- rected for using SIMPL.
(a)
(b) -0.514 log- likelihood
poooooooo SIMPL4D (hand position velocity)
SIMPL20
posilion)
(test)
SIMPL2D (hand velocity)
-0.528
epoch
(c)
behaviour (hand position only)
SIMPL
behaviour (hand velocity only)
SIMPL
sod now
earr. 0.59
corr
0.41
time
time [sec]
3.3 SOMATOSENSORY
CORTEX DATA DURING A HAND-REACHING TASK
behaviour (hand position... .ww.............................
4D
behaviour hand velocity)
SIMPL
To test SIMPL beyond naviga- tional/hippocampal datasets we ran on macaque somatosensory cor- tex dataset Chowdhury et al. (2020). During this recording monkey made series of reaches to target in one of 8 directions, 4 On half of the als the reach was "active" whereby the monkey moved the manipulan- dum towards the target by itself. On the other half, the reach was "pas- sive". whereby the monkey's 's hand was bumped in the direction of one of the targets by force applied to the manipulandum, forcing the monkey correct and return the cursor the centre. We binned the data (N = 65 neurons, 37 mins, 10 spikes) at 20 Hz and ran SIMPL models on its tirety (i.e. active and passive reaches, well as the inter-trial intervals) for 10 epochs.
CEBRA 4D
MAW Way corr. 0.74
ww mnll carr. @ 0.5T
law Manp many why War Arod
time [sec]
time [sec]
Figure SIMPL applied somatosensory cortex data. macaque performs centre-out reaches;
65 somatosen- sory neurons are recorded. (b) bog-likelihood curves for the three SIMPL models in pancls c-c. (c) SIMPL trained with 2D latent initialised from hand position. Top-left: raw behaviour, averaged across trials aligned to movement onset; top-right: after SIMPL. Middle: 40 s of behaviour (yellow) and latent (green). Bottom: exemplar tuning curves before and after SIMPL. (d) As in c. but initialised with hand velocity. (e) As in c, but with 4D latent ini- tialised to hand- position (dims and 2) and velocity (dims and 4). Inset: 2D visualisation of 4D latent embedding from CEBRA trained on hand position, adapted from Schneider et al. 2023)
First SIMPL was with 2D latent initialised to the monkeys measured and y-hand position (Fig. 4p). Afterwards, the latent trajectory- here averaged across trials with the same direction, aligned to movement onset- had diverged from, but remained correlated with, initial hand-position

---

Published as a conference paper at ICLR 2025
(correlation 0.59). Despite improvement in likelihood over the behavioural initialisation, latent trajectories for distinct directions substantially overlapped with one another, indicating an cient dimensionality to capture the full complexity of the data. A similar result was obtained when initialising to hand-velocity (Fig. 41).
insuffi-
We then trained SIMPL with 4D latent space. Two of the dimensions were initialised with hand position and the other two with hand velocity. This model performed better than either 2D model, converging a higher likelihood. The latent dimensions initialised hand-position remained highly correlated with hand-position (corr. = 0.74) after optimisation as did the velocity dimensions (corr. 0.57). The latent trajectory was also more structured, with distinct and less overlapping motifs for each trial type. We visualised two-dimensional slices of the four-dimensional tuning curves for each neuron and found that they had well-defined receptive fields, similar to place fields in the hippocampus, which were visibly sharper after optimisation. 'hese results suggest that the somatosensory cortex neurons encode a complex and high-dimensional latent, closely correlated to hand position and velocity, which can be partially recovered by SIMPL.
3.4 PERFORMANCE IS IMPROVED BY INITIALISING AT BEHAVIOUR
Latent variable models trained with EM can experience two issues that usually complicate the scien- tific interpretability their results. The first concerns the quality of the solution; does the algorithm converge a good model of the data which predicts the spikes well? The second issue concerns identifioiility; even if the recovered latent trajectory and tuning curves (f(e), x(+) are of high qual- ity, they may differ from the true ones (f x*) by some invertible "warp" & in a way that does not affect the overall goodness-of-fn of the model. These warps could include innocuous rotations and symmteries or, more concerningly the exact structure of the tuning curve is a quantity of interest, stretches or fragmentations. Here we show that behavioural initialisation drastically minimises the severity of both of these issues for SIMPL.
(a) ground truth
(b) initialised: noisy ground truth (fig. 3)
(c) randomly
exemplar
epoch 0
epoch
epoch 10
epoch 10
tuning
curve
(x((o)
0.41
0.83
0.98
00.7
(x;
warp maps
reference
warp distance
0.126
0.050
0.039
0.898
Figure 5: Latent manifold analysis: (Top) Examplar tuning curve in (a) the ground truth latent space, (b) the latent space discovered by -nriroorrlll-iittillie------- discovered by STMPL initialised with random latent trajectory. Inset scatter plots show the and predicted firing rates of all neurons across all times as well as their correlation values ("accurate" models have higher correlations). (Bottom) The warp mappings from each latent space to the "closest" location in ground truth as measured by the distance between the tuning curves population vectors.
after 0, and 10 epochs and (c) the latent space
To do so, we first assess the absolute goodness-of-fit of SIMPL by computing, for neurons, the correlation between the estimated instantaneous firing rates f. (x(*)) quantity invariant warping) and the true firing rates f (xī). SIMPL converges to a highly accurate model (r=0.98)
under behavioural initialization, but less accurate, though still quite accurate, model (r 0.87) when initialised with a random trajectory uncorrelated to the true latent. Next, we estimate, quantify and visualize the warp map é between SIMPL'S 's estimates (f(e) x(e) and the ground truth (f* x*). We obtain this by finding, for every location in the warped space, the position in the true latent space where the tuning curves most similar (o(x) arg miny ||f" (y) f(() (x)]l2). We then quantify the "warpness" " of this mapping as the average distance between x and Φ(x) across the environment, normalized by its characteristic length scale (1 m). This warp-distance should be for totally
8

---

Published as a conference paper at ICLR 2025
warped models and (1)) for heavy warps. In addition to perfectly fitting the data, the solution found by SIMPL under behavioural initialization is minimally warped (warp dist 0.050). In contrast, the good (but imperfect) solution found by SIMPL under random initialization is very heavily warped (warp dist. = 0.498) in fragmented manner. These results (Fig. 5) strongly motivate the use of behavioural initializations in latent variable models an effective means encourage convergence towards latent spaces which are both accurate and un-warped with respect to the ground truth.
3.5 BENCHMARKING SIMPL AGAINST EXISTING TECHNIQUES
We compared SIMPL to four popular methods for latent variable extraction: Vi-VAE (Zhou & Wei 2020), CEBRA (Schneider al. 2023) (which use neural network function approximators), GPLVM (Lawrence 2003) and GPDM (Wang et al.) 2005) (which use Gaussian processes). Cru- cially, and like SIMPL, none of these methods make restrictive linear assumptions about the struc- ture of the tuning curves.
To match SIMPL,, we initialise the latent variable estimates of GPLVM and GPDM to be- haviour (pi-VAE and CEBRA handle behaviour natively by using it to condition a prior over the latent or as contrastive label). All models were trained for their default number of it- erations/epochs. After training we aligned the discovered latents to behaviour and visualised them on top of the ground truth (Fig. (+). All models successfully uncovered a latent tra- jectory closer to the ground truth than behaviour (Fig. GD). SIMPL performed better than
the other models, achieving a final error of 4.2 cm, half that of VAE (8.4 cm).
We posit that Pi-VAE, CEBRA and GPLVM may suffer from the lack of an explicit dy- namical systems component their generative models while GPDM may suffer from the data- subsampling we were required to do to cap the
(8) compute time
(b) latent error
(d)
Ground
SIMPL 140 secs
4.2 cm
behaviour
pi-VAE
10.4 mins
8.4 cm
23.7
8.2
GPLVM
11.6
GPDM
ming
training time to less than two-hours. SIMPL converged in 40 seconds, over 15 times quicker than the next fastest (pi-VAE, 10.4 minutes, Fig. (p)) Except for GPDM, which required GPU, all techniques were run and timed CPU. Only SIMPL was able to recover sharp and accurate grid fields close to the ground
[CPU-mins]
70
ccm]
(c)
ground truth
SIMPL
DA-VEE
truth.
CEBRA
4 RELATED WORK
GPI VM
pos [m]
Probabilistic inference in neural data modu- lated by latent variables has been a major topic
GPDM
time [min]
of study for decades
see, c.g. Tipping
Figure 6: Comparison pi-VAE, CEBRA, GPLVM and
& Bishopp | Yu et al. (2006) 2008b a);
GPDM on the synthetic grid cell dataset. (a) Compute time. (b) Final error in the latent, (c) Alignment of the discovered latent to the ground truth. (d) Exemplar tuning curves constructed using bernel-based estimation the latent (i.e. an "M-step").
Macke al. (201); Mangion et al. (2011)
Park et al. ) Gao et al. (201)) Hernan-
dez et al. (2018): Dong et (2020) Zhou &
Weil (2020): Gondur et al. (023)) Bjerke et al.
however not all methods were designed for the kind of data considered in this work. Many methods contain model complex latent space dynamics but combine these with simplistic tuning curves which restrict firing rates to (exponential- Jlineer functions of the latent (Smith
2023
Brown,
2003 Yu et al. [0088a] Macke et al.) 2011 Duncker et al. 2019: Linderman et al. [016]] Pandari-
nath et al.] 2018, Zoltowski et al. 2020 Sani et al.] 2021 Hurwitz al., 2021 Kim et al. 2021
Gondur et al.) (2023) so cannot interpretably account for the representations (place cells, grid cells)
considered here. Other methods do not/cannot use behaviour to aid latent discovery (Gao et al.
2016 Nam 2015 Hernandez et al.] 2018 Gondur al., 2023 Bjerke et al.) (023) instead taking a fully "unsupervised approach (meaning they can be applied to spike data without an obvious behavioural correlate) at the expense of complexity and identifiability.
Algorithms that both don't restrict to simplistic linear tuming curves and exploit behaviour form a small set of relevant alternatives SIMPI.. Behaviour-informed latent discovery tools have become

---

Published as a conference paper at ICLR 2025
popular in recent years due to the explosion of large neural datasets taken from behaving animals and the observation that behaviour can explain substantial variance in the neural dynamics.
Gaussian process latent variable models (GPLVMs), Lawrence (2003): Wang et al. (2005) form a family of methods that learn smooth, non-linear tuning curves by placing GP priors on them and performing approximate marginal log-likelihood optimisation the latent variable. Popular implementations leave the initial condition of this optimisation user-defined and therefore compat- ible with the behaviour-informed initialisation used here. However, most such models were intro- duced outside of the neuroscience literature thus use Gaussian (instead of Poisson) emission models (Lawrence 2003) Wang et al. 2005; Jensen al.) 2020). or do not make smoothness assumptions on the latent trajectory (Jensen et al.) [020 Lawrence| 2003). P-GPLVM, which employs Pois- son emissions and GP prior on the latent trajectory. is an exception, but its cubic scaling with time points makes it impractical for hour-long datasets. In contrast, available GPLVM implementa- (inghhm et al., 2018) use inducing point approximations to achieve linear time complexity.
CEBRA (Schneider et al. 2023 learns a deterministic neural network mapping from spikes to la- tents using behaviour- or ume-guided contrastive learning. Unlike most methods, CEBRA does not natively learn generative model nor tuning curves, which are of primary interest in our setting. CE- BRA also treats each data point independently instead of modelling whole-traeectories preventing it from taking advantage of the temporal smoothness inherent in many underlying latent codes.
pi- VAEV (Zhou & Wei 2020) uses a variational autoencoder (Kingma & Welling, 2014) to infer the latent trajectories and learn tuning curves using neural network function approximators. pi- VAE places a learnable prior, conditioned on behaviour, to the latent variable in order to obtain a model with provable identifiability properties. However, Pi-VAE suffers from the same limitation as CEBRA in that it treats each data point as an i.i.d observation instead of a part of a whole trajectory.
The properties of large scale neural datasets suggest five desiderata on the algorithms used to analyse them. These are (1) the absence of restrictive tuning curve assumptions, (2) modelling smooth latent dynamics, (3) the presence of spiking component (e.g. Poisson emissions), (4) the ability to exploit behaviour (including as an initial condition) and (5) scalability to large datasets. None of the methods described in our literature review satisfy all five desiderata. In Appendix we provide a table comparing all methods discussed in this section and more with respect to these desiderata.
5 DISCUSSION
We introduced SIMPL, tool for optimizing tuning curves and latent trajectories using technique which refines estimates obtained from behaviour. It hinges on two well-established sub-routines tuning curve fitting and decoding which are widely used by both experimentalists and theorists for analysing neural data, By presenting SIMPL as an iterative application of these techniques, we aim to make latent variable modelling more accessible to the neuroscience community.
SIMPL could be seen as an instance of a broader class latent optimization algorithms. In principle any curve fitting procedure and any decoder (which uses those tuning curves) could be coupled into a candidate algorithm for optimizing latents from neural data. Our specific design choices, while attractive due to their conceptual and computational simplicity, will come with limitations. For example, we predict SIMPL's kernel-based estimator won't scale well to very high dimensional latent spaces (Györf al. 2006) where parametric models, e.g. a neural networks, are known to perform better (Bach) 2017), potentially at the cost of compute time.
Our synthetic analysis focused on settings where behaviour and the true latent differed only in an unbiased manner. It would be interesting to determine if SIMPL's performance extends to more complex perturbations. Fast, non-local and asymmetric perturbations are common in the brain; for instance "replay" (Carr et al. 2011) where the latent jumps to another location in the environment. Likewise, during theta sequences (equrr agent. This forward-biased discrepancy could theoretically induce a backward-biand havioural place fields, even if the true tuning curves remain unskewed. If this the case, proper latent dynamical analysis-via tools like SIMPL place field tuning curves (Stachenfeld et al, [017]] Fang et al. 2023) Bono et al. 2023 George et al. 2023), similar to how it reduced the asymmetry in place field sizes further from walls (Fig. 3).
Maurer et al. 2006). the encoded latent moves away from the
skew in be-
help reinterpret the predictive nature of
10

---

Published as a conference paper at ICLR 2025
REFERENCES
Afsheen Afshar, Gopal Santhanam, M Yu Byron, Stephen Ryu, Maneesh Sahani, and Krishna V Shenoy. Single-trial neural correlates of arm movement preparation. Neuron, 2011.
Francis Bach. Breaking the curse of dimensionality with convex neural networks. Journal of Ma- chine Learning Research, (19):1-53, 2017.
Patrick Billingsley. Statistical methods in markov chains. The annals of mathematical statisties, 1961.
Eli Bingham, Jonathan P. Chen, Martin Jankowiak, Fritz Obermeyer, Neeraj Pradhan, Theofanis Karaletsos, Rohit Singh, Paul Szerlip, Paul Horsfall, and Noah D. Goodman. Pyro: Deep Univer- sal Probabilistic Programming. Journal of Machine Learning Research, 2018.
Martin Bjerke, Lukas Schott, Kristopher T Jensen, Claudia Battistin, David A. Klindt, and Ben- jamin Adric Dunn. Understanding neural coding on latent manifolds by sharing features and di- viding ensembles. In The Eleventh International Conference on Learning Representations, 2023.
URL httpttss///peeereeeess net wawpeee//oouuueeeeeeeee wwwummidd-uuuuOUUe
Rafal Bogacz, Eric Brown, Jeff Moehlis, Philip Holmes, and Jonathan D Cohen. The physics of optimal decision making: formal analysis of models of performance two-alternative forced- choice tasks. Psychological review, 2006.
Jacopo Bono, Sara Zannone, Victor Pedrosa, and Claudia Clopath. Learning predictive cognitive maps with spiking neurons during behavior and replays. eLife, pp. e80671, 2023.
James Bradbury, Roy Frostig, Peter Hawkins, Matthew James Johnson, Chris Leary, Dougal Maclaurin, George Necula, Adam Paszke, Jake VanderPlas, Skye Wanderman-Milee and Qiao Zhang. JAX: composable transformations of Python+NumPy programs. 2018. URL http:
//github.com ax ml/jax
Ralph A Bradley and John J Gart. The asymptotic properties of ml estimators when sampling from associated populations. Biometrika, 1962.
Federico Carnevale, Victor de Lafuente, Ranulfo Romo, Omri Barak, and Néstor Parga. Dynamic control of response criterion in premotor cortex during perceptual detection under temporal un- certainty, Neuron, 2015.
Margaret F Carr, Shantanu Jadhav, and Loren Frank. Hippocampal replay the awake state: a potential substrate for memory consolidation and retrieval. Nature neuroscience, 2011.
Raeed H Chowdhury, Joshua Glaser, and Lee Miller. Area 2 of primary somatosensory cortex encodes kinematics of the whole arm. eLife, pp. e48198, 2020.
John Cunningham and Byron M Yu. Dimensionality reduction for large-scale neural recordings. Nature nearoscience, 2014.
Arthur P Dempster, Nan M Laird, and Donald B Rubin. Maximum likelihood from incomplete data via the em algorithm. Journal of the royal statistical society: series B (methodological), 1977.
Eric L Denovellis, Anna K Gillespie, Michael E Coulter, Marielena Sosa, Jason E Chung, Uri T Eden, and Loren M Frank. Hippocampal replay of experience at real-world speeds. eLife, pp- e64505, 2021.
Zhe Dong, Bryan Seybold, Kevin Murphy, and Hung Bui. Collapsed amortized variational inference for switching nonlinear dynamical systems. In International Conference on Machine Learning, PMLR, 2020.
pp.
Lea Duncker, Gergo Bohner, Julien Boussard, and Maneesh Sahani. Learning interpretable continuous-tin time models of latent stochastic dynamical systems. In International conference on machine learning. PMLR, 2019.
11

---

Published as a conference paper at ICLR 2025
Ching Fang, Dmitriy Aronov, LF Abbott, and Emily L Mackevicius. Neural learning rules for generating flexible predictions and computing the successor representation. eLife, pp. e80680, 2023.
André A Fenton and Robert U Muller. Place cell discharge is extremely variable during individual passes of the rat through the firing field. Proceedings of the National Academy of Sciences, 1998.
Ronald Aylmer Fisher. Theory of statistical estimation. In Mathematical proceedings of the Cam- bridge philosophical society. Cambridge University Press, 1925.
Yuanjun Gao, Evan W Archer, Liam Paninski, and John P Cunningham. Linear dynamical neural population models through nonlinear embeddings. Advances in Neural Information Processing Systems, 2016.
Tom M George, William de Cothi, Kimberly L Stachenfeld, and Caswell Barry. Rapid learning of predictive maps with stdp and theta phase precession. eLife, pp. e80663, 2023.
Tom M George, Mehul Rastogi, William de Cothi, Claudia Clopath, Kimberly Stachenfeld, and Caswell Barry. Ratinabox, toolkit for modelling locomotion and neuronal activity in continuous environments. eLife, 2024a.
Tom M George, Kimberly L Stachenfeld, Caswell Barry, Claudia Clopath, and Tomoki Fukai. A generative model of the hippocampal formation trained with theta driven local learning rules. Advances in Neural Information Processing Systems, 2024b.
Apostolos P Georgopoulos, Andrew Schwartz, and Ronald E Kettner. Neuronal population coding of movement direction. Science, 233(4771):1416-144) 1986.
Joshua Glaser, Ari S Benjamin, Raeed H Chowdhury, Matthew G Perich, Lee E Miller, and Kon- rad Kording. Machine learning for neural decoding. eNeuro, (4), 2020.
Rabia Gondur, Usama Bin Sikandar, Evan Schaffer, Mikio Christian Aoi, and Stephen L Keeley. Multi-modal gaussian process variational autoencoders for neural and behavioral data. arXiv preprint ar arXiv:2310.03111, 2023.
László Györfi, Michael Kohler, Adam Krzyzak, and Harro Walk. A distribution-free theory of nonparametric regression. Springer Science & Business Media, 2006.
Torkel Hafting, Marianne Fyhn, Sturla Molden, May-Britt Moser, and Edvard I Moser. Microstruc- ture of a spatial map in the entorhinal cortex. Nature, 2005.
Christopher Harvey, Philip Coen, and David W Tank. Choice-specific sequences in parietal cortex during virtual-navigationn decision task. Nature, 2012.
Daniel Hernandez, Antonio Khalil Moretti, Zigiang Wei, Shreya Saxena, John Cunningham, and Liam Paninski. Nonlincar evolution via spatially dependent linear dynamics for electrophysiol- ogy and calcium data. arXiv preprint / 1.0459, 2018.
Øyvind Ame Hoydal, Emilie Ranheim Skytøen, Sebastian Ola Andersson, May-Britt Moser, and Edvard Moser. Object-vector coding in the medial entorhinal cortex. Nature, 2019.
Amber Hu, David Zoltowski, Aditya Nair, David Anderson, Lea Duncker, and Scott Linderman. Modeling latent neural dynamics with gaussian process switching linear dynamical systems. arXiv preprint N.Xi::2408.03330 2024.
Cole Hurwitz, Akash Srivastava, Kai Xu, Justin Jude, Matthew Perich, Lee Miller, and Matthias Hennig. Targeted neural dynamical modeling. Advances in Neural Information Processing Sys- tems, 0229379-29392 2021.
Aapo Hyvärinen and Petteri Pajunen. Nonlinear independent component analysis: Existence and uniqueness results. Neural networks, 1999.
Kristopher Jensen, Ta-Chu Kao, Mareo Tripodi, and Guillaume Hennequin. Manifold gplvms for discovering non-cuclidean latent structure in neural data. Advances in Neural Information Pro- cessing Systems, 2020.
12

---

Published as a conference paper at ICLR 2025
James Jun, Nicholas A Steinmetz, Joshua H. Siegle, Daniel Denman, Marius Bauza, Brian Barbar- its, Albert Lee, Costas A Anastassiou, Alexandru Andrei, Cağatay Aydın, et al. Fully integrated silicon probes for high-density recording of neural activity, Nature, 2017.
Rudolph Emil Kalman. A new approach to linear filtering and prediction problems. 1960.
Timothy D Kim, Thomas Z Luo, Jonathan W Pillow, and Carlos D Brody. Inferring latent dy- namics underlying neural population activity via neural differential equations. In International Conference Machine Learning, pp. 5551-5561. PMLR, 2021.
Diederik P. Kingma and Max Welling, Auto auto-Encoding Variational Bayes. In 2nd International Conference on Learning Representations, 2014.
Dmitry Kobak, Wieland Brendel, Christos Constantinidis, Claudia E Feierstein, Adam Kepecs, Zachary F Mainen, Xue-Lian Qi, Ranulfo Romo, Naoshige Uchida, and Christian K Machens. Demixed principal component analysis of neural population data. eLife, 2016.
Rahul G Krishnan, Uri Shalit, and David Sontag.
Deep kalman filters. arXiv preprint
arXiv:1511.05/21, 2015.
Neil Lawrence. Gaussian process latent variable models for visualisation of high dimensional data. Advances in neural information processing systems, 16, 2003.
Colin Lever, Stephen Burton, Ali Jeewajee, John OeKeef, and Neil Burgess. Boundary vector cells in the subiculum of the hippocampal formation. Journal of Neuroscience, 2009.
Itay Lieder, Vincent Adam, Or Frenkel, Sagi Jaffe-Dax, Maneesh Sahani, and Merav Ahissar. Per- ceptual bias reveals slow-updating autism and fast-forgetting dyslexia. Nature neuroscience, 2019.
Scott W Linderman, Andrew C Miller, Ryan P Adams, David M Blei, Liam Paninski, and Matthew J Johnson. Recurrent switching linear dynamical systems. arXiv preprint
2016.
Francesco Locatello, Stefan Bauer, Mario Lucic, Gunnar Raetsch, Sylvain Gelly, Bernhard Schölkopf, and Olivier Bachem. Challenging common assumptions in the unsupervised learn- ing of disentangled representations. In International Conference on Machine Learning, 2019.
Ryan Low, Sam Lewallen, Dmitriy Aronov, Rhino Nevers, and David W Tank. Probing variability in cognitive map using manifold inference from neural dynamics. BioRxiv, 2018.
Jakob H Macke, Lars Buesing, John P Cunningham, Byron M Yu, Krishna V Shenoy, and Maneesh Sahani. Empirical models of spiking in neural populations. Advances Neural Information Processing Systems, 2011.
Emily L Mackevicius, Andrew H Bahle, Alex H Williams, Shijie Gu, Natalia I Denisenko, Mark Goldman, and Michale Fee. Unsupervised discovery of temporal sequences in high-ddmensional datasets, with applications neuroscience. el.ife, 2019.
Andrew Zammit Mangion, Ke Yuan, Visakan Kadirkamanathan, Mahesan Niranjan, and Guido San- guinetti. Online variational inference for state-space models with Neural Computation, 2011.
observations.
Valerio Mante, David Sussillo, Krishna V Shenoy, and William T Newsome. Context-dependent computation by recurrent dynamics in prefrontal cortex. Nature, 2013.
Andrew P Maurer, Stephen L Cowen, Sara N Burke, Carol A Barnes, and Bruce McNaughton. Organization of hippocampal cell assemblies based on theta phase precession. Hippocampus, 2006.
Leland McInnes, John Healy, and James Melville. Umap: Uniform manifold approximation and projection for dimension reduction. arXiv preprint 180.vv8802.0:00 2018.
Bruce L McNaughton, Carol A Barnes, and JJEBR O' O'Keefe. The contributions of position, direc- tion, and velocity to single unit activity in the hippocampus of freely-moving rats. Experimental brain research, 1983.
13

---

Published as a conference paper at ICLR 2025
Mayank R Mehta, Carol A Barnes, and Bruce L McNaughton. Experience-dependder asymmetric expansion of hippocampal place fields. Proceedings of the National Academy of Sciences, 1997.
Robert U Muller and John L Kubie. The firing of hippocampal place cells predicts the future position of freely moving rats, Journal of Neuroscience, 1989.
Hooram Nam. Poisson extension of gaussian process factor analysis for modeling spiking neural populations. Master's thesis, Department of Neural Computation and Behaviour, Max Planck Institute for Biological Cybernetics, Sribingen[Goooee Scholar), 2015.
John O' Keeff and Jonathan Dostnovsky. The hippocampus from unit activity in the freely-moving rat. Brain research, 1971.
spatial map: preliminary evidence
Chethan Pandarinath, Daniel J O'Shea Jasmine Collins, Rafal Jozefowicz, Sergey D Stavisky, Jonathan C Kao, Eric M. Trautmann, Matthew T Kaufman, Stephen Ryu, Leigh R Hochberg, et al. Inferring single-trial neural population dynamics using sequential auto-encoders. Nature methods, 2018.
EunHye Park, Dino Dvorak, and André A Fenton. Ensemble place codes in hippocampus: Cal, ca3, and dentate gyrus place cells have multiple place fields in large environments. PloS one, 2011.
Mijung Park, Gergo Bohner, and Jakob H Macke. Unlocking neural population non-stationarities using hierarchical dynamics models. Advances in Neural Information Processing Systems, 2015.
Karl Pearson. Liii. on lines and planes of closest fit to systems of points space. The London, Edinburgh, and Dublin philosophical magazine and journal of science, 2(11):559-572, 1901.
Michael L Platt and Paul W Glimcher. Neural correlates of decision variables in parietal cortex. Nature, 1999.
Herbert E Rauch, F'unng and Charlotte Striebel. Maximum likelihood estimates linear dynamic systems. AIAA journal, 1965.
Honi Sanders, César Rennó-Costa Marco Idiart, and John Lisman. Grid cells and place cells: integrated view of their navigational and memory function. Trends in neurosciences, 2015.
Omid G Sani, Hamidreza Abbaspourazad, Yan Wong, Bijan Pesaran, and Maryam M Shanechi. Modeling behaviorally relevant neural dynamics enabled by preferential subspace identification. Nature Neuroscience, 24(1):140- 10-1149, 2021.
Steffen Schneider, Jin Hwa Lee, and Mackenzie Weygandt Mathis. Learnable latent embeddings for joint behavioural and neural analysis. Nature, 2023.
Anne C Smith and Emery N Brown, Estimating a state-space model from point process observations. Neural Computation, 2003.
Hugo J Spiers and Elcanor A Maguire. Thoughts, behaviour, and brain dynamics during navigation in the real world. Neuroimage, 2006.
Larry Squire, Anna S van der Horst, Susan GR McDuff, Jennifer C Frascino, Ramona O Hopkins, and Kristin N Mauldin. Role of the hippocampus in remembering the past and imagining the future. Proceedings of the National Academy of Sciences, 2010.
Kimberly L Stachenfeld, Matthew Botvinick, and Samuel Gershman. The hippocampus as predictive map. Nature neuroscience, 2017.
Sander Tanni, William De Cothi, and Caswell Barry. State transitions in the statistically stable place cell population correspond to rate of perceptual change. Current Biology, 2022.
Jeffrey S Taube, Robert Muller, and James B Ranck. Hcad-direction cells recorded from the post- subiculum in freely moving rats. description and quantitative analysis. Journal of Neuroscience, 1990.
Michael E Tipping and Christopher M Bishop. Probabilistic principal component analysis. Journal of the Royal Statistical Society Series B: Statistical Methodology, 61(3):61 1-622, 1999.
14

---

Published as a conference paper at ICLR 2025
Laurens Van der Maaten and Geoffrey Hinton. Visualizing data using t-sne. Journal of Machine Learning Research, 2008.
Aad W Van der Vaart. Asymptotic statistics, volume 3. Cambridge university press, 2000.
Jack Wang, Aaron Hertzmann, and David J Fleet. Gaussian process dynamical models. Advances in neural information processing systems, 18, 2005.
James CR Whittington, Timothy Muller, Shirley Mark, Guifen Chen, Caswell Barry, Neil Burgess, and Timothy EJ Behrens. The tolman-cichenbaum machine: unifying space and relational mem- ory through generalization in the hippocampal formation. Cell, 2020.
Alex Williams, Anthony Degleris, Yixin Wang, and Scott Linderman. Point process models for sequence detection in high-dimensional neural spike trains. Advances in Neural Information Processing Systems, 2020.
Brian A Wilt, Laurie D Burns, Eric Tatt Wei Ho, Kunal K Ghosh, Eran A Mukamel, and Mark J Schnitzer. Advances in light microscopy for neuroscience. Annual review of neuroxcience, 2009.
Anqi Wu, Nicholas A Roy, Stephen Keeley, and Jonathan W Pillow, Gaussian process based non- linear latent structure discovery in multivariate spike train data, Advances in Neural Information Processing Systems, 2017.
Byron M Yu, Krishna V Shenoy, and Maneesh Sahani. Expectation propagation for inference in non- linear dynamical models with poisson observations. In 2006 IEEE Nonlinear Statistical Signal Processing Workshop. IEEE, 2006.
Byron M Yu, John P Cunningham, Gopal Santhanam, Stephen Ryu, Krishna V Shenoy, and Ma- neesh Sahani. Gaussian-proces factor analysis for low-dimensional single-trial analysis of neural population activity. Advances in Neural Information Processing Systems, 2008a.
Byron M Yu, John P Cunningham, Krishna V Shenoy, and Maneesh Sahani. Neural decoding of movements: From linear to nonlinear trajectory models. In Neural Information Processing: 14th International Conference, ICONIP 2007, Kitakyushu, Japan, November 13-16, 2007, Revised Selected Papers, Part 14. Springer, 2008b.
Yuan Zhao and Il Memming Park. Variational latent gaussian process for recovering single-trial dynamics from population spike trains. Neural Computation, 2017.
Ding Zhou and Xue-Xin Wei. Learning identifiable and interpretable latent models of high- dimensional neural activity using pi-vae. Advances in Neural Information Processing Systems, 2020.
Petr Znamenskiy and Anthony M Zador. Corticostriatal neurons in auditory cortex drive decisions during auditory discrimination. Nature, 2013.
David Zoltowski. Jonathan Pillow, and Scott Linderman. general recurrent state space framework for modeling neural dynamics during decision-making In International Conference on Machine Learning, pp. 11680-11691 1191.. PMLR, 2020.
15

---

Published as a conference paper at ICLR 2025
Supplementary Material for "SIMPL: Scalable and hassle-free optimisation of neural representations from behaviour"
A BACKGROUND
A.1 EXPECTATION MAXIMIZATION
Expectation Maximization (EM. Dempster al. 1977) is widely used paradigm to perform statis- tical estimation in latent variable models. The goal of EM is to maximise the Free Energy, lower bound on the lig-liilihood log p(s; f) of the data, given by (following the notations of Section 2.1
F(f, == (x)
s; f)] - B,cx) v(x)] < log p(sp f),
where g is some probability distribution on the latent variable x. Importantly F is maximised, and the lower bound becomes "tight", at q" p(x|s; i.e. the posterior distribution of the latent variable given the s and f. Moreover, for a fixed q, the only f-dependent term in F is Exy()) [log p(x, f)]. To maximise g) and thus also increase the log- ligllikelhhood EM produces sequence (e())))20 of parameters p(e) by invoking, at each step (or "epoch") e, two well known subroutines:
E-step: Define g(e) == p(x|s; f f)=1))); compute F > E
Egie) s
M-step: Compute f(e) arg maxe F(f, = arg maxt Egees s(
with the property that p(; f(+) 2 log p(s; f(e-1)) for all e, grounding the use of EM to max- imise the likelihood of the data. In context it is important note that, due Gaussianity assumption, calculating the expectation in the E-step requires estimating the mean (and variance) of the posterior p(x|s: ; f(e which can be treated as a point estimate of the latent trajectory. iee. a "decoding" of the latent from the spikes. Thus, in the context of neural data, EM offers frame- work to both estimate intensity functions via maximum likelihood, and also to decode the variable encoded by the neurons.
Impossibility of Exact EM for Gaussian-Mddlaated Poisson Processes The E-step of the EM algorithm requires computing a function defined as an expectation w.r.t p(x's ss; f(e In the case intractable to compute in closed form, unless
of Hidden Markov Models, such expectations the latent variable x is discrete (i.e. numerical estimation), or both the transition and the emission probabilities are Gaussian (with mean and variance depending linearly our particular case, exact inference in the model described in Section [2. is impossible because the emission probabilities are Poisson with mean given by non-linear function of via each neurons tuning curve.
x, (Rauch et al. (965)). In
In order perform statistical inference for our spike train model and avoid resorting numerical estimation which is computationally expensive SIMPL makes a set of approximations which we detail below. At high level the goal convert the non-linear, non- non-Gaussian spiking observations, into a variable which is linear and Gaussian with resepct to the latent, thus EM can be performed exactly using a Kalman smoother.
A.2 LINEAR GAUSSIAN STATE SPACE MODELS AND KALMAN SMOOTHING
Linear Gaussian State Space Models (LGSSM) are dynamical systems of the form:
2 Fizz + Ct; G ~~ Q.)
(5)
H,2; + òz, 8, ~ R₂).
where E Rd, x R" Q. E Rdxd H. E Rexd and R E m LGSSMs can be used as latent variable models given some observed data x, where z is treated as a latent vari- able. While these models are limited in their expressiveness, their benefits that inference (here, the "E-steps") can be done very efficiently: not only the posterior p(z)...., Zr|x],.... .x) Gaussian distribution (of dimension Td), but all of its marginals and pairwise marginals
16

---

Published as a conference paper at ICLR 2025
p(z, (1: .x..... XT) ). p( Ze, Zt+1 X1. . .xr) (crucially, the only distributions needed for learning the
parameters of LGSSM via EM ) can be computed jointly in OTT) time using an efficient technique
known Kalman Smoothing (Kalman 1960 Rauch et al. 1965).
Such a scaling contrasts with naive numerical binning-basce alternatives for inference continuous, non-Gaussian State Space Models, which require maintaining estimate of each bin vector of size n (no. bins) where n grows exponentially with the dimension of the latent space, as used in e.g. Denovellis et al. [222]] Instead, for LGSSMs, the Gaussianity means only the mean and covariance of the marginal posterior distributions - of size d and d² respectively - need to be stored. This is not memory intensive and, perhaps more importantly, the Kalman Filter proceeds to compute them in combined (TT) time. In our experiments, we found the cost of the Kalman Filter to be negligible relative to the kernel evaluations which are the main computational bottleneck of S SIMPL.
Note from here onwards we will return to using x denote the latent variable in the LGSSM, and X for observations.
B SIMPL AS AN APPROXIMATE EM ALGORRTTH
B.I MLE-BASED APPRONIMATE E-STEP
Instead of (0)) = p(x|s; f()-1)), SIMPL computes an approximation to g(o) s g( p(x[x; ; f(e-1)) where X is the Maximum Likelihood Estimate (MIF) of given the observations and the current tuning curves f(e-1) defined as:
x arg max log p(s f(e-x))(()))) arg max ΣΣ log |
[x, f(e-1))
N
X arg max X
Σ
log p(ss) Xr: 1(e-1)).
As defined, computing the MLE returns point estimate of the true trajectory that led the observed spike train however we seek a posterior. In particular, MLE does not use the prior knowledge encoded by pxx).
To find the approximate posterior we note that, as function of s, X is itself a random variable. In the many neurons limit, and under certain regularity assumptions, the distribution of this random variable converges to Gaussian, fact known as asymptotic normality. In other words; though (conditioned on x) is a non-Gaussian random variable, X (a deterministic function of s) is approxi- mately Gaussian in the many neurons limit and thus satisfied the conditions of the LGSSM.
We restate a formal statement of this asymptotic normality result in the case of independent, but non-identically distributed observations originally established in Bradley & Gart (1962), and formulated using the notations of the model at hand. For simplicity, we will consider the case where only P distinct intensity functions f,, fi,...,fe assumption.
...,,f fp exist, although versions of this result exist without this
Theorem B.1 (Asymptotic Normality of the MLE ). Let x E Rd. Let s (sit...., SNt) be independent random variables with probability densities p(sul x(;f(())) where t(i) E (l...., P} is the index of the intensity function faty that generated the spike train St.. For p E 1...., denote n. the number of times the intensity function fp appeared in the sequence f(())- Assume that the MLE &, exists and it is unique. Then, under mild regularity conditions, we have:
VN (R:
xi)
d ((O,I(x))
N
where I(x;) = Σ Hoss /(s: xx: ) is the Fisher Information matrix of the model
at x 9 means convergence in distribution, and we defined pl == lim N.
ng
'The i.idd established in Fisher 1925
17

---

Published as a conference paper at ICLR 2025
The asymptotic Gaussianity of the MLE in the many neurons limit suggests performing approximate inference a surrogate Hidden Markov Model, with the same transition probabilities p(X 1 (x)) as the original ones, but where the observations
replaced by X. Leveraging Theorem B.1 SIMPL
approximates the emission prohabilities px, xx)) by the Gaussian distribution ((x,, Ex) where
Σ, = (NI();)) x (NI(x:)) -1, the Fisher information of the spikes.
Temporarily ignoring the x; x,-dependence the covariance matrices Σ, and treating them as deter- ministic (discussed below), the variables (xr;X.) then form the following latent variable system with hidden variables X, and observed variables X2 given by:
X
X o &
(6)
X
Nx E.)
This model is precisely instance of Linear Gaussian State Space Models defined in Equation 5 and the four matrices set to:
F =
(constant)
H, = I
(constant)
Q: I
(constant)
R = Σ,
(time-varying).
This correspondence allows STMPL to compute an approximation of the marginal posterior distribu- tions p(x,ls) Ra p(x,) X) using Kalman Smoothing (Kalman 1960: Rauch et al. 1965). Importantly, the MLE estimates X, can be obtained in parallel for all f; the only sequential procedure remaining being the Kalman Smoothing step. The trajectory (()) of SIMPL's E-step then set the mean of gle).
B.2 SPIKE SMOOTHING AS GENERALIZED M-STEP
SIMPL'S M-step computationally cheap and interpretable. However, it differs from the step of the EM algorithm, which recall is given by:
(c)) arg max F(C, arg max E,co
(7)
Below, we reconcile the two approaches by showing that SIMPL M-step can be seen as instances of more general "model fitting" M-step.
To see why, note the objective function of standard M-step equals (up to a constant in the nega- tive KL divergence between the joint distribution qite) (x(s)p(s) (observable through samples) and the model p(s, x; f). Thus, a standard M-step can be understood as minimizing this KL divergence approximately, by replacing the expectation over p(s) by empirical average over the data s, an approximation which is asymptotically consistent in the large number of time steps limit under suitable ergodicity conditions (Billingsley 1961).
Similarly, SIMPL's M-Step also fits the model p(s, x, f) to the "data" distribution g(e) (x)s) However, instead of doing so by minimizing the KL divergence between the data and the model, does so using a kernel-based estimate | Thus, Both M-steps can be understood as having the same goal, simply differing in their solution to solve it. In that sense, SIMPL's 's M-step is indeed a generalized M-step.
C IMPLEMENTATION DETAILS
Below provide some implementation details that were important to maximise the computational efficiency of the method.
'We denote g² (x) by g² (a]s) to highlight the dependence between æ and s. AAdiitionally, replaces the expectation over (x()) (x|s) by one-sample estimate of it through *
18

---

Published as a conference paper at ICLR 2025
C.I MAXIMIZING SIMPL's COMPUTATIONAL EFFICIENCY
C.I.1 COMPUTATIONAL BOTTLENECKS IN SIMPL
A single evaluation of the log- likelihood log (s, x₂) requires evaluating the kernel-based rate map estimates given Equation This takes (T) time since it involves a sum across all timesteps. Moreover, this calculation will be repeated itself 7'times for each step of the Kalman smoother in order to (1) compute the MLEs X3 (which naively require gradient ascent on log p(s:(x()) and (2) evaluate the MLE variance Σg (NI((R)) implementation of SIMPL's E-step would have quadratic ( time complexity, which would be prohibitively slow for long datasets. Moreover, the second-order differentiation needed to compute I(x,) is also computationally expensive (formally, in introduces a large constant factor in front of the (T²)) term).
(NH, p(s(
All in all, an exact
In the next sections, we describe additional approximations which allow SIMPL to estimate the MLE and its variance in (TT) time and without differentiating the maps.
C.1.2 LINEAR-TIME MLE ESTIMATION
Naive gradient-based solution The naive way to calculate the MLE X₂ is evaluate all N tun- ing curves (recall each evaluation costs O(T)) for some location x, use these to establish the log- likelihood log ((s))x), calculate the gradient of this log-likelihood w.r.t x, and then take, for ex- ample, & gradient descent steps to find the MLE. This process is repeated for each timestep t in the
Kalman smoother, leading a quadratic time complexity of (kNTT).
SIMPL's approach To compute the MLE in linear time SIMPL bypasses the need recalculate the tuning curves at each time step by, instead, binning them onto a discretised grid of points once the start of each iteration.
Formally SIMPL computes n evaluations the tuning curves ř ==
(f(g))..... ((..),
on grid of n points g (gi,... . gn). This has time complexity O(NnT). We uniform rectangular grid of points (the smallest rectangle containing the full observed behavioural variable) of spacing dr. For example, in of 50 points (n 2500).
environment with dr 0.02 m, this would yield grid
Then, given f. SIMPL then discretizes the log- likeliioddo functions log p(s//x) over that same grid:
N log pl 81, (1t)1:))
N
fu
T. == log logp(s:(gg)
Σ
log
8tj!
j=l
(8)
N fo + S1j log F.g log sty! Σ j=1
where we noted fo == [ggg)] Finally, given such evaluations, SIMPL sets its approximation of the MLE to be
x == arg max log p(s = arg max
gEg
This way of calculating the MLE has linear time complexity yielding an improvement for kT.
C.1.3 LINEAR-TIME DERIVATIVE-FRES MLE VARIANCE ESTIMATION
HH A similar strategy could be employed to also compute I(X,) (log which appears in Σ. Here H, is the Hessian operator defined as (z))(z) = V=f(z). To do so, one could compute the Hessian of the rate maps and their logarithm on that grid, from which any H. ((og)(s(&.)))() at the grid-point-based MLE obtained above can be evaluated as Σラ= H (f))())) Sty H, (log (g.). This would be linear in time
H. (log p(,, (g.) however, found that differentiating f could be very slow.
Instead SIMPL takes an entirely different approach and produces an estimation of E, by instead estimating the variance of the posterior distribution p(x+ls:) p(x))p(s,))) x,) pSs, x.). The
19

---

Published as a conference paper at ICLR 2025
posterior variance and the MLE variance expected to closely match, discussed in theoret- ical justification above. Moreover, as this posterior available analytically up to the normalizing constant p(s;). its variance can be approximately computed by binning p(x, (s)) onto the same grid g introduced above, yielding the following fast estimator for Σg
Σ R Cov p(xe(se) Dulg (,)(g H.)
Σ giPe Di Dit
p4
(9)
Σ, Pis
where Pis exp((1)) p(s)) Intuitively, this is equivalent to fitting multivariate Gaussian to the binned likelihood map. The covariance matrix this Gaussian then used as an approximation of the MLE variance. We provide theoretical argument justifying the validity of this formula below.
Theoretical Justification Equation 9is justified by the Bernstein Von Mises theorem, which states that the difference between the posterior distribution and the distribution of the MLE vanishes in the many neurons limit. We restate this theorem using the notations of our paper, assuming a unique rate map, and without stating some of the required regularity assumptions for simplicity. We refer the reader to (Van der Vaart 2000 Theorem 10.1, p.141-144) for the full version.
Theorem C.1 (Bernstein-we von Mises), Let x; E Rd. Ler = (s1e...., be iid random vari- ables with probability density p(sr|xt;; f). Assume that the MLE &, exists and it is unique. Then,
under mild regularity conditions, for any prior p on Xt, we have;
[I(xx(ss) - -N,,, (NI(x))) l|
p(s)) 0 N-000
p(s) where + denotes convergence in probability, and ||lv denotes the Total Variation norm on
bounded measures.
From this theorem, thus have that the (random) posterior distribution behaves (in total varia- tion) as a Gaussian whose covariance matrix is precisely the asymptotic variance of the MLE. Note however that convergence in total variation does not a priori imply convergence of variances. Fur- ther work could examine under which assumptions such a convergence of variances may hold. In practice, found that this approximation yielded satisfying trade-off between performance and accuracy.
C.2 ITERATIVE LINEAR REALIGNMENT OF THE TRAJECTORIES
To improve the identifiability properties and the numerical stability of SIMPL, we also transform
the decoded latent trajectory at each iteration using linear mapping which maximally aligns it
with behaviour defined as xx() <- Mx(")) + where M, c arg min |xg (Mx (e) + c)ll-
This approach ensures the scale, orientation and centre of the optimised latent trajectory are tied to behaviour, preventing accumulation of linear shifts/rotations across iterations and allowing us to interpret the latent relative to, and in the same units behaviour. We suspect that performing this alignment on all iterates after the optimisation would yield similar results. Because the trans- formed latent necessarily has similar scale to the behaviour which was used to set the size of the discretised environment - we can reuse the same discrete grid for the latent avoiding the need to rediscretize the environment at each iteration.
C.3 HYPERPARAMETRA SETTINGS
SIMPL has two model hyperparameters:
v: the diffusion rate for Kalman smoothing, which sets a prior over expected velocity of the latent variable. Units in ms-
σ: the bandwidth of the kernel used in the M-step to smooth spikes. Units are in m.
Additionally there are some impleementaiion-secciiee parameters:
dæ: the bin size for the variance estimation of the MLE. Units in
20

---

Published as a conference paper at ICLR 2025
dt: the time step of the discretization of the latent variable. Units in
E: the number of iterations of the EM algorithm.
Finally, in all simulations used fraction of 10% and held out 'speckled" data segments of length I second to evaluate the performance of the model. We provide in Table [ the value of these hyperparameters for the Artificial Grid Cell Dataset and the Real Hippocampal Dataset.
Table I: Hyperparameters settings
Dataset
dr
dt
E
Artificial Grid Cell Dataset (Fig.
2
0.4 ms 1
0.02 m
0.02 m
0.1
10
Real Hippocampal Dataset (Fig.
3
1.0 ms -1
0.1 m
0.04 m
0.2 s
10
Motor task dataset (2D) (Fig.
4 Acd
1.0
0.1
0.02
0.05
10
Motor task dataset (4D) (Fig.
4c)
1.5
0.09
0.1
0.05 s
10
C.4 SYNTHETIC DATA GENERATION WITH THE RATINA BOX PACKAGE
All synthetic grid cell data were generated using the RatInABox package (George al. 2024a
In this model, an agent moves through I m by I m environment following a smooth, continuous random motion policy (details can be found in the original paper)under default parameters whereby the agent's 's mean speed is 0.08 m/s. Specifically, the RatInABox model was used to generate the true latent trajectory. denoted as x". This trajectory was then "noised" to produce the trajectory used the initial condition for SIMPL, denoted x("), which represents the animal's measured position. The noise, or "discrepancy" vector Ax (0) -x*, was generated by sampling veloeity trajectory from a 2D Ornstein-Uhllnbbed process with zero mean and a coherence time scale seconds. This velocity trajectory was then integrated to obtain the discrepancy. Finally, the agent's behaviour (0) was additionally influenced by the same environmental forces implemented in the standard RatInABox model, i.e. the agent smoothly drifts away from walls to avoid crashing. The adjusted that the
scale of the Krnstein-Uhlenber process latent and observed trajectories was 20 cm in pen-space (i.e. away from walls).
discrepancy between the
Grid cells were modelled as the thresholded sum of three cosine plane waves (see RatInABox meth- ods) with width ratio-ratio of field width to inter-field distance-of of 0.55. Each grid cell is as- signed wave direction 0,, gridscale A; and 2D phase offset d, = [6i, Φ2] T Specifically, N=225 grid cells were divided into 3 modules of N=75 grid cells. Within each module all cells had the same grid scale (0.3, 0.5 and 0.8 m) and wave direction (0, 0.1 and 0.2 rad) but random phase offsets, ap- proximately matching grid cells in the brain. Grid cell firing rates are all scaled to maximum of 10 Hz.
Grid cells firing rates are determined by the latent position time f, i.e. f((t) = p0o from which spikes are sampled according to an inhomogeneous Poisson process Nopikes (t, + dt) Poi(fi(() di)
C.5 TEST-TRAIN PARTITIONING
To assess performance we partition the spike data matrix, s, into testing and training sets, Siest Swain- Inference is performed solely the training set and we then track the log-likelihood of data in both sets (Fig. 1 left), e-g. ((+) | E(i,t)-sua log p({)] r(e)). This partitioning requires careful consideration: entire time intervals cannot be withheld for testing without impairing the model's ability to infer the latent over this period. Likewise, entire neurons cannot be withheld without impairing the model's capacity to estimate their tuning curves. Instead, we adopt a speckled train-test mask previously used in latent variable modelling set-ups (Williams et al., 2020) which withholds for testing extended chunks of time bins arranged in an irregular "speckled" pattern across the data matrix (totalling 10% of the data).
21

---

Published as a conference paper at ICLR 2025
C.6 BENCHMARKING DETAILS
In section 3.5 we benchmarked SIMPL against four comparable methods on the synthetic grid cell datasets. For all techniques, tuning curves were visualised in the same way as for SIMPL: by extracting the latent trajectories after optimisation and using kernel smoothing to construct the rate maps, i..., single "M-step." Default parameters were used for all methods. CEBRA constrains its N-dimensional latent space to lie on an (N-1) (e-1)-dimensional hypersphere. We disabled this constraint for the 2D grid cell dataset. For GPLVM, we used a variant that exploits induction points to improve scalability with the amount of data, performed a grid-search across the number induction points, and reported the best result. For GPDM, which does not feature induction points and thus has a cubic complexity w.r.t T, restricted training for this technique to 5 minutes of data (compared to 60 minutes for the other methods) to keep computation times under I hour, For pi VAE, we set the "task variable" Ur to the behaviour z; (0) and the latent εe to be the true position Tą. Additionally, set the distribution p(za) (e) (e.g. the position given behavior mapping) to a Gaussian distribution
centered u zzz (0) and fixed variance.
D ADDITIONAL RESULTS
D.1 TOY MODEL OF A DISCRETE LATENT VARIABLE TASK
Before testing SIMPL on large temporally continuous dataset we constructed a smaller dataset akin to a discrete two-alternative forced choice task (2AFC) (Fig. 7 widely studied decision- making paradigm (Platt& & Glimcher, 1999 Bogacz et al. 2006 Znamenskiy & Zador 2013 Lieder et al. 2019). The true latent states x E {0, 1J are binary and have no temporal structure (here series of random "left" "right" choices (Fig. 7p). This latent state is stochastically encoded by a population of neurons with random tuning curves giving the Bernoulli emission probabilities under each latent state:
subscript indexes trials not time), analagous
(x = [o0 ~ u(0,1) 1) x 0, (fi ~ u(0, 1)
x= 1,
x (ernoulli(0.5) and 815 X (er)oulli([(x()))
Data is then sampled for 7 50 trials and N = 15 neurons as shown in Fig. 7 Initial conditions, x, are generated from the true latent by randomly resampling a fraction of trials p = 0.5 (Fig. [p). This partial resample represents an initial discrepancy between the behavioural measurement and the true internal state of the agent.
We perform inference this dataset using reduced version of the model (SIMPL-R). In the M-step, tuning curves were fitted by calculating the average activity of a neuron across each latent condition (c.g. f(()(x) (x) E. 0xx()())) x) ,(x ((x)") xx), conceptually similar to ker- nel smoothing). For the E-step, each latent was the decoded according to the maximum like- lihood estimate under the observed spikes and tuning curve estimates from the previous epoch:
x; (e+1) = arg maxx E, log x, )(()) (there is no time dependence between latents, thus no Kalman smoothing). This process repeated for 5 epochs and, with high reliability, converged on the true latents after upproximately two (Fig. 77 & d. distributions show repeat for 1000 randomly seeded datasets, dotted lines show ceiling performance on model perfectly initialised with noise- less (0) x"). We repeated this experiment for various values of p: latent recovery was almost perfect when P was small (i.e. when the initial conditions were close to the true latent), dropping off as o approached 1. At p 1 when the conditions recover latent space that is either perfectly correlated or perfectly anti-correlated ("left" ++ "right") with the true latent (Fig. 77 right), isomorphic solution.
completely random, the model was biased to
D.2 HYPERPARAMETER SWEEP
We swept over the two hyperparameters v (the velocity prior) and (the kernel bandwidth) assess how sensitive SIMPL is to these hyperparameters, as shown in Figure [8] For this we used the same synthetic grid cell dataset used in Fig. 2 Notably, SIMPL's performance (measured in terms
22

---

Published as a conference paper at ICLR 2025
(a)
?
data, s p(-xx", f*)
(c) correlation
log-likelihood
final corr.
(d) tuning curves
15
-20
10 100-++ f: 1
neurons
Neurons
460
trials
50
epoch,
epoch,
0
(b) true latent
x*
"measured behaviour" x
"optimised" latent
O0O
OO O0o0
LA LA LR
trials
50
trials
50
P(spike)
Figure 7: A two-alternative forced choice task (2AFC) toy-model. (a) Data generation: Spikes sampled
from simple generative model. For each of T=50 independent trials random binary latent analogous
"left" or "right" choice is encoded by a population N=15 neurons with randomly initialised tuning curves.
(b) Model performance: Starting from a noisy estimate (yellow) of the true latent (black) where a fraction 0.5 of trials are resampled, SIMPL-R recovers the true latent variables (green) with high accuracy. (c)
Left: Correlation between x(e) and x". Middle: Log- Log-likelihood, log p(s[x(e), f(o)). Right: Final correlation between x5) and a function of initialization noise p. Violin plots show distributions over 1000 randomly
seeded datasets, dotted lines show ceiling performance of perfectly initialised model (x =x) (d) Tuning
curves.
of the final error, see panel b) is relatively stable across wide range of hyperparameters: kernel bandwidths between 0.1 cm and 5 em and velocity priors between 0.2 m/s and m/s all yield similar performance. When the tuning curves are confirmed that kernel bandwidth has significant effect on their appearance. Broader kernels give smoother tuning curves eventually blurring the individual grid fields together whilst narrower kernels give sharper tuning curves eventual leading to overfitting where individual spikes are resolved.
(a)
(b)
1.0
40
20
CPU-
o, [ms
[ms
speed prior 0.4
speed priof
0.2-
80
0.06
N
[cm] kemel bandwidth
0.2
0.06
A
[cm] kernel bandwidthh
Figure 8: Performance of SIMPL. on the synthetic grid cell dataset as function of the hyperparameters (speed prior) and a (kernel bandwidth). (a) Tuning curves. (b) Final error between the latent and ground truth (colour) and total compute time (size).
D.3 NON-CONTINUOUS HIPPOCAMPAL REPLAY DATASET
Since SIMPL. places explicit prior on latent trajectories which are smooth and continuous we tested whether it could be used to model a dataset where the latent variable is non-continuous. For this we simulated synthetic "replay" dataset from N 225 small Gaussian place cells. In this dataset the latent variable and behaviour perfectly match except for for regular, brief periods of "replay" where the latent variable jumps to text we found that SIMPL was able to recover the latent variable, capturing (or "decoding") the replay events with high accuracy (Fig. 9), despite its smoothness prior.
location. Using the same hyperparameters in the main
23

---

Published as a conference paper at ICLR 2025
(a)
Ground truth
Epoch (behaviour)
(b)
o [m]
pos [m]
Pop
Ground truth
Epoch 10
time [mins]
Figure 9: A synthetic hippocampal "replay" dataset. (a) One minute of trajectory, x-coordinate in solid line, y-coordinate in dashed. The behaviour (light-green, top panel) is smooth, actually matching the latent most of the time except when the latent takes regular, brief discontinuous jumps reminiscent of hippocampal replay events. After optimisation SIMPL able to recover the latent glark,green, bottom panel) and capture the replay events with high accuracy. (b) Spike raster plots; spikes plotted against the hehaviour, optimised latent and ground truth latent.
D.4 AUTOMATIC PLACE FIELD DETECTION
In Fig. it was shown that the tuning curves of place cells in the hippocampus undergo statistically significant changes when optimised using SIMPL. For this analysis, individual place fields were automatically identified from the binned rate maps isolated regions of elevated activity within cells tuning curve. This was done by thresholding the activity of each neuron at I Hz and identifying contiguous regions of activity with a peak firing rate above Hz and total area less than half that of the full environment, similar to approaches taken in previous work (Tanni et al. 2022).
E SUMMARY TABLE OF RELATED METHODS
Here we summarize some of the most relevant LVM and dimensionality reduction techniques in the context of our five key desiderata described in the related work section. These are:
Complex tuning curves: Does the model learn/infer non-lineat tuning curves as opposed to lineareepoonntiiallineenn lineerrrel- tuning curves.
2. Smooth latent dynamics: Does the model impose smooth temporal dynamics on the latent space (e.g. by assuming linear dynamical system. Gaussian process using an RNN), as opposed to treating each time point independently.
3. Spike-friendly: Was the method designed for spiking data. For probabilistic models, this refers whether the generative noise model is Poisson as opposed to, say, Gaussian.
4. Exploit behaviour: Does/can the model use behaviour (as observation, contrastive loss- target, initialisation, or otherwise) to guide latent discovery.
5. Scalable: Can the model scale to datasets long duration. Specifically, available open- source implementations of the method does training/inferecee have near-linear time com- plexity. Note this does not mean that compute time necessarily fast in an absolute sense, just that scaling is linear.
24

---

Published as a conference paper at ICLR 2025
Model
Complex tuning curves
Smooth latent
Spike- friendly
Exploit be- haviour
Scalable
dynamics
SIMPL (Our method)
Yes
Yes
Yes
GPLVM Lawrence
2003
Yes
N/S
No
Yes
Yes"
P-GPLVM Wu
al.
2017)
Yes
Yes
Yes
Yes
No
M-GPLVM
Jensen
al.
Yes
No
No
No
Yes"
2020
faeLVM
Bjerke
2023
Yes
Yes
Yes
No
Yes
PILDS
Gao et al. 2016
Yes
Yes
Yes
No
Yes
VIND
Hernandez et al.
Yes
Yes
Yes
No
Yes
2018)
pi-VAE Zhou & Wei
2020
Yes
No
Yes
Yes
Yes
CEBRA
Schneider et al.
N/A
No
Yes
Yes
Yes
2023
MIND et al. 2018
Yes
No
No
No
Yes
LFADS
Pandarinath et al.
No
Yes
Yes
2018)
TNDM
(urwwitz et al. 2021
No
Yes
Yes
Yes
GP-SDEs
Duncker al.
No
N/S
N/S
2019)
sSLD
Linderman al.
No
Yes
Yes
2016)
8PSLDS
Hu
2024
No
Yes
Yes
GPDM Wang
2005
Yes
Yes
No
Yes
No
MM-GPVAE Gondur al.
No
Yes
Yes
No
Yes
2023
PSID
Sani
2021
No
Yes
No
Yes
Yes
GPFA
Yu al.
2008a)
No
Yes
No
No
Yes
P-GPFA
Nam 2015)
No
Yes
Yes
No
No
SSMDM
Zoltowski et al.
No
Yes
Yes
Yes
Yes
2020)
PLNDE
Kim et al.
2021
No
Yes
Yes
Yes
Yes
GLDS
Kalman
1960
No
Yes
No
No
Yes
DKF
Krishnan
2015
Yes
Yes
No
No
Yes
PLDS
Macke al.
2011)
No
Yes
Yes
No
Yes
UMAP
McInnes al.
2018
N/A
No
No
No
No
TSNE
Van der Maaten &
N/A
No
No
No
No
Hinton
2008
PPCA Pearson
1901
Tip-
No
No
No
No
Yes
ping & Bishop 1999)
dPCA Kobak et al. 2016
No
No
No
Yes
Yes
table of comparable models and their properties. N/A means the criterion is not applicable to the model. N/S means the criterion is not specified or may be dependent implementation specifies. Techniques bold compared
Table 3:
benchmark dataset in Fig. [6]
#. Scalable if using an implementation making induction point approximations.
Algorithm could be initialised at behaviour.
Assumes neurons in given ensemble have a shared tuning curve structure (c.g. Gaussian) with neuron-
specific transformations (c.g. shift and scale)
25


Reviews:
[{"review_id": "I2jJYiVGKJ", "reviewer": "Reviewer_6vN3", "summary": "The paper introduces SIMPL, an EM-style algorithm for refining latent variables and tuning curves from spiking neural data. The key idea is using behavior as initialization and combining kernel density estimation with Kalman smoothing in an iterative optimization framework. The authors validate their approach on synthetic datasets and real hippocampal recordings, showing improvements in both computational efficiency and biological interpretability compared to CEBRA (a recent popular deep learning approach for learning latent embeddings from neural and behavioural data)", "questions": "- The authors mention improved identifiability of the proposed approach. Is that solely because of the behavioral initialization?\n- The authors say: \"we see SIMPL as a specific instance of a broader class of latent optimization algorithms, ...\" and then move on to the idea of replacing KDE with neural networks. Can the authors elaborate on that?\n- Can the authors provide the intuition for choosing the kernel bandwidth in the KDE step? This seems critical for the method's performance.", "limitations": "", "rating": 6, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "- The paper is a fresh departure from current latent variable models relying on deep neural networks, with a number of practical advantages (e.g., fast optimization, minimal hyperparameters (only 2), simple implementation based on standard components (Kernel Density Estimation and Kalman smoothing))\n- The paper reveals scientific insights through the use of the proposed method (finer structure in place field representations, new interpretation of place field size distribution, relationship between behavioral uncertainty and neural encoding)\n- The authors do a good job presenting the approach clearly and acknowledging the limitations of the proposed method.", "weaknesses": "The proposed approach relies on strong assumptions that seem to limit the application of the method for more complex cases beyond the datasets used here (e.g., settings where behavioral recordings are missing, higher dimensional latent states, latent spaces not dominated by behavior, non-smooth latent trajectories). It would be great to demonstrate the applicability of the approach on commonly used datasets in the current literature on LVMs in neuroscience (e.g., Neural Latents Benchmark).", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "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": "4gMGGHh0G1", "reviewer": "Reviewer_Kcef", "summary": "The paper introduces a new method they call SIMPL, an EM-style algorithm, aiming to recover low-dimensional, time-evolving latent variables from high-dimensional neural activity in large neural datasets. Their  approach fits tuning curves to observed behaviour and iteratively refines these through a two-step process (EM like). The originality of SIMPL lies in its novel combination of expectation-maximization (EM) techniques with behavior-driven initialization, providing a straightforward yet effective approach to refine neural representations. It’s main advantage is that it is a scalable and fast approach compared to existing popular methods like CEBRA.", "questions": "How does the method perform with smaller neural populations? Can SIMPL perform well when there are limited neurons recorded. Some downsampling analysis on neural data can be helpful to check if the method is sensitive to neuron number.\n\nHow much data is required to optimize SIMPL? How does the mother deal with short recordings? For example, would a short half hour recording be sufficient to fit the model?\n\nThe author mentioned that they trained CEBRA on the synthetic grid cell data using out-of-the-box hyperparameters training for the default 10000 iterations. How much of the performance inferiority of CEBRA VS SIMPL comes from not carefully searching through the hyperparameter space? Is the comparison fair? Would increasing the iterations improve the performance of CEBRA?\n\n\nThe author mentioned that SIMPL finds a modified latent space with smaller, more numerous, and more\nuniformly-sized place fields, suggesting the brain may encode space with greater resolution than previously thought. How could the author disentangle the possibility that this observation is a result of artifacts by their method versus reflecting true representation feature of the hippocampus?", "limitations": "", "rating": 8, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "SIMPL has only two hyperparameters. Minimal hyperparameters make it practical for experimentalists.\n\nRigorous comparison with existing methods (CEBRA).It outperforms CEBRA (whose latent embedding was noisier and had larger final error) and  is over 30 faster.\n\nThe capability to achieve results over 30 times faster than comparable methods while maintaining accuracy makes SIMPL particularly attractive for large-scale neuroscience research\n\nCareful experiments on synthetic as well as real data.", "weaknesses": "Limited Evaluation on Non-Spatial Tasks: While the paper demonstrates impressive results in the domain of spatial navigation, it would benefit from evaluation on other tasks, such as the ones tested in the CEBRA paper so that the two methods can be directly compared in more diverse settings.\n\nScalability Concerns for High-Dimensional Latents: The authors mention that kernel density estimation (KDE) may not scale well to high-dimensional latent spaces. This potential limitation needs more elaboration. What alternatives might be suitable when dealing with truly high-dimensional data.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "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": "cU5rdDtRSB", "reviewer": "Reviewer_gB66", "summary": "The authors propose a new latent variable model for neuroscience that uses kernel density estimation to learn tuning curves from an unobserved latent variable for population spiking data. The model is fit using a closed form EM algorithm and is scalable to large datasets. The authors validate the model on simulated and real neural data.", "questions": "None", "limitations": "", "rating": 5, "confidence": 4, "soundness": 1, "presentation": 1, "contribution": 1, "strengths": "The validation on synthetic and real neural data is a nice addition. The simplicity of the approach has potentially some appeal to experimentalists who want to avoid some of the more bespoke or complex latent variable models for the neural spiking data.", "weaknesses": "The authors provide a very brief overview a very broad field of latent variable models used in neuroscience to motivate their approach. They claim existing methods do not scale well, or have other shortcomings compared to SIMPL. To robustly demonstrate that this model indeed scales better than most or all of these existing approaches while retaining accurate latent and neural tuning identification would require far more comparisons, discussion and evaluation. \n\nIn particular, the authors model is motivated in a similar way to the GPLVM, which is discussed but not compared to. As the authors point out GP based models often do have issues with scalability. However,  there have been many steps toward improving the scalability of these models in recent years -- note that 1 and 2 uses inducing points, a well-known approach to improve scalability in GP  models. These models however have the same 'tuning curve' interpretations that the authors posit SIMPL has, and so it would be important to directly compare to them for both scalability and neural tuning identification. \n\n1 \"Manifold GPLVMs for discovering non-Euclidean latent structure in neural data\" \n2. Learning interpretable continuous-time models of latent stochastic dynamical systems\n\nThere is also no mention of switching linear dynamical systems models, which often scale better than GP based methods due admitting an EM based approach and utilizing forward-backward passing (see e.g. 3, 4, but there are many others). \n\n3 Bayesian learning and inference in recurrent switching linear dynamical systems.\n4. A general recurrent state space framework for modeling neural dynamics during decision-making.\n\nStill there are others in this space that scale well and admit flexible non-linear latent characterizations. See for example 4, 5. and of course LFADS, as the authors discuss. The authors also mention Pi-VAE (and there are of course now many other VAE based models used in neuroscience), but they don't provide a principled comparison or discussion of these approaches. \n\n4 Inferring Latent Dynamics Underlying Neural Population Activity via Neural Differential Equations\n5  Collapsed amortized variational inference for switching nonlinear dynamical systems.\n\nThe choice of CEBRA as the only benchmark is not well motivated. CEBRA is a fundamentally different model it uses behavioral information for contrastive learning, and it was designed to visualize a behaviorally-informed latent space. In this sense it isn't an unsupervised model used purely on spiking data (like SIMPL as well as the ones cited above, among others) and it does not permit a tuning curve interpretation -- which the authors emphasize is an important component of SIMPL and existing in many of these other LVMs. Comparison to the GPLVM, as well as switching LDS, or other nonlinear dynamical models would be far more appropriate then CEBRA. \n\nIn it's current state, it is unclear exactly how much better SIMPL scales than any of the existing approaches, and how it's latent identification or tuning identification may compare to these many state-of-the-art methods.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "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": "gwP1MFRnij", "reviewer": "Reviewer_5Qja", "summary": "The paper proposes a new method, SIMPL for estimating neural representation regarding behaviors. The method is based on the Expectation and Maximization method where during E-step, a model produces estimates of latent trajectories using Kalman smoothing and during M-step, the model uses KDE for fitting the intensity function.", "questions": "* What the arrow of epoch 1 -> 10 in Figure 3 is different from others? Does it have a specific meaning?\n\n* What does “data now shown in L466 mean? Does it mean that the authors intentionally did not include the results since it is insignificant? If so, please include it in the appendix.", "limitations": "", "rating": 6, "confidence": 4, "soundness": 3, "presentation": 2, "contribution": 2, "strengths": "* The paper is well-written. Especially, the references in the introduction can be a good guide to readers who are not familiar with neuroscience, especially, the hippocampus and MEC although the authors simply listed 10 papers in the actually related work section.\n\n* The figure is generally helpful in understanding the method and the results.\n\n* The proposed method is a very simple E-M-based algorithm but it is more efficient and has better performance than CEBRA in the grid cell simulation with the RatsInABox simulator.", "weaknesses": "1. The author claims general but limited to specific regions of brain place cells and grid cells. I highly recommend authors compare with various datasets such as the macaque dataset and the mouse visual cortex datasets used in Schneider et al. (2023).\n\n2. The authors misuse big-O notation. For example, the authors mentioned O(1 hour, 200 neurons, 10^6 spikes) ~ O(1 min) in L 120 and O(hours) in L219. This causes confusion, especially to computer scientists who have been trained in big-O notation for measuring time or space complexity. If the authors want to use big-O notations, please change them to match the function of neurons/spikes and so on. However, I believe, just changing to plain text might be better.\n\n3. Except for Section 4.2, there are no comparisons with existing methods or ablation studies. Although the results appear promising, it remains challenging to discern whether this is due to the simplicity of the task or the efficacy of the proposed method. Moreover, I am not sure whether CEBRA is the only method to compare with.\nThe authors should provide more details on the training and testing splits. L305–306 and 416–417 lack information: are the splits based on spatial segmentation (using specific parts of the box for training) or trial separation?\n\n4. Although experiments with Tanni et al. (2022) represent the main result in this paper, the authors scarcely describe the task and the significance of each plot (particularly in Figure 6c), simply referring back to the original paper.\n\n5. The authors compare SIMPL with CEBRA, a machine learning-based method running on CPUs. If they provide a runtime comparison showing that SIMPL on a CPU is faster than CEBRA on a GPU, it would further support their efficiency claims.\n\n6. Lines 477–479 are ambiguous as to whether the authors discuss place cell remapping or another concept. If remapping is the intended topic, I recommend explicitly stating this and including a citation to aid readers from the machine learning community who may be unfamiliar with this concept.\n\n7. The writing can be improved a lot:\n\n7.1 It is possible that the readers do not know what the “turning curve” is since it is used in neuroscience literature and ICLR is generally a machine learning community. It is good to define what it is in the introduction.\n\n7.2 It would be better to move the Related Work section to Before Method (standard ML conference styles) or Discussion (some ICML styles). Currently, it disconnects methods and results.\n\n7.3 Please choose between British English and American English. Currently, it is used together (e.g., In L477, optimization, behaviour).\n\n7.4 Consider that readers may not be familiar with “tuning curve,” a term more common in neuroscience literature. Since ICLR has a general machine learning audience, defining it in the introduction would be helpful.\n\n7.5 Moving the Related Work section to a position before the Method section (per standard ML conference styles) or to the Discussion (in line with some ICML styles) may improve flow, as it currently separates the methods and results sections.\n\n7.6 Maintain consistency in English dialect, as the paper currently alternates between British and American English (e.g., \"optimization\" vs. \"behaviour\" in line 477).\n\n7.7 In line 131, specify “Appendix A” instead of “Appendix,” even if there is only one section.\n7.8 Ensure the font for v is consistent in lines 142 and 143.\n\n7.9 Since space permits, it may be preferable to adjust the placement of Figure 5 so that it doesn’t span pages 7 and 8, avoiding the empty space near lines 399–405.\n\n7.10 IN L816, 1, …, P => \\{ 1, …, P \\}.\n\n7.11 Correct citation formatting:\n* In line 116, “Tanni et al. (2022)” should use citep.\n* In lines 486, 493, and 790, change citep to citet.\n* In line 522, switch citet to citep.\n\n7.12 It would also be beneficial to unify the reference formatting:\n* Journals like Nature and Nature Neuroscience do not list volume numbers, while JMLR and most eLife papers do (except in lines 583–585). Additionally, both “elife” and “Elife” are used inconsistently, as are “Nature” and “nature.”\n* Clarify “cite” with a red background in L608.\n* While “ICLR” is abbreviated, other conferences like NeurIPS and ICML are listed in full; consider unifying this approach.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "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:

Decision: Accept (Poster)