High-dimensional neuronal activity from low-dimensional latent dynamics: a solvable model
Authors: Valentin Schmutz
Ali Haydaroğlu
Shuqi Wang
Yixiao Feng
Matteo Carandini
Kenneth D. Harris
Venue: NeurIPS 2025 oral ()
Keywords: recurrent neural networks
neuronal recordings
visual cortex
latent variable models
PCA
eigenvalue decay
mean-field limit
Primary area: neuroscience_and_cognitive_science
TL;DR: We show that high-dimensional neural activity can arise from low-dimensional latent dynamics, both in RNNs and in the brain.

Abstract: Computation in recurrent networks of neurons has been hypothesized to occur at the level of low-dimensional latent dynamics, both in artificial systems and in the brain. This hypothesis seems at odds with evidence from large-scale neuronal recordings in mice showing that neuronal population activity is high-dimensional. To demonstrate that low-dimensional latent dynamics and high-dimensional activity can be two sides of the same coin, we present an analytically solvable recurrent neural network (RNN) model whose dynamics can be exactly reduced to a low-dimensional dynamical system, but generates an activity manifold that has a high linear embedding dimension. This raises the question: Do low-dimensional latents explain the high-dimensional activity observed in mouse visual cortex? Spectral theory tells us that the covariance eigenspectrum alone does not allow us to recover the dimensionality of the latents, which can be low or high, when neurons are nonlinear. To address this indeterminacy, we develop Neural Cross-Encoder (NCE), an interpretable, nonlinear latent variable modeling method for neuronal recordings, and find that high-dimensional neuronal responses to drifting gratings and spontaneous activity in visual cortex can be reduced to low-dimensional latents, while the responses to natural images cannot. We conclude that the high-dimensional activity measured in certain conditions, such as in the absence of a stimulus, is explained by low-dimensional latents that are nonlinearly processed by individual neurons.

High-dimensional neuronal activity from low-dimensional latent dynamics: a solvable model
Valentin Schmutz' University College London WCIE 6BT London, UK v. s.schmutzuul..c ac.uk
Ali Haydaroğlu"
University College London WC1E 6BT London, UK
ali hww.aaooou.uu...... 20@ucl.000cc...mm ac. uk
Shuqi Wang École Polytechnique Fédérale de Lausanne 1015 Lausanne, Switzerland shuqi. wang@epfl. ch
Yixiao Feng
Shanghai Jiao Tong University
200240 Shanghai, China y yf2887
Matteo Carandini University College London WCIE 6BT London, UK m. carandini@ucl..... uk
Kenneth D. Harris
University College London WCIE 6BT London, UK
kenneth, www.ett..arristccl.cco ac.
Abstract
Computation in recurrent networks of neurons has been hypothesized to occur at the level low-dimensional latent dynamics, both artificial systems and the brain. This hypothesis seems at odds with evidence from large-scale neuronal recordings mice showing that neuronal population activity is high-dimensional. To demonstrate that low-dimensional latent dynamics and high-dimensional activity can be two sides of the same coin, we present an analytically solvable recurrent neural network (RNN) model whose dynamics can be exactly reduced to low- dimensional dynamical system, but generates an activity manifold that has high linear embedding dimension. This raises the question: Do low-dimensional latents explain the high-dimensional activity observed in mouse visual cortex? Spectral theory tells us that the covariance eigenspectrum alone does not allow us recover the dimensionality of the latents, which can be low or high, when neurons are nonlinear. To address this indeterminacy, we develop Neural Cross-Eneoder (NCE), interpretable, nonlinear latent variable modeling method for neuronal recordings, and find that high- gigh-dimnnsional neuronal responses to drifting gratings and sponta- neous activity in visual cortex can be reduced to low-dimensional latents, while the responses to natural images cannot. We conclude that the high-dimensional activity measured in certain conditions, such in the absence of stimulus, is explained by now-dimensional latents that are nonlinearly processed by individual neurons.
Introduction
The mammalian cortex comprises large number neurons, which, in principle, should allow it tigh-dimensional neural code to represent sensory, motor, and cognitive information. Nevertheless,
recordings in nonhuman primates 4 have suggested that cortical
"quul contribution
39th Conference on Neural Information Processing Systems (NcurIPS 2025).

---

populations perform computations by approximating low-dimensional dynamical systems [S [5 with neuronal firing rates lying on a low-dimensional "neural manifold" 77]. In support of this hypothesis, low-dimensional dynamics have been inferred from multi-neuronal recordings through wide variety of methods [3]: they spontaneously emerge in recurrent neural networks (RNNs) trained to solve behavioral tasks (2 [24344 ; and they appear in several theoretical models of noise-robust neuronal population dynamics 335-38]. A result that might at first sight challenge the low-dimensional dynamical systems hypothesis is that visual cortical population activity in mice has high linear dimension 33 40] with shared neuronal covariance having a heavy-tailed eigenspectrum (see also
41 and 22 for recordings in cerebellum and across cortex, respectively). In particular, the shared covariance eigenspectrum has a power-law tail with an exponent close to 1 (a z 1.04) [39] for responses to natural images and exponent of R 1.14 for spontaneous activity [40]. Are these two views on the dimensionality of population activity compatible? Namely, can a low-dimensional dynamical system produce neural manifold that has high linear embedding dimension?
Here, we first construct solvable RNN model that reconciles the low- and high- high-dimensional perspectives on population activity by carefully disambiguating the linear dimension of the system before and after the neurons' nonlinearity, which we refer to as the pre- and post-activation dimension, respectively. This dichotomy refines the usual distinction between linear and "intrinsic" dimension [39 43 4.), since the intrinsic dimension of a system is the same before and after any continuous, injective nonlinearity. Using the notions of pre- and post-activation linear dimensions, we show that our RNN can be exactly reduced to a low-dimensional dynamical system in the space of pre- activations, making the pre-activations cow-dimensional. Then, show that these latent dynamics generate high-dimensional post-activation activity that has a power-law covariance eigenspectrum. (In this work, dimension will always refer linear dimension, unless stated otherwise.)
Before analyzing experimental recordings, we revisit the spectral theory of infinite-width neural networks (random feature kernels) 45 to quantitatively relate the pre-activation dimension, the neuronal activation function, and the post-activation covariance eigenspectrum. This three-way relationship tells us that high-dimensional activity is consistent with both low- and high-dimensional pre-activations To uncover the pre-activation dimension of high-dimensional activity in visual cortex, we perform two-photon calcium recordings of tens of thousands of neurons from mouse visual cortex, and infer the pre-activn rreaactivation dimension using the Neural pooss-Encoder (NCE), an interpretable. nonlinear latent variable modeling method which models the activity of each neuron as a simple imnar-nonlinear readout of low-dimensional latents. NCE reveals that both the responses to drifting gratings and spontaneous activity can be well approximated by low-diment dimeeniinio antvaaiiv,, but that responses to natural images cannot. This suggests that the encoding of natural images in visual cortex is already high-dimensional in the space of pre-activations.
2 Solvable RNN Model
To demonstrate how high-dimensional post-activations can arise from low-dimensional pre-activation dynamics, we first present solvable RNN model whose autonomous dynamics is low-dimensional in the space of pre-activations, but nigh-dimensional in the space of post-activations, with the post-activations producing a power-law covariance eigenspectrum.
We consider an RNN consisting of N rate-units (neurons). The pre-activation 2; of neuron evolves according to
1
= + N Wyyb(zy),
where Way denotes the synaptic weight from neuron neuron i, and R - Rzo is nonlinear activation function converting the pre-activations into post-activations (firing rates). To define the weights Wij. we randomly place neurons on a ring [48-50] by assigning to each neuron i an independent and uniformly distributed angle 0; E [0, 2T) (Fig. [4]. The weights Wij are then given by the following shifted cosine function:
Wij := cos(0: 8, - A).
The shift A in Eq. (2) makes the weights asymmetric, with neurons sending their strongest excitatory output to neurons located at an angle A counter-clockwise (Fig. [ To make the model solvable,

---

we assume that the activation function Φ is the Heaviside step function θ i.c.,
1, i 2
((x)
(
0, if <.0. 0.
A
RNN on a ring
B
shifted cosine weights
C
rank-2 factorization
2x
Vr
W
W
U
2m
post-activation eigenspectrum
D
latent dynamics
E
F
10"
norm. eigevaalue
K1
10-5p
a
Kg
D
16
10'
10
10'
10'
rank
Figure 1: Solvable RNN model. (A) Schematic the RNN model. (B) Shifted cosine function defining the Desition-dependen synaptic weights Wij- (C) Factorization of weight matrix W as the product of 2-column matrix U and two-row matrix VT. (D) Vector field of the dynamics of the latent variables, K1 and K2, in the large-network limit. The latent dynamics produces a stable limit cycle on the unit circle. (E) Graphical representation of the RNN's effective dynamics. (F) Post-activation covariance eigenspectrum. The eigenvalues follow a power-law with decay exponent 2. They are normalized such that the largest eigenvalue is The activation function used is the step function defined in Eq. (3). Shown in D and F are theoretical values given by Eqs. (5) and (8) respectively, which closely match simulations of large networks (see Appendix
2.1 Low-dimensionall pre-activation dynamics
The RNN model defined above has low-dimensional pre activation dynamics because the weight matrix W has rank 2. Indeed, using an elementary trigonometric identity P W can be factorized as the outer product W UVT of two N 2-matrices (Fig. (C). with
cos(O sin(O.)
cos(O + +
sin(O +
U
and V
cos(f
sin(
cos(AN A) sin + A)
Then, following Beiran dimensional system describing the dynamics of the latent variables k Ulx, where denotes the pseudoinverse and x the N W-dimensional vector of activvivas (z1...., The dynamics of the latent variables follows
[5T] we can reduce the N-dimensional system, Eq. (T). to
A = - - N V (UR),
"cos(a- (3) cos(a) cos(3) + sin(a) sin(
3

---

where ø is applied element-wise to the N-dimensional vector of pre-activations UR = x. In the
equation above, the vector (UUx = ((x) represents the joint post-activations (firing rates) of the N neurons.
Taking the number of neurons N oc yields neural field limit becomes an integral over the ring,
52
where the sum over neurons
22t v(0) (1))) K)
d0
k =
44)
6=
2T
and v(0) == J (cos(0 A), sin ( + A)) Equation 4 describes the dynamics of the latent variables K. as the solution to 2-dimensional dynamical system whose vector field involves an integral over the "circuit structure" [52] (the ring).
with u(0) :=
Since e is the step function, the integral over the ring in Eq. 14 can be solved, and we obtain the
solvable 2-dimensional dynamical system,
= - + (1
(5)
||+|||
when J: TV2 and A == x/4 (derivation presented in Appendix Aj. Equation 5 generates a stable limit eycle over the unit circle (Fig. [D], which implies that the latent variables k. will eventually rotate on the unit circle indefinitely. In Appendix E, we provide other examples of low-rank RNNs for which the latent dynamics can be expressed in a tractable form similar Eq. (
Neuronal activity, modeled here as the post-activations (UU)), are simple linear-nonlinear readouts of latent variables к (Fig. [E). Hence, we have effectively reduced the dynamics of the RNN, Eq. (T. to a 2-dimensional latent dynamical system. We will say that neuronal activity is a linear- nonlinear function of latent variables if it is given by the composition of a linear mapping (U) and an wlss nondecreasing nonlinear mapping (Ø).
2.2 Post-activations produce a power-law eigenspectrum
Since the dynamics of the latent variables is solvable in the aarge-networf limit, and rotates on the unit circle, we can compute the correlation between the post-activations of two neurons. For any pair of neurons i and j. with positions 0: and f, on the ring, respectively, the correlation of their post-activations Cij is, in the long-recording limit, given by
Cij 2
10 0,1 1,
(6)
where 18, 10,,1,1 - := cos cos-ooossaoo (cos(0)--)))) 0,)) ) is the absolute angle difference 6, and 8, (derivation presented in Appendix B)
We can find the eigenvalue spectrum of post-activations by noting that as the number of recorded neurons M + oo, the eigenvalues of the M A correlation matrix (lM) defined by Eq. (6) converge to the eigenvalues of an integral operator. Since the angles 0, are independently and uniformly sampled on circle [0, 2T), the correlation matrix ((M) is a so-called Euclidean random matrix (53), that is, matrix whose entries are given by the pairwise distances between randomly
sampled points in given space. Writing X (M) >/2
(M) 2 ... 2 X for the ranked eigenvalues of tells us that, as M > co, the scaled eigenvalues
(M)
the matrix C(M), random matrix theory 55
55
(M) /M) M converge (in lg sense) to the eigenvalues of the integral operator, (x(1)/
do'
f(0)
n - 1--000 01) - f(0')
2π
0-th are cosine kemel
The eigenvalues of this integral operator can be computed analytically. In Eq. 7 we have highlighted the presence of the O-th arc-cosine kernel ko(0, := 10 0¹ of Cho and Saul [5] which is
well-known in machine learning, and whose eigenvalues have been computed in [46]. In short, by the rotational invariance of ko. have that, for any positive integer m, the functions 0 -> cos(m)f and
0 + sin(m are orthogonal eigenfunctions of the operator Eq. ( sharing the same eigenvalue,
2
(-1)m
(I 0) cos(mm))d
T² 0=0
nT2
m²

---

Using this result, we obtain that the ranked eigenvalues 22>, 2 of the operator Eq. 7 are given by
An 77 (
4 1)/2] (2 1) Vn E
(8)
that is, eigenvalues come in identical pairs that decay exactly power law with decay exponent = (Fig- TT) Hence, post-activation are high-dimensional in the sense that their covariance eigenspectrum has heavy tail [39]. Although the model presented assumed, for simplicity, a cosine lateral connectivity on the ring. Eq. ( similar results can be derived for more general lateral connectivity; see Appendix D for example.
In summary, this solvable model shows that www-dimeniional dynamics in the space of pre-activations can generate high-dimensioni post-activations. The heavy tail of the covariance eigenspectrum implies that post-activations are confined to any inite-dimensional linear subspace. Formally, the smallest vector space containing the post-activations generated by our model has the same size as the feinite-dimensiona reproducing kernel Hilbert space associated with the kernel kg. We stress that, in this model, the heavy tail of the post-activation eigenspectrum is not due to noise, since we used a deterministic, non-chaotic RNN. Also, all the results presented above remain exact if the rate- units in Eq. ① are replaced by Ahnea-nomlinear-Poisse. neurons, as spike noise cancels out in the limits we consider 38 52
3 Post-activation eigenspectrum depends on pre-activation dimension and activation function
To shed light on the relationship between the post-activation eigenspectrum, pre-activation dimension. and the activation function é. we now turn to more general setup, which allows us to relax some of the strong assumptions of the solvable model (Fig. +). First, we allow the number of latent variables d to be greater than 2, assuming that the latent variables, henceforth denoted by z (instead of k)), are uniformly distributed on the unit sphere gd-1 in Rd We assume that the pre-activations of the network are determined by passing the latent activity through N d feedforward weight matrix U with i..d. standard normal entries. this setup (Fig. [ we call d the pre-activation dimension, as it sets the linear dimensionality of the pre-activations. In the solvable model of Sec. 2.2 for example, the pre-activation dimension was d = 2 (Fig. TE). Finally, we replace the step function, Eq. (3), by the general rectified power activation function
(p.c(x) = [max(0, c)]²,
where the activation parameter p E R220 is a nonnegative real value and the bias E RR. (By convention, (x) == O(x c).)
This setup can be analyzed within the framework of random feature kernels (see 156] Sec. 9.5]). Denoting Pdd1 the uniform probability measure on the sphere Sd-1 let us take T independent latent variable samples 21,....24 Zy from 14d-1, and define the N x T post-activation matrix (N,T) ==
(B..,,U (Uzr)). In the limits N and T -> taken successively, the covariance eigenspectrum of A (N.2) converges (when properly scaled) the eigenvalue spectrum of the integral operator
f(z) + [
(10)
where Kp,add 8d-1 x gd-1 - IR is the random feature kernel
Kp Kp.ca(z,z') == Ee (a..llll))pppll
z)0p. z')]
(11)
(see
46 [6 or Appendix For more details).
Drawing intuition from Fourier analysis, the smoothness of a function (here the kernel) should be related to the decay rate of its Fourier transform (here the eigenspectrum) the smoother the function, the faster the decay rate of its Fourier transform. Known results on the eigenvalues of random feature kernels for the cases p 0 and p 1, with O. confirm this intuition and show how it extends to general integers d [46]. Extrapolating those results to any nonnegative p and any real e, we get the following conjecture.
5

---

A
B
Eigcuspectrum decay rate
C
Effect of the bias
16
simulation thoury
0.0
10"
norm. eigenvalue
c 0.2
decny exposiens
10
10
F2
4
10-3
10-
3
10
10
0.0
1.5 activation perameter p
0.5
1.
10'
10
10'
Tank
Figure 2: Eigenspectrum of the random feature kernel (II). (A) Schematic of the feedforward
setup (in the case of three input variables). (B) Comparison between simulations of the eigenspectrum decay rate a and the theoretical value predicated by Conjecture I for [0, 2] and d 2, 4. Up to finite-size effects (number neurons = 2 10 and number of inputs T = 10*), simulations match the theoretical prediction a = 1 + 2p+11 (In this plot, c = 0.) (C) Simulations of the eigenspectrum for various bias parameter e, while keeping and fixed. (In this plot, p = and d = 2, which gives = 2. Eigenvalues are normalized as in Fig. ?
d-T
Conjecture 1. For any p E R20, e R, and any integer 2 the ranked eigenvalues A1 of the integral operator 10) obey the following power-law decay:
A2 2
with a a= 1=1 2p +1
An a
(12)
d
where an by means lim an/bn C (0, +00).
To the best of our knowledge, Conjecture T is not straightforward consequence of any existing result in theoretical machine learning [47 57 58 or harmonic analysis (59H62), hence our presentation of Eq. (12) as conjecture. Note that when the activation parameter is an integer, Φp.c is p-times weakly differentiable, that is, the first weak derivatives of op.c are all locally integrable. This, and the fact that the bias does not affect the decay rate, suggest further extension of the conjecture to more general activation functions, with replaced by the weak differentiability of the activation function.
We tested Conjecture nnmerrially by performing PCA on large post-activation matrices A (N.T).
The linear and continuous dependence of the decay exponent a the activation parameter p predicted by Eq. (2) was confirmed in simulations (Fig.23) Simulations also confirmed that the bias c of the
activation function does not affect the decay rate (Fig. 2C), fact already mentioned in 63 64]
To summarize, the spectral theory of random feature kernels suggests a three-way relationship between the power-law tail exponent of the post-activation eigenspectrum, the pre-activation dimension, and the activation function. This relationship should hold when we can consider the neurons as linear- nonlinear functions of the latent vector, with weights that vary randomly and independently between neurons. The relationship suggests that, when high-dimensional neuronal activity (modeled here as cost-activations) is observed [39 40), two scenarios possible: high- nigh-dimensional activity could arise from nonlinear transformation of low-dimensional latent states, or it could reflect pre-activations that are already high-dimensional. To distinguish these two scenarios, propose, in what follows, method for inferring the pre-activation dimension of neuronal activity in experimental recordings.
Latent Variable Modeling of Neuronal Recordings
To estimate the input dimensionality of neuronal activity in mouse visual cortex (as defined in Sec. 3) we developed the Neural Cross-Encons (NCE), nonlinear generalization of Reduced Rank Regression. Using NCE, we show that high-dimensional neuronal responses drifting gratings
'The k-th weak derivative a function f R + IR is the defined as the function g E Ti that satisfics fa (())g(x)dr (-1)# Ja (z))) for all ip C (R).

---

are well-approximated by linear-nonlineee nonlinear readout of low dimensional latent variable, whereas responses natural images are not. Finally, we apply NCE high-dimensional spontaneous dynamics in the cortex and find that they are well-approximated by linear-nonlinear readout of low-dimensional latents.
4.1 Experimental Data
We conducted large-scale volumetric two-photon microscopy on awake, adult mice during visual stimulation and spontaneous activity. We targeted primary and higher visual cortices with Light Beads Microscope [65]. and extracted deconvolved activity traces for 19,223 = 2.948 neurons using Suite3D [66] as described in Appendix[H] Recordings were performed (1) responses to 320 full-field drifting grating stimuli with 2-14 repeats each; (2) responses to 1866 natural images with repeats each; (3) spontaneous activity in the absence of stimuli for 10-15 minutes.
three stimulus conditions:
4.2 Neural Cross-Encoder (NCE)
The Neural
(NCE) divides neurons randomly into two sets: source set and target set. It predicts the activity b, of the target set from the source set a, via non-linear readout of set of latents, 2₁ (Fig. A). NCE uses a multi-layer feedforward encoder & that ends in bottleneck layer whose activity 2, = E(a)) represents a low-dimensional latent state estimated from the source neurons. The reason to use this rather than an autoencoder, which predicts one set of neurons from themselves, is to discard variability that is not shared across neurons. The NCE we used here has single power-ReLL output layer, matching the setup of section 3 and a 3-layer encoder allowing flexible estimation of latent variables, so that the number of latent variables can be readily interpreted as the pre-activation dimension. We train NCE with stochastic gradient descent on source-target activity pairs described in Appendix I When all nonlincarities
removed, NCE becomes equivalent to Reduced Rank Regression [67]. When predicting stimulus-driven activity, we pair the activity of source and target neurons on different repeats of the same stimulus, to also discard shared variability that is not related to the stimulus 33
Linear-nonlineat readout. The recorded activity of a set of target neurons at time t. b, E R is modeled weighted sum of the latent variables, z, Rd passed through nondecreasing nonlinearity:
b- (p.epp(U +r.
(13)
Here, Pn.c the rectified power activation function, defined in Eq. with a power parameter p that is constant across neurons, pre-activation bias that varies across neurons encoded by an N-dimensional vector c, and post-activation added bias r to account for non-zero baseline firing rates. The decoder parameters {p. c, r are learned alongside the encoder parameters of E. The fact that the decoder, Eq. (I3), has single-layer crucial as it allows us to interpret the latents (z,) as linear factors of the observed neurons' ire-activations (Uz,). It is this constrained decoder that allows to infer the pre-activation dimension of neuronal activity; in comparison, a multi-layer decoder as used in
13 23 would infer something closer to the intrinsic dimension of neuronal activity, which is not our
goal.
4.3
Results
NCE identifies the latent dimensionality of simulated data. To validate that NCE can identify the pre-activation latent variables. we test it on simulated data generated from the toy model in Sec. 2 with a ReLU readout (p = 1) and d = 2. NCE recovers the true latents up to a scaling and shift (Fig B). Moreover, NCE can explain all of the variance the population with only two pre- pre-activation dimensions, while the corresponding linear model (Reduced Rank Regression) requires more dimensions 3C).
Pre-activation dimension is low for grating responses, high for natural image responses. We next consider the pre-activation dimensional of visual stimulus responses of visual cortex neurons. To ensure that the NCE focused on the stimulus responses, and not correlated ongoing activity such
7

---

A.
B.
C
E())
9,-
Ground Trtth
1.0
Fracton of Varience Expleinei (norm.)
Interred (NCE)
Länear
Pre-Activation Dim (4)
BOtrсe Deurous
pre-activation
target neurons
lstents
Figure Validation of the Neural Cross Encoder (A) Schematic of the NCE. B Simulated pre-activation latents of a two-dimensional toy model (top), and the inferred pre-activation latents using NCE (bottom). C Reconstruction score function of pre-activation dimensionality for Linear (RRR) and NCE models on simulated data.
as spontaneous activity or encoding of movements, the activity of the target cells and the source cells were taken from different repeats of the same stimuli. In the case of drifting gratings, for which we know there is low-dimensional latent variable (the grating orientation), an NCE model with low-dimensional pre-activations accurately predicts neuronal responses (Fig. (Fig.BA, (A,)). NCE requires fewer dimensions (5.5 1.2, mean + std) than the corresponding linear model (13.9 = 3.0) to predict 95% of the explainable variance (defined as the maximum variance explained across all d in both models). On the other hand, NCE models with low pre-activation dimension are not sufficient to predict responses to natural images (Fig. DC.D. requiring 93.9 + 6.0 dimensions to reach the threshold, suggesting that natural images produce high-dimensional representations in the space of pre-activations (but see Limitations below). Linear models only account for smaller fraction of the total variance (Fig. 4DD), and therefore undcrestimate the dimensionality natural image responses (48.0 + 13.1).
Spontaneous activity has low pre-activation dimension. Spontaneous activity is well predicted by NCE with low pre-activation dimension (Fig. 4E.F). Across all recordings, spontaneous activity of 1000 target neurons has an estimated pre-activation dimension d of 7.0 1.0 (mean + std.), somewhat larger than grating responses but substantially lower than natural images responses (Fig. G,ig- G.H The linear model finds similar dimension (7.5 + 2.0). though its performance deteriorates at high pre- activation dimensions due to overfitting, while the NCE performance remains consistent (Fig. 4
These results indicate that visual cortex activity can be modeled as a linear-nonlinear transformation of a latent vector, which is now-dimensional for grating responses and spontaneous activity, but high-dimensional for natural image responses. In the case of grating responses, we find latents that resemble the sine and cosine of the stimulus angle (Fig. 4/A)this this is what one would expect find neurons follow the canonical model of simple cells in visual cortex [68]. On the other hand, during spontaneous activity, the dynamies of the latent variables are correlated with the running speed the mouse, and are perhaps related to its arousal state (see Appendix
5 Summary of Technical Contributions and Previous Works
Latent dynamics of low-rank RNNs
Low-rank RNNs are tractable models of how the brain can
perform computations through low-dimensional population dynamics [36] [5 69 722 particular,
the dynamics of certain low-rank RNNs, in the large-network limit, reduce to that of "effective
circuits", i.e., dynamical systems describing the evolution of the latent variables [51 52 73 A
limitation of these effective dynamical systems is that the expression of the vector field involves an
integral over the distribution of weights (the "circuit structure" (52)). making them somewhat opaque
and costly solve numerically in general. In this work (Sec. 2. and Appendix
we prove that, in

---

A.
Drifting Gratings
C.
Natural Imagpes
EL
Spontaneous
protiction
(8 00:3)
pooilctin
E8 23
PreAActhattin Dis (4
+t
lateato
www
ansivity
grating
spons
unsting
H.
direction
sdcutity
Brortion Vottanee
Explanze (atma
B.
D.
F.
Eiveson Vadlare Explinal jootth.
Pruc Vesteros lalted
Pastinn Villare Esphinal foort..
Pre-Acttiatie Dim
Lanwar
NCB
20
512
2
22
0LS
22
GLB
Pro Act-roottooonn Dom (4)
Pre-Activatoon Dim (4)
Pro- Cro-Arttattonn Den
Figure 4: Pre-activation dimension in visual cortex depends on stimulus condition. (A) Rasters of true (bottom) and predicted (top) neuronal activity, and estimated latents (middle) for 1000 example neurons using NCE with pre-activy rreaactivttion dimension d 3. for drifting grating responses. (B) Fraction variance explained as a function of pre-activation dimension for an example session of grating responses. Plotted values the Reduced Rank Regression (green). Shaded regions are the standard deviation across source/target selections. All experiments contain 500 source and 1000 target neurons, and all scores are computed on a held-out test set. (C,D) Same as A,B for natural image responses. (E,F) Same as A,B capture 95% of
across source/target selections for NCE (purple) and
for spontaneous activity. (G) Estimated pre-activation dimension d sufficient explainable variance in the NCE across the three conditions. Each point is the mean for a single imaging session, error bars indicate the min/max across source/target selections.
(H) Fraction of variance explained by NCE models with varying pre-activation dimension (log scale) across the three conditions. Each line represents a single imaging session.
several special cases that are beyond the case treated in (74) the integral over the weight distribution can be solved, yielding simple exact equations for the latent dynamics.
Eigenvalue decay of random feature kernels In the infinite-width limit, two-layer neural networks with random input weights behave like random randomfeature kernels that depend on the distribution of input weights and the activation function of the neurons in the hidden layer [56] Sec. 9.5]. This functional perspective can be generalized to deep networks [45] (75) and constitute the basis of the Neural Tangent Kernel formalism for studying learning dynamics (76). When the activation function is the ReLU, the decay rate the eigenvalues has been proven to be polynomial
[46
47
58
64
77 78),
when inputs not assumed to be uniformly distributed on the sphere [58 for general results
on dot-produee product kernels, see [57 61 63
79 80). In this work, we propose a simple formula that links the power-law exponent of the eigenvalue decay rate, the power of the rectified-power activation function, and the input dimensionality. This formula, which goes beyond known results [46 47 presented conjecture that we test in simulations.
Latent variable modeling of neuronal activity While most latent variable models of neuronal activity were originally developed for electrophysiological recordings [3, some are tailored for calcium recordings [1]. 83]. These models vary in their mechanistic interpretability: The inferred latents are either abstract variables, for example when the model's mapping from latents to neuronal activity involves a multi-layer neural network [3 2
they can be interpreted as linear factors of the neurons' pre-activations, as 99 TO). With nonlinear dimensionality reduction methods such as CEBRA (84) or Rastermap [85], latent variables also abstract as there is no explicit mapping going

---

from the latents to neuronal activity. In this work (Sec. D. we developed NCE, latent variable model for calcium recordings that models neuronal activity as an interpretable of latent variables. NCE also uses a cross-encoding scheme, which allows it to discard variability not shared across neurons. We demonstrate that NCE capable of identifying a low-dimensional pre-activation space even when the recorded neuronal activity has high linear dimension.
readout
6 Discussion
Dimensionality of neural systems The solvable RNN model we proposed produces eyclic pop- ulation dynamics that is low-dimensional in the space of pre-activations, and nigh-dimensional in the space of post-activations (firing rates). Thus, an RNN can produce trajectories that are simulta- neously low-dimensional and high-dimensional, depending on the variables being considered. In this work, focused on the notion of linear dimensionality and adopted Hilbert space formalism borrowed from kernel methods to characterize the linear dimensionality of firing rate trajectories in the large-netwond limit. Of course, the "intrinsic" dimension of neuronal activity is equal to 1, since the dynamics is periodic; this highlights the important distinction between intrinsic and linear (or "embedding") dimension of neuronal activity (see also
44) The type of 4
high-dimensional activity our model produces is computationally relevant: It can be exploited by a downstream readout neuron to represent arbitrary periodic functions (see also [86 ) or, following
the random readout approach of [87 88). it can be used to represent a Gaussian process prior over periodic functions (see Appendix G
Two definitions of high-dimensional neuronal activity have been studied in neuroscience. Throughout this work, we define high- digennimeniional neuronal activity mean covariance cigenspectrum with a heavy tail that decays strictly faster that 1/n [39]. This definition is well-suited for systems generating activity whose pairwise correlations do not converge to zero in the limit of large network size. In contrast, random chaotic RNNs-solvabll models whose pairwise correlations do converge to zero zero-eerooooee proo-per different form of high-dimensional activity, where the eigenspectrum decays slower than 1/n [89]. reflecting noise that is not shared between neurons. While a comprehensive comparison of these two types of high-dimensional activity is beyond the scope of this work, we mention that the latter is relevant when one wants to study the noisiness of neuronal responses [90].
Pre-activations and subthreshold membrane potentials We developed NCE to disentangle the linear dimension of neuronal activity and of the neurons' pre-activations. We show that even when activity is high-dimensional, it can be well-explained with low -imeenmnniooal pre-activations in the case of grating responses and spontaneous activity. From a biological point of view, how should we interpret the pre-activations inferred by NCE? If one assumes that the link between synaptic integration and neuronal firing is well approximated by simple nonlinear activation function in cortical neurons, as NCE does, one could argue that pre-activations represent estimates of the neurons' synaptic inputs subthreshold membrane potentials. This is an experimental prediction that large-scale voltage imaging of neuronal populations 91 may make testable in the near future.
Limitations Training NCE, which is a non-convex optimization problem, can be challenging when neural data is limited. The duration of imaging experiments is limited to 3.5 h, which yields only a few thousand training examples per session for tens of thousands of neurons: thus, optimization can get stuck in local minima. To facilitate training, we limit to fitting only 1000 highly responsive neurons at a time, and use tools such data augmentation and pretraining as described in Appendix !] The fact that we were not able to accurately predict neuronal responses to natural images with dow-dimensional NCE model does not necessarily exclude the possibility for such a model to exist, and one could possibly find it with a larger training set or a different model class/hyperprammter configuration. Note that this work did not revisit the problem of how to estimate the tail of the shared
covariance eigenspectrum from neuronal recordings, which is discussed in (39 40 92 93
Acknowledgments and Disclosure of Funding
The authors thank Louis Pezon for useful discussions; Michael Krumin, Bex Terry and Charu B. Reddy for experimental support; Kimberly Ren for feedback on the manuscript. We also thank the four anonymous reviewers who have helped improve the manuscript. This work was funded
10

---

by UKRI (Frontier Award EP/X022366/1 to MC), BBSRC (grant BB/W019884/1 to MC), the MC). the Wellcome Trust
National Institutes of Health BRAIN initiative (grant UOINS (Investigator Award 223144/7/2/7. to MC and KDH), and the ERC (101097874 to KDH), MC holds the GlaxoSmithKline Fight for Sight Chair Visual Neuroscience. AH is supported by a studentship from the Gatsby Charitable Foundation (GAT3755) and the Wellcome Trust (219627/2/19/2), VS is supported by a Royal Society Newton International Fellowship from the Swiss National Science Foundation (grant no. 222150).
and a fellowship
References
[1] Mark M Churchland, John P Cunningham, Matthew T Kaufman, Justin Foster, Paul Nuyu- jukian, Stephen Ryu, and Krishna V Shenoy. Neural population dynamics during reaching. Nature, 487(7405):51-56, 2012.
[2] Valerio Mante, David Sussillo, Krishna V Shenoy, and William T Newsome, Context-dependent computation by recurrent dynamics in prefrontal cortex. Nature, 503(7474):78 84, 2013.
[3] Abigail A Russo, Sean R Bittner, Sean M Perkins, Jeffrey S Seely, Brian M London, Antonio H Lara, Andrew Miri, Najja Marshall, Adam Kohn, Thomas M Jessell, et al. Motor cortex embeds muscle-like commands in an untangled population response. Neuron, 97(4):953-966, 2018.
[4] Evan D Remington, Devika Narain, Eghbal A Hosseini, and Mehrdad Jazayeri. Flexible sensorimotor computations through rapid reconfiguration of cortical dynamies. Neuron, 98(5): 1009-1019, 2018.
[5] Krishna V Shenoy, Maneesh Sahani, and Mark Churchland. Cortical control of arm move- ments: dynamical systems perspective. Annual Review of Neuroscience, 36(1):337-359, 2013.
[6] Saurabh Vyas, Matthew D Golub, David Sussillo, and Krishna V Shenoy, Computation through neural population dynamics. Anmual Review of Neuroscience, 43(1):249-211 2020.
[7] Juan A Gallego, Matthew G Perich, Lee E Miller, and Sara A Solla. Neural manifolds for the control of movement. Neuron, 94(5):978-955 2017.
[8] M Yu Byron, John P Cunningham, Gopal Santhanam, Stephen I Ryu, Krishna V Shenoy, and Maneesh Sahani. Gaussian- factor analysis for low-dimensional single-trial analysis of neural population activity. Journal of Neurophysiology, 102(1):614 635, 2009.
[9] Jakob H Macke, Lars Buesing, John P Cunningham, Byron M Yu, Krishna V Shenoy, and Maneesh Sahani. Empirical models of spiking neural populations. In Advances in Neural Information Processing Systems, volume 24, 2011.
[10] Chethan Pandarinath, Daniel O' Oheaa Jasmine Collins, Rafal Jozefowicz, Sergey D. Stavisky, Jonathan C. Kao, Eric M. Trautmann, Matthew T. Kaufman, Stephen I. Ryu, Leigh R. Hochberg, Jaimie M. Henderson, Krishna V. Shenoy, L. F. Abbott, and David Sussillo. Inferring single-trial neural population dynamics using sequential auto-encoders. Nature Methods, 15(10):805-815, 2018. doi: 10. 10.1038/s4159) 100444192--111 00022-01--10000 018-0109-9.
[11] Lea Duncker, Gergo Bohner, Julien Boussard, and Maneesh Sahani. Learning interpretable continuous time models of latent stochastic dynamical systems. In International conference on machine learning. pages 1726-1734, PMLR, 2019.
[12] Joshua Glaser, Matthew Whiteway, John P Cunningham, Liam Paninski, and Scott Linderman. Recurrent switching dynamical systems models for multiple interacting neural populations. In Advances Neural Information Processing Systems, volume 33, pages 14867-14878, 2020.
[13] Ding Zhou and Xue-Xin Wei. Learning identifiable and interpretable latent models of high- dimensional neural activity using pi-vae. Advances Neural Information Processing Systems, volume 33, pages 7234-7247, 2020.
11

---

[14] Timothy D Kim, Thomas Z Luo. Jonathan W Pillow, and Carlos D Brody. Inferring latent dynamics underlying neural population activity via neural differential equations. In International Conference Machine Learning, pages 5551-5561 PMLR, 2021.
[15] Marine Schimel, Ta-Chu Kao, Kristopher T Jensen, and Guillaume Hennequin. iLQR- control-based learning of input-driven dynamics with applications to neural data. In Interna- lional Conference on Learning Representations, 2022.
[16] Mohammad Reza Keshtkaran, Andrew R Sedler, Raeed H Chowdhury, Raghav Tandon, Diya Basrai, Sarah L Nguyen, Hansem Sohn, Mehrdad Jazayeri, Lee E Miller, and Chethan Pandari- nath. A large-scale neural network training framework for generalized estimation of single-trial population dynamics. Nature Methods, 19(12):1522-157; 2022.
[17] Adrian Valente, Jonathan W Pillow, and Srdjan Ostojic. Extracting computational mechanisms from neural data using low-rank RNNs. In Advances in Neural Information Processing Systems, volume 35, pages 24072-24086, 2022.
[18] Shuqi Wang, Valentin Schmutz, Guillaume Bellec, and Wulfram Gerstner. Mesoscopic modeling of hidden spiking neurons. In Advances in Neural Information Processing Systems, volume 35, pages 23566 23579, 2022.
[19] Mikhail Genkin, Krishna V Shenoy, Chandramouli Chandrasekaran, and Tatiana A Engel. The dynamics and geometry of choice in the premotor cortex. Nature, pages 1-9, 2025.
[20] Timothy Doyeon Kim, Thomas Zhihao Luo, Tankut Can, Kamesh Krishnamurthy, Jonathan W. Pillow, and Carlos D Brody. Flow-field inference from neural data using deep recurrent networks. In Forty-second International Conference on Machine Learning, 2025.
[21] Aditi Jha, Diksha Gupta, Carlos Brody, and Jonathan W Pillow. Disentangling the roles of distinct cell classes with cell-type dynamical systems. In Advances in Neural Information Processing Systems, volume 37, pages 33668-33690 2024.
[22] Christopher Langdon and Tatiana A Engel. Latent circuit inference from heterogeneous neural responses during cognitive tasks. Nature Neuroscience, pages 1-11, 2025.
[23] Adam Gosztolai, Robert Peach, Alexis Amaudon, Mauricio Barahona, and Pierre Van- dergheynst. MARBLE: interpretable representations of neural population dynamics using geometric deep learning. Nature Methods, pages 1-9, 2025.
[24] David Sussillo and Omri Barak. Opening the black box: low-dimensional dynamics in high- dimensional recurrent neural networks. Neural Computation, 25(3):626-64)) 2013.
[25] David Sussillo, Mark M Churchland, Matthew T Kaufman, and Krishna V Shenoy. A neural net- work that finds naturalistic solution for the production of muscle activity. Nature Neuroscience,
1 1025)3305) 2015.
[26] Guillaume Hennequin, Tim P Vogels, and Wulfram Gerstner. Optimal control of transient dynamics in balanced networks supports generation of complex movements. Neuron, 82(6): 3394-1406, 2014.
[27] 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 uncertainty. Neuron, 86(4): 1067-10777 2015.
[28] Jonathan A Michaels, Benjamin Dann, and Hansjörg Scherberger. Neural population dynamics during reaching are better explained by a dynamical system than representational tuning. PLoS Computational Biology, 12(11)ee
2016.
[29] Alexander Rivkind and Omri Barak. Local dynamics in trained recurrent neural networks. Physical Review Letters, 11 118(25):288001 2017.
[30] Joao Barbosa, Rémi Proville, Chris C Rodgers, Michael R DeWeese, Srdjan Ostojic, and Yves Boubenec. Early selection of task-relevant features through population gating. Nature Communications, 14(1):6837, 2023.
12

---

[31] Christopher Langdon, Mikhail Genkin, and Tatiana A Engel. A unifying perspective on neural manifolds and circuits for cognition. Nature Reviews Neuroscience, 24(6):363-377, 2023.
[32] Daniel Durstewitz, Georgia Koppe, and Max Ingo Thurm. Reconstructing computational system dynamics from neural data with recurrent neural networks. Nature Reviews Neuroscience, 24 (11):693-710, 2023.
[33] Laura N Driscoll, Krishna Shenoy, and David Sussillo. Flexible multitask computation in recurrent networks utilizes shared dynamical motifs. Nature Neuroscience, 2777):149-1136, 2024.
[34] Harsha Gurnani, Weixuan Liu, and Bingni W Brunton. Feedback control of recurrent dynamics constrains learning timescales during motor adaptation. bioRxiv, pages 2024-02 05, 2024.
[35] Ashok Litwin-Kumar and Brent Doiron. Slow dynamics and high variability in balanced cortical networks with clustered connections. Nature Neuroscience, (((1)):14881155 2012. doi: 10. 10.1038/m..3220.
[36] Francesca Mastrogiuseppe and Srdjan Ostojic. Linking connectivity, dynamics, and computa- tions in low-rank recurrent neural networks. Neuron, 99(3):609-623, 2018.
[37] Brian DePasquale, David Sussillo, LF Abbott, and Mark M Churchland. The centrality of population- factors to network computation is demonstrated by versatile approach for training spiking networks. Neuron, 111(5):631-44) 2023.
[38] Valentin Schmutz, Johanni Brea, and Wulfram Gerstner. Emergent rate-based dynamics in duplicate-free populations of spiking neurons. Physical Review Letters, 134(1):01 8401, 2025,
[39] Carsen Stringer, Marius Pachitariu, Nicholas Steinmetz, Matteo Carandini, and Kenneth D Harris. High-dimensional geometry of population responses in visual cortex. Nature, 571(7765): 361-365, 2019.
[40] Carsen Stringer, Marius Pachitariu, Nicholas Steinmetz, Charu Bai Reddy, Matteo Carandini, and Kenneth D Harris. Spontaneous behaviors drive multidimensional, brainwide activity. Science, 364(64437)eaaxv799) 2019.
[41] Frederic Lanore, Alex Gajic, Harsha Gurnani, Diccon Coyle, and R Angus Silver. Cerebellar granule cell axons support high-dimensional representations. Nature Neuroscience, 148)):1442113: 1150, 2021.
[42] Jason Manley, Sihao Lu, Kevin Barber, Jeffrey Demas, Hyewon Kim, David Meyer, Fran- cisca Martínez Traub, and Alipasha Vaziri. Simultancous, cortex-wide dynamics of up to I million neurons reveal unbounded scaling of dimensionality with neuron number. Neuron, 112 (10):1694-1709, 2024.
[43] Mehrdad Jazayeri and Srdjan Ostojic. Interpreting neural computations by examining intrinsic and embedding dimensionality of neural activity. Current Opinion in Neurobiology, 70:113-120, 2021.
[44] Mark D Humphries. Strong and weak principles of neural dimension reduction. Neurons, Behavior; Data analysis, and Theory, 1-2):1-28, 2021.
[45] Youngmin Cho and Lawrence Saul. Kernel methods for deep learning. In Advances in Neural Information Processing Systems, volume 22, 2009.
[46] Francis Bach. Breaking the curse of dimensionality with convex neural networks. Journal of Machine Learning Research, 18(19):1-53) 2017.
[47] Alberto Bietti and Francis Bach. Deep equals shallow for ReLU networks in kernel regimes. In International Conference on Learning Representations, 2021.
[48] William Skaggs, James Knierim, Hemant Kudrimoti, and Bruce McNaughton. A model of the neural basis of the rat's sense of direction. In Advances in Neural Information Processing Systems, volume 7, 1994.
13

---

[49] Rani Ben- Yishai, R Lev Bar-Or, and Haim Sompolinsky. Theory of orientation tuning in visual cortex. Proceedings of the National Academy of Sciences, 92(9):3844-388) 1995.
[50] Kechen Zhang. Representation of spatial orientation by the intrinsic dynamics of the head- direction cell ensemble: theory. Journal of Neuroscience, 16(6):2112-2126 1996.
[51] Manuel Beiran, Alexis Dubreuil, Adrian Valente, Francesca Mastrogiuseppe, and Srdjan Os- tojie. Shaping dynamics with multiple populations in low-rank recurrent networks. Neural Computation, 33(6):1572-16155 2021.
[52] Louis Pezon, Valentin Schmutz, and Wulfram Gerstner. Linking neural manifolds to circuit structure in recurrent networks. bioRxiv, pages 2024-02, 2024.
[53] Mare Mézard, Giorgio Parisi, and Anthony Zee. Spectra of euclidean random matrices. Nuclear Physics B, 559())669-771,, 1999.
[54] Vladimir Koltchinskii and Evarist Giné. Random matrix approximation of speetra of integral operators. Bernoulli, 6(1):113-16)) 167, 2000. ISSN 1550-7265,1573--77 doi: 10.2307/3318636,
[55] Charles Bordenave. Eigenvalues of euclidean random matrices. Random Structures Algo- rithms, 3355)55 2008. doi: http://doo..rr/11..100///s/sa
[56] Francis Bach. Learning theory from first principles. MIT press, 2024.
[57] Meyer Scetbon and 7aid Harchaoui. A spectral analysis of dot-product kernels. In International Conference on Artificial Intelligence and Statistics, pages
PMLR, 2021.
[58] Yicheng Li, Zixiong Yu, Guhan Chen, and Qian Lin. On the eigenvalue decay rates of class of neural-network related kernel functions defined on general domains. Journal of Machine Learning Research, 25(82):1-44 2024.
[59] Thomas Kühn. Eigenvalues of integral operators with smooth positive definite kernels. Archiv der Mathematik, 49:525-534, 1987.
[60] M Castro and V Menegatto. Eigenvalue decay of positive integral operators on the sphere. Mathematics of Computation, 81 (3(280))203-2231, 2012.
[61] Douglas Azevedo and Valdir Antonio Menegatto. Sharp estimates for eigenvalues of integral operators generated by dot product kernels on the sphere. Journal of Approximation Theory, 177:57-68, 2014.
[62] Thaís Jordão and V Menegatto. Estimates for fourier sums and eigenvalues of integral operators via multipliers on the sphere. Proceedings of the American Mathematical Society, 144(1): 269-283,, 2016.
[63] Blake Bordelon and Cengiz Pehlevan. Population codes enable learning from few examples by shaping inductive bias. Elife, 11:e78606, 2022.
[64] Marjorie Xie, Samuel P Muscinelli, Kameron Decker Harris, and Ashok Kitwin-Kumar. Task- dependent optimal representations for cerebellar learning. Elife, 12:e82914, 2023.
[65] Jeffrey Demas, Jason Manley, Frank Tejera, Kevin Barber, Hyewon Kim, Francisca Martínez Traub, Brandon Chen, and Alipasha Vaziri. High-speed, cortex-wide volumetric recording of neuroactivity at cellular resolution using light beads microscopy. Nature Methods, 18(9):
1103-1111, 2021.
[66] Ali Haydaroglu, Tinya Chang, Andrew Landau, Michael Krumin, Sam Dodgson, Liad J. Baruchin, Maria Cozan, Jingkun Guo, David Meyer, Charu Bai Reddy, Jian Zhong, Na Ji, Sylvia Schröder, Kenneth D. Harris, Alipasha Vaziri, and Matteo Carandini. Suite3D: Volumetric cell detection for two-photon microscopy. bioRxiv, page 2025.03.22 233.26666 26.645628. 2025.
(67] Alan Julian Izenman. Reduced-rant regression for the multivariate linear model. Journal of Multivariate Analysis, 5(2):248-24 2244444,444444 1975.
14

---

[68] Matteo Carandini. What simple and complex cells compute. The Journal of Physiology, 577(Pt 2):463 466, December 2006. ISSN 0022-3751. doi: 10.1113/pphysio02200.11117
[69] Friedrich Schuessler, Alexis Dubreuil, Francesca Mastrogiuseppe, Srdjan Ostojic, and Omri Barak. Dynamics of random recurrent networks with correlated low-rank structure. Physical Review Research, 2(1):013111, 2020.
[70] Alexis Dubreuil, Adrian Valente, Manuel Beiran, Francesca Mastrogiuseppe, and Srdjan Ostojic. The role of population structure in computations through neural dynamics. Nature Neuroscience, 25(6):783-794, 2022.
[71] Yue Wan and Robert Rosenbaum. High-dimensional dynamics in low-dimensional networks. arXiv preprint arXiv:2504.13727, 2025.
[72] Francesca Mastrogiuseppe, Joana Carmona, and Christian Machens. Stochastic activity in low-rank recurrent neural networks. bioRxiv, pages 2025-04, 2025.
[73] Romain Veltz and Olivier Faugeras. Local/global analysis of the stationary solutions of some neural field equations. SIAM Journal on Applied Dynamical Systems, 9(3):55--99) 2010.
[74] Matthijs Pals, Jakob H Macke, and Omri Barak. Trained recurrent neural networks develop phase-locked limit cycles in a working memory task. PLOS Computational Biology, 20(2): c1011852, 2024.
[75] Jachoon Lee, Jascha Sohl-dickstein, Jeffrey Pennington, Roman Novak. Sam Schoenholz, and Yasaman Bahri. Deep neural networks gaussian processes. In International Conference Learning Representations, 2018
[76] Arthur Jacot, Franck Gabriel, and Clement Hongler. Neural tangent kernel: Convergence and generalization in neural networks. In Advances in Neural Information Processing Systems, volume 31, 2018.
[77] Amnon Geifman, Abhay Yadav, Yoni Kasten, Meirav Galun, David Jacobs, and Basri Ronen. On the similarity between the laplace and neural tangent kernels. In Advances in Neural Information Processing Systems, volume 33, pages 1451-1461, 2020.
[78] Alberto Bietti and Julien Mairal. On the inductive bias of neural tangent kernels. In Advances in Neural Information Processing Systems, volume 32, 2019.
[79] Alex Smola, Zoltán Ovári, and Robert C Williamson. Regularization with dot-product kernels. In Advances in Neural Information Processing Systems, volume 13, 2000.
[80] Ha Quang Minh, Partha Niyogi. and Yuan Yao. Mercer's theorem, feature maps, and smoothing. In International Conference on Computational Learning Theory, pages 154-168. Springer,
2006.
[81] Luke Y Prince, Shahab Bakhtiari, Colleen J Gillon, and Blake A Richards. Parallel inference of hierarchical latent dynamics in two-photon calcium imaging of neuronal populations. bioRxiv, pages 2021-03, 2021.
[82] Feng Zhu, Harrison A Grier, Raghav Tandon, Changjia Cai, Anjali Agarwal, Andrea Gio- vannucci, Matthew T Kaufman, and Chethan Pandarinath. A deep learning framework for inference of single-trial neural population dynamics from calcium imaging with subframe temporal resolution. Nature Neuroscience, 25(12):1724-1734 2022.
[83] Tze Hui Koh, William E Bishop, Takashi Kawashima, Brian B Jeon, Ranjani Srinivasan, Yu Mu, Ziqiang Wei, Sandra Kuhlman, Misha B Ahrens, Steven M Chase, et al. Dimensionality reduction calcium-imaged neuronal population activity. Nature Computational Science, 3(1): 71-85. 2023.
[84] Steffen Schneider, Jin Hwa Lee, and Mackenzie Weygandt Mathis. Learnable latent embeddings for joint behavioural and neural analysis. Nature, 617((960):360-300 -6866000-66666 May 2023.
15

---

[85] Carsen Stringer, Lin Zhong, Atika Syeda, Fengtong Du Maria Kesa, and Marius Pachitariu. Rastermap: discovery method for neural population recordings. Nature Neuroscience, 28(1): 201-212, 2025.
[86] Soledad Gonzalo Cogno, Horst A Obenhaus, Ane Lautrup, R Irene Jacobsen, Claudia Clopath, Sebastian O Andersson, Flavio Donato, May-Brit Moser, and Edvard Moser. Minute-scale oscillatory sequences in medial entorhinal cortex. Nature, 625(7994):338 344, 2024.
[87] Radford M Neal. Priors for infinite networks. Bayesian learning for neural networks, pages 29-53, 1996.
[88] Christopher Williams. Computing with infinite networks. In Advances in Neural Information Processing Systems, volume 1996.
[89] David G Clark, LF Abbott, and Ashok Litwin-Kumar. Dimension of activity in random neural networks. Physical Review Letters, (11(11):1 18401, 2023.
[90] Gengshuo John Tian, Ou Zhu, Vinay Shirhatti, Charles M Greenspon, John E Downey, David J Freedman, and Brent Doiron. Neuronal firing rate diversity lowers the dimension of population covariability. bioRxiv, pages 2024-08, 2024.
[91] Jian Zhong, Ryan G. Natan, Qinrong Zhang, Justin S. J. Wong, Christoph Miehl, Krishnashish Bose, Xiaoyu Lu, François St-Pierre, Su Guo, Brent Doiron, Kevin K. Tsia, and Na Ji. Faced 2.0 enables large-scale voltage and calcium imaging in vivo. page 0025003.06.6417888 2025.
[92] Dean Pospisil and Jonathan W. Pillow. Revisiting the high high-dimensioonal geometry of popula- tion responses in visual cortex. bioRxiv, 2024. doi: 10.1101/2024.0..666 16.552222....11.88776.
[93] Marius Pachitariu, Lin Zhong, Alexa Gracias, Amanda Minisi, Crystall Lopez, and Carsen Stringer. A critical initialization for biological neural networks. bioRxiv, pages 2025-01, 2025.
[94] Tsai-Wes Chen, Trevor J. Wardill, Yi Sun, Stefan R. Pulver, Sabine L. Remninger, Amy Baohan, Eric R. Schreiter, Rex A. Kerr, Michael B. Orger, Vivek Jayaraman, Loren L. Looger, Karel Svoboda, and Douglas S. Kim. Ultrasenstivve fluorescent proteins for imaging neuronal activity. Nature, (99(444)))59))25530000 2013.
[95] Geoffrey Hinton, Oriol Vinyals, and Jeff Dean. Distilling the knowledge in neural network. arXiv, (arXiv:1503.0253)) 2015. arXiv:1503.02531.
[96] Diederik P. Kingma and Jimmy Ba. Adam: A method for stochastic optimization. arXiv, (arXiv:1412.6980), 2017. arXiv:1412.6980.
[97] Kaiming He, Xiangyu Zhang, Shaoqing Ren, and Jian Sun. Delving deep into rectifiers: Surpassing human-level performance imagenet classification. In 2015 IEEE International Conference on Computer Vision (ICCV), page 1026-1034. IEEE, 2015.
16

---

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

---

Justification: Our nontrivial mathematical statements are either fully proved or clearly stated as conjectures (e.g. Conjecture in the main text).
Guidelines:
The answer NA means that the paper does not include theoretical results.
All the theorems, formulas, and proofs in the paper should be numbered and cross- referenced.
All assumptions should be clearly stated or referenced the statement of any theorems.
The proofs can either appear in the main paper or the supplemental material, but they appear in the supplemental material, the authors are encouraged to provide short proof sketch to provide intuition.
Inversely, any informal proof provided in the core of the paper should be complemented by formal proofs provided in appendix supplemental material.
Theorems and Lemmas that the proof relies upon should be properly referenced.
4. Experimental result reproducibility
Question: Does the paper fully disclose all the information needed to reproduce the main ex- perimental results of the paper to the extent that affects the main claims and/or conclusions of the paper (regardless of whether the code and data are provided or not)?
Answer: [Yes]
The architecture, training procedure and hyperparameters for the Neural Cross- Encoder are described in detail. While the text descriptions are sufficient to reproduce the results, we also provide the code for in silica experiments. We also share the preprocessed neuronal datasets. The methods for in vivo experiments and the preprocessing of imaging data are described clearly to enable reproducibility. Guidelines:
Justification:
The answer NA means that the paper does not include experiments.
If the paper includes experiments, a No answer to this question will not be perceived well by the reviewers: Making the paper reproducible is important, regardless of whether the code and data are provided or not.
If the contribution is dataset and/or model, the authors should describe the steps taken to make their results reproducible or verifiable.
Depending on the contribution, reproducibility can be accomplished in various ways. For example, if the contribution is novel architecture, describing the architecture fully might suffice, or if the contribution a specific model and empirical evaluation, it may be necessary to either make it possible for others to replicate the model with the same dataset, or provide access to the model. In general. releasing code and data is often one good way accomplish this, but reproducibility instructions for how to replicate the results, access to a hosted model (e.g., in the case of large language model), releasing of a model checkpoint, other means that appropriate to the research performed.
also be provided via detailed
While NeurIPS does not require releasing code, the conference does require all submis- sions to provide some reasonable avenue for reproducibility, which may depend on the nature of the contribution. For example
(a) If the contribution is primarily a new algorithm, the paper should make it clear how reproduce that algorithm. (b) If the contribution is primarily new model architecture, the paper should describe the architecture clearly and fully.
(c) the contribution is a new model (c.g., a large language model), then there should either be way to access this model for reproducing the results or way to reproduce the model (e.g., with an open-source dataset or instructions for how to construct the dataset).
(d) We recognize that reproducibility may be tricky in some cases, in which case authors are welcome describe the particular way they provide for reproducibility. In the case of closed-source models, it may be that access to the model is limited in some way (c.g., to registered users), but it should be possible for other researchers to have some path to reproducing verifying the results.
18

---

5. Open access to data and code
Question: Does the paper provide open access to the data and code, with sufficient instruc- tions faithfully reproduce the main experimental results, described in supplemental material?
Answer: [Yes]
Justification: All code required reproduce the main experimental results anonymized and submitted with the supplementary materials. Preprocessed neural datasets, comprised of deconvolved spontaneous activity and averaged, deconvolved stimulus responses, is also anonymized and submitted with the supplementary materials.
Guidelines:
The answer NA means that paper does not include experiments requiring code.
Please see the NeurlPS code and data submission guidelines (https://nips... cc/
for more details.
While we encourage the release of code and data, we understand that this might not be possible, so "No" is an acceptable answer. Papers cannot be rejected simply for not including code, unless this is central to the contribution (e.g., for new open-source benchmark).
The instructions should contain the exact command and environment needed to run to reproduce the results. See the NeurlPS code and data submission guidelines (https
Sooeeees//ooSSSSmmcmcaa (wwoodeSumiissiooooooooooooo ssionlolicy) for more details.
The authors should provide instructions on data access and preparation, including how to access the raw data, preprocessed data, intermediate data, and generated data, etc.
The authors should provide scripts to reproduce all experimental results for the new proposed method and baselines. only subset of experiments are reproducible, they should state which ones are omitted from the script and why.
At submission time, preserve anonymity, the authors should release anonymized versions (if applicable).
Providing as much information as possible in supplemental material (appended to the paper) recommended, but including URLs to data and code permitted.
6. Experimental setting/details
Question: Does the paper specify all the training and test details (c.g.. data splits, parameters, how they were chosen, type of optimizer, etc.) necessary to understand the results?
Answer: [Yes]
Justification; All relevant details of the training and testing of presented models are described in the main text and appendices, including the dataset generation, hyperparameters, and optimization procedure.
Guidelines:
The answer NA means that the paper does not include experiments.
The experimental setting should be presented in the core of the paper to level of detail that is necessary to appreciate the results and make sense of them.
The full details can be provided either with the code, in appendix, or as supplemental material.
7. Experiment statistical significance
Question: Does the paper report error bars suitably and correctly defined or other appropriate information about the statistical significance of the experiments?
Answer: [Yes]
Justification: On Fig, error bars reported and clearly defined. The factors of variability they represent, including variability from subselection of neurons within a session and variability across session, are elearly stated. On Fig. 2 error bars number of points representing single simulations (> 20 per line) already visually indicates the variance of the simulations.
reported the
Guidelines:
19

---

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

---

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

---

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

---

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

---

Appendix
A Derivation of the latent dynamics Eq. (5) from Eq. (4)
To go from Eq. to Eq. we only have to solve the integral
2x
d0
v(0) K.)
2T
with (d)) = (cos f sin 0)T and v(0) = ( + A), i + A)) (The symbol above denotes the dot product.) Using the polar change of coordinate,
cos
K =
||^|||
sin 8
(14)
we have
| J v(0))(u(0)
2x
x)
2π
db
=
A=0
coss + A)
sin(O A)
à ||ril| sin 0
cos
cos
sin
d0
2m
cos(@ A)) sin (0 + A)
d0
Φb [[|x cos(1 0))
0=
2x
Recalling that @ is the Hcaviside step function O we have
coss( + A) J 0= (8 0 + A) & cos(( - -0))
d0
cos(( A4) + (||x=|| cos(# - 00) 2m
2n
J
0=
sin(O
AA)
04a/2
cos(( + A)
sin + -) d0.
2
0=0-x/2
The integral above can be solved, which yields
J 04-/22
cos(( + A) sin(O + A) d0
J sin(O + A)
-2+0/2
J
sin + 2 + A)
2T
-0=00//2
2x cos(0 (0 + A) 0=0-x/2
π - cos(() + 3/4 + A)
Using the sum and difference trigonometric identities,
si _ + A) sin( - + A) cos 2 + cos( ( + A) sin 8,
(o 10 - 2 + A) = cos(? 2 + A) cas 2 + sin( 3 + A sin 0,
we can write
sin 1 + + A)
sin x + cos
cos + +A) )
cos sin 1
cos( 1 + h A))
sin
=x/||m]
where we revert the polar change of coordinate, Eq. obtain
14 Plugging in J = / and A = T /4, we
J
sin + cos
A) T + +
cos I +
A) V2
k
sin 3n
cos
3x R
sin + A |
cos( 3
sin
3m
||<
B Derivation of pairwise correlations Eq.
66
First, we compute the time-average and the standard deviation of the post-activation of neuron i, ((c0s/)
0, sin 0,
.x)
assuming that the latent variables k rotate at constant speed the unit circle (which is the asymptotic
behavior of the dynamical system Eq. (5). Using again the polar change of coordinate Eq. (14) and recalling that Φ & (the Heaviside step function), we have
((cm0) 0g 0,
x = ((st))
[cos 0: cos 1/ sin ssin0;))
sin 0) )) (ccos(O: - -)).
24

---

The time-average of the post-activation of the neuron i is
(0cos(0, -d)) du
0.47/2
di
Fy
(15)
0=0
2m
xx/2 0=01-
2T
2
Similarly, the variance the post-activation of neuron i is
(ccos(0; - 0))2 do
01tn/2 d0
of =
r?
(16)
-=
2x
0==6- - 2n
4
4
Hence, the correlation Ciy between any pair of neurons i and j js
du
Ci = - (ccos(0, - ())s(( - +))
(17)
2π
Developing the product in the integrand above, we see that only remains to compute the integral
d0
(ccos(0, - ()))((c ())
2m
Using the change of variable 8 0: + =
0, we have
2
2x
Ocoos(0, - 0))0((osO0 - 0)) d0
2
dē
(ccos(1 0))0(( + +0 - 0; - ())
0=0
2n
2m
(ccos 0 0 0: 2)) dē 2T
1=0
Without loss of generality, we assume that 0, then distinguish two cases: If 02 02 n
02 2 0 (otherwise, permute the indices j and i). We
(cos( @ +0; 0. }) dū h 1/
dØ
(nT 0j 02).
2T
0=0;-00
2T
2x
if 0 - 0; > n,
du
0,-0i-5 dd
1
(cos(( + 0, - 0:
7))
(0 - -)).
2n
J0=0
2TT
2T
However, by the definition of the absolute angle difference 10; - 031 = cos -1 (coos (0: 02,)). we know that
π - 0j 0;, if 0 0, S n
π - 10 - -0,1 =
=00 0; if 0, 0, >
Hence,
2m
do
Ocoos (0: 0) 0)))((c 0))
(Λ 10, 021),
(18)
2
d=0
and verify that the equality above is invariant upon permutation of the indices and Using Eqs. (15). (16), and (18) in Eq. (17). we obtain
Oz 2 - 0,1) 1.
C Simulation of an RNN with finite number of neurons
To verify that the theory for infinite-size networks presented in Sec. gives good approximation of the activity RNNs with a large but finite number of neurons, we simulated the dynamies of RNN with N 1000 neurons. The simulation shows that the latent variables & are attracted to limit cycle (Fig. 5 and covariance the network's 's activity has power-law eigenspectrum (Fig. 5 as predicted by the theory for the infinite-size network.
25

---

A
B
norn. eigenvalue
10
10
10
10
10
K1
rank
Figure 5: Simulation of an RNN with N = 1000 neurons. (A) Same as Fig. [D, but with the addition of the trajectory (red curve) of the latent variables k obtained from a simulation of finite-size RNN. (B) Same Fig. D but with the addition of the PCA eigenspectrum (red dashed curve) obtained from the simulated activity the finite-size RNN (only the first 400 eigenvalues are shown). The theory derived for infinite-size networks gives good approximation of the activity of an RNN with N 1000 neurons.
D Solvable model with more general lateral connectivity
To illustrate that the analytical results obtained in Sec. 2 can be generalized models with lateral connectivity that are more general than the rank-2, cosine connectivity prescribed by Eq. (2), we present here solvable model where the weight matrix W is of rank (inspired by model studied in [52]. Let us consider an RNN with N neurons and with weight matrix given by
W e J -0---- 03- A) + 8 c2) 0,)- -)],
(19)
where the parameter 2 O fixes the importance of the higher frequency component, cos(2(0:
0z) A). of the connectivity. As in Sec.
2 the activation function á is assumed to be the Heaviside
step function e, Note that when 3 = 0, we recover the rank-2 RNN of Sec. 2 Following the same steps as in Sec. 2.1 1. and defining == U x, the 4-dimensional vector of latent variables, obtain
that in the limit + 00, the dynamics of the latent variables ks follows
d0
K
v(0)
K)
(20)
6=
2T
where
u(f) = (cos(0), sin cos( (28), sissinnn))) (20))
v(0) J (cos(( A), sin + A), B (2s A), sin(20 + A))
If we assume that at time t 0, K3 0) Kg (0) i..., the third and forth latent variables are null, we can use a polar change of coordinate similar to that used in Appendix A and write
u(f) K <x||
cos(1)
0)-
Then, recalling that = O, we can solve the integral in Eq. 20 and, choosing J = TV2 and
A // as in Sec. 2 we find that
cos(O + A)
2x
J
00
+ A
θ [[]==| 8))
d0
1
K2
9=0
cosl 20 A)
2
(20 A)
0
From the equality above, we deduce that xn (f), (2(() (1) follow the same dynamics as the rank-2 network of Sec. 2 and (t( sa(t) 0, for all 0. i.c., the third and forth latent variables
26

---

1.0
0.5
1
K2
0.0
Kig
K4g
0..5
-1.0
10 20
30
40
50
Figure 6: Simulation of the latent dynamics of an RNN with connectivity given by Eq. (9). In this simulation, N = 10a, 3 = 0.25, and all latent variables are initialized at at time 0. The simulation shows that K2 and K2 converge limit cycle, whereas K3 and remain close 0.
remain null at all times. Although we do not perform here a formal stability analysis the limit cycle in the plane (K1, K2, 0,0), verified numerically that the latent variables of RNN with N 10* neurons follow a limit cycle confined to the plane (Mi, f2z O, 0), when 3 is not too large (Fig.
Since the generalized model we have proposed can produce a limit cycle in the space or pre-activations that identical to that the rank-2 model of Sec. 2 we know that, when this limit cycle, the post-activations of the neurons produce power-law eigenspectrum identical to that of the rank-2 model. Hence, we have shown that the strong assumption of a rank-2, cosine weight matrix, Eq. (2). not strictly necessary for an RNN to produce 2-dimensional latent dynamics and high-dimensional activity in the space of ppst-activasions.
E Random low-rank RNNs with tractable latent dynamics
The model analyzed Sec. 2 is not the only model for which a simple expression for the latent dynamics can be obtained. To illustrate this fact, we re-use the general definition of RNN dynamics, Eq. (T), but consider here the class of low-rank RNNs with rank-d weights matrices W given by
W UMUT
(21)
where U an x d random matrix with i.i.d. standard normal entries and M E Rdxd is an arbitrary "overlap" matrix [5I]. As explained in [51], thanks to the fact that the "patterns" U Gaussian, the dynamics the d-dimensional latent variable vector & == U can be described, in the large-network limit N xo, by the dynamical system
By~A((, (Yo)) M
(22)
|||]||
In the present work, we go one step further and ask whether the expectation Ey (Y [|+[])] be solved, in order yield simpler expressions for the latent dynamics. The following proposition shows that the answer is positive when the activation function ø is either (i) step function, (ii) ReLU function, or (iii) Gaussian cumulative distribution function
Proposition 1. (i) For any bias b E R. if d is the step function = -(x + b), Eq.
22
reduces to
k = K V2
exp
122
M
K
|<]
(ii) For any bias E RR if s is the Rel.U function @ max(0. xo + b). Eq.
reduces to 22
1
k = ー
+rf (vain))) Mr.
(23)
2
27

---

(ii)) For any mean E R and variance x² 2
0, if Φ is the Gaussian c.d.f. 6(x)
<2z0 Ji=- exp (- [t-n)²
dl, Eq. 22 reduces to
20
e 122/0220++1
M K V2T(02+ + ||xl)3) |||]
A
(24)
The proposition shows the nontrivial effect of the activation function has on the latent dynamics. For example, let take again d = 2. When & is the step function with bias b 0, and if
-1
M = 22x
we recover the latent dynamics Eq. ( which generates stable limit cycle.
In contrast, when é is the ReLU function with bias O, Eq. becomes the linear system = K + MR, which can generate rotational dynamics but cannot generate a stable limit cycle. Note that Eq. (24) generalizes a previous result obtained for the erf activation function, used as an approximation of tanh [74].
Proof. The prove these statements, suffices to solve the integral,
Ey~N((.1) (Y = y=-oc
( e dy.
2x
for the three different functions d.
For (i).
1 f y | + b)
/2
2 dy
1
2m
dy
2TT
exp
x|
6
y=
For (ii),
2 dy
100
y max w/1/|| b)
s(y||< +)
dyy
x=-x
2xx
2
W=-
and we integrate by parts,
6)
4 dy
22T
INT
-wlllacl + 6)
222
+||x||
+
2 dy
/2m
2
TAT
TAT
=0
1 (1 erf (v ||||||-
2
For (iii), we first integrate by parts,
[ x=- y {
1
22 dy
c o dy
d'(y|[=/]]
=
22x
22
==-
exp (- [w]leel]
4 dy.
2TTO
20
2m
y=
28

---

Then, we complete the square,
1
(11|=|1 2)²
$ dy
V2Ta
exp
20
V2T | y22
2p0 11y
112 + 02y²
exp
dy
V
V2TT
20²
y=-00 00
80
02 + ||>|||
11 11||es1
1 10/11
112
2m
VTTa
cxp
2x²
σ² + ||=||3
+
202 + ||x||33)
20
dy
y=-0
1
0²
1 exp
0² + ||=11
11 |1||
2
y/
dy
V /0 + 1|+112)
V2To
20²
-² + ||||
10=00
112
x exp
2( + ||0||3)
F From the eigenvalues of the activity matrix A (N,2) to the eigenvalues of the random feature kernel operator Eq. (IO)
First, we recall that the eigenvalues \ ( 2 X1 ((NT) 2 of the "covariance matrix",
C(N.T) = A (M.T) ( (M
T
are related to the eigenvalues 1/1
(N,T)
2 1/2
(N.T)
2 of the "Gram matrix",
G (N.T)
Indeed, using the singular value decomposition of A (N,T), one can casily verify that
(N.T)
(N,T)
Th
Vn E min,
N
T
By the law of large numbers, in the limit N - oc, the entries of the Gram matrix (NNTT) converge to an expectation corresponding to the random feature kernel Eq. (11): For any S 8, T,
(N.T)
0.5.
BE-- 22). (p.)p. K Ky (Zs, 2t),
N
N-+00
where u; above denotes the i-th row of the N d random matrix U. Hence, N - co, the eigenvalues (N.T) of the Gram matrix (NNT)) converge the eigenvalues 1/(²² of the kernel matrix
( = (KnealZs, IA))
Furthermore, we have that
A
(N.T)
(T)
Mp
lim N00
T
Finally, by a theorem from Koltchinskii and Giné [54 Theorem 3.1], the scaled eigenvalues y/ (T) /T of the of the kernel matrix (T) converge (in &a sense) to the eigenvalues X, of the integral operator Eq. (10). In conclusion, have that, for any given n,
X (N.T)
lim
lim
An.
T-300
N+000
N
29

---

G Gaussian processes from RNNs
As in Appendix B
we start from infinite-size network where the latent variables rotate at constant speed on the unit circle. Using again the change of coordinate Eq. [14.. the post-activation (firing rate) of neuron with location 0g is
(cos(B, d)), where i is the angle of the latent variable k. Also, from Eq. ( one can easily deduce that the norm of the velocity vector is always /cos(102 + sin (i)( = 1, which implies that the cycle time period of the system is 2m.
Taking M independently and uniformly sampled neurons on the circle, we consider the random projection readout
M
(M)
VM Σ w. (0cos(0, -)).
where the readout weights W, u,..... ÜMM are i.i.d. standard normal variables that are also indepen- dent from the random angles 0g. Since the readout (M) is weighted sum of 2 per-oddiodi functions, g(M) is itself -rrppoiooic function.
As the number of randomly sampled neurons M tends to infinity, the random readout g converges in law to a Gaussian process on the circle S!. To see this, we first compute the mean and covariance of the 27-periodic function ()M), as M - oc. For the mean, we have
M
1
E
प (cos(() 0))
E L E
(Ocos((,
0)}] 0,
for any M. For the covariance, we have, we have, for any pair of angles 1 and d',
M
E
@,ooss - -0)
1
Σ
87))
M
M
1 M
Σ
Occos(O: 0))o(cos((1 - - )))),
M
using that El E(wū = 1 and F( = when + j. By the law of large numbers,
1 M
M Σ O(cos(0, (aonn(;;,))e0ii0 (cos(() )')) M+
B - 0)0(0 8')]
and using the same line of arguments we used to derive Eq. [0-2(((0.==) (0(0- 0)0(0 0'I Then, by the multivariate central limit theorem, have that, for any finite sequence of angles 01, 02..... 0m
(18) in Appendix B we obtain
1 0' (1), 8').
(M)) (M) (B2).
((M) (0m)))
c X ~ N(0, E),
(25)
M-000
where
π 101 02l
π
01 1
Σ = 02
01
π
π
02
1
π 10m 01 10 m 021
By the Kolmogorov extension theorem, there exists a Gaussian process on S g ~ GP ko). such that X can be replaced by (g(01), (12))......) ... g ) in Eq. (25 for any finite sequence of angles 0,, 0z..................... ... 0m. Hence, we have shown that the random projection readout (MM) defines random 2T perpooii function that has the same law as a Gaussian process as M + oc.
In summary, the above arguments show that our solvable RNN model can be used to draw samples from periodic Gaussian process. In other words, it can represent a Gaussian process prior over periodic functions. This indicates that the Gaussian process framework originally developed for infinite-width feedforward neural networks with recurrent dynamics.
also applies, in certain cases, to large networks
30

---

H In vivo experimental methods
Neuronal Recordings Transgenic mice expressing the GCamP6s calcium indicator [94] in all excitatory cortical neurons (CamKIl x Ail62) were implanted with a mm imaging window over visual cortex. Mice were headfixed under a custom Light Beads Microscope [65], modified with an optical stabilization module. A volumetric field of view of 2 x x 0.3 mm over primary and higher visual areas was imaged with 2-3 pm lateral and 20 pm axial voxel spacing, with volume rates between 4.1.4.5 Hz. During recordings, mice were not anesthetized and were free to run on a low-resistance treadmill.
Visual Stimulation Visual stimuli were delivered via three screens that covered approximately 270" horizontally and 70" vertically of the visual field of the mouse. During drifting grating experiments, full-field drifting gratings with a temporal frequency of 2.0 Hz and spatial frequency of 0.04 cycles/degree were presented for 2s, with an inter-stimulus interval (ISI) drawn from a uniform distribution between 2-3 s. sequence of 720 stimuli was presented with random directions drawn from uniform distribution of integer angles ranging from 09 to 359". The same set of 720 stimuli was repeated in a new random order in the second half of the experiment. For natural image experiments, 1866 unique stimuli were presented, each composed of a mosaic of three distinct natural images. The images were chosen from the same set presented in [39]. Each stimulus was presented for 0.8s with random ISI between 0.9-1.3 Stimuli were presented in two shuffled sequences such that each stimulus was presented twice. For sessions of spontaneous activity, no stimuli were displayed and the screens were grey. Mice were awake, free to run (as in the other sessions) for 10 to 15 minutes of imaging.
Data preprocessing To detect cells and extract neuronal activity from volumetric fluorescence movies, we used Suite3D (66), an analysis pipeline for notion-correecton, cell detection and signal extraction from volumetric two-photon data. After from 19,223 + 2,948 (mean +: std) neurons per experiment, across 10 experiments (4 spontaneous activity, 3 natural images, 3 drifting gratings, across three mice). The extracted cell fluorescence was deconvolved after neuropil substraction. Deconvolved traces were re-sampled 5 Hz and aligned to stimulus onset events. For stimulus response sessions, we projected out the activity along the top 30 principal components (PCs) of spontancous activity (computed on a spontancous session recorded consecutively from the same mouse), similar to what was done in (39]. To compute stimulus responses, neuronal activity within a response window was averaged (0.5s for images, 2.0 for gratings) to produce a single value per neuron per stimulus presentation. Deconvolved spontaneous activity was resampled using linear interpolation at 5 Hz.
curation, we analyzed activity
The recordings covered part of primary visual cortex, higher visual areas AM, PM, LM, and RL, with occasional inclusion of some somatosensory areas. We confirmed the recorded area via retinotopic mapping using sparse noise stimuli. To determine stimulus responsiveness, we used the signal-related variance metric to determine what fraction of a neuron's 's variance could be explained by stimuli (described in Stringer et al. Nature 2019). In typical recording with 22,000 neurons, the top 1000 cells had 49% stimulus-relates variance, the next 1000 had 22%, 14%, 10% and so on, where the 10,000th cell had 3%. The responsive neurons were distributed across primary and higher visual cortices, and the selection of source and target neurons was agnostic to location.
I Neural Cross-Encoder
L1 Datasets
For all three data conditions (gratings, natural images and spontaneous activity), we compared the performance of the (nonlinear) NCE model to (linear) Reduced Rank Regression on the same datasets. For results in Fig. 4 we used datasets with 500 source neurons and 1,000 target neurons. For stimulus- driven conditions (gratings and images), the source neurons were chosen to be the 500 neurons with the most reliable stimulus-evoked responses, quantified by the fraction of stimulus-related variance defined in Ref. (39]. 1,000 target neurons were chosen from the 2,000 non-source neurons with the highest stimulus- related variance. For the spontaneous condition, the target and source neurons were selected at random from the full population. We repeated the model fit for each condition and session
31

---

five times with random initializations and neuron selections. In all cases, source and target neuron sets were non-overlapping.
To reduce the impact of non-stimulus related activity in stimulus-driven conditions, we used the fact that each stimulus was presented multiple times, and trained the model to predict the activity of the target neurons on one repeat of a stimulus from the activity of the source neurons on a different repeat, using all possible permutations. For example, if the same stimulus was shown at times 1ı 12, the network was trained to minimize the sum of the squared Euclidean distances [lb:, - 1 2 + l - b,, 112, where b is the recorded activity and b is the NCE reconstruction.
For the spontaneous condition, capture the temporal structure in the latent dynamics, used multiple timepoints of the source neurons to estimate the latent variable at a single timepoint. To predict z,, we took window of size L. timepoints and concatenated the source activity within this time window, increasing the input dimensionality of the encoder: 2, = E(a ...... ay.... ........444 -11))-
The number and size of the hidden layers were unchanged. For the experiments presented, the time window was 1.2s s (L= timepoints).
In all three conditions, the paired source-target datasets were split into vain-validation-tee sets with 50%-20%%30%% split. For stimulus-driven datasets, the split was done using stimulus identity to ensure that all permuted pairs of repeats of a single stimulus were in the same set after splitting. For spontaneous datasets, instead of assigning each timepoint to one of three sets independently, separated the time series into 10s chunks with 2 buffers between each chunk, and assigned each chunk one of three sets. This procedure prevents the training from contaminating the test set through the slow temporal autocorrelation of calcium signals.
1.2 Model Fitting
Linear models were fit on the training set using the closed-form solution of Reduced Rank Regression with ridge penalty. The ridge penalty (which varied between 10-1-10% giving the best performance in the validation set was selected.
For the NCE, the encoder is a 3-layer feedforward network with ReLU activation functions. The hidden layers contain (500, 250, 100) units. The parameters of the encoder & and the parameters of the decoder, i.e., the readout weights and biases U E Bxd and e E RB, the activation parameter and the baseline firing rates ro BBB were learned by minimizing the mean-squared error (MSE) where b, = &p.e p (U)( + r (b, is the predicted activity of the target neurons, a, is the recorded activity of the source neurons, and & represents the encoder). NCE models were trained in two phases: pretraining and fine-tuning. In the pretraining phase, we used a form of knowledge distillation [95]. Instead of predicting the activity of target neurons directly. the NCE was trained on the predicted activity of target neurons from the best-fitting linear model, one that potentially had a higher pre-activation dimension, d, than the NCE. In the fine-tuning phase, the NCE was trained to predict the true neuronal activity.
loss T El_ |b; b,l ?
In both phases of training, we augmented the training set by adding noise to the activity of source neurons. First, we added Poisson-like shot noise present in two-photon recordings by generating zero-mean Gaussian noise with variance equal to the activity of each neuron at a given timepoint, scaled by a factor of 0.05. Next, we added independent Gaussian noise with mean of 0 and variance of 0.05. Finally, we incorporated multiplicative noise, multiplying the activity of each neuron at each timepoint with an independent Gaussian random variable with a mean of and variance of 0.05. The original source dataset was concatenated with three noise-augmented datasets to produce the final training set.
In each phase, parameters were optimized through stochastic gradient descent using Adam 96) and a batch size of 2,048, with a learning rate of 10 momentum parameters Bn = 0.85, 82 = 0.95, and an L2 penalty of 10-5, Model parameters were randomly initialized with Kaiming initialization 197]. Within each phase, we iterated between optimizing the readout parameters, encoder parameters, and all parameters. At each iteration, the model with the best validation performance, evaluated using the coefficient of determination, was selected. Once models were fitted, we report the fraction of explained variance on the test set. To normalize this value to only account for the "explainable" variance, we computed the coefficient of determination of each model, and normalized by the maximum of the best model across all pre-activation dimensionalities for the given set of neurons.
32

---

1.3 Hyperparameter Selection
The architecture of the encoder, the numbers of source and target neurons, well as the training parameters were selected to balance computational efficiency with performance. The encoder is feedforward network with parametrized by the number of layers (depth) and the number of units in the first layer (width), with the width each successive layer reduced by a factor of two. We conducted grid searches example sessions across range of values of and found that the resulting performance is only mildly sensitive to the encoder of the architecture (see Tables 12 for representative examples).
Table Normalized R² for an example grating session with d
Layer size # Layers
3
1
1000
0.97
0,97
0.96
500
0.97
0.97
0.94
250
0.95
0.96
0.94
Table 2: Normalized R² for an example natural image session with d
Layer size # Layers
3
1
1000
0.54
0.55
0.52
500
0.55
0.56
0.53
250
0.56
0.54
0.51
Table Normalized R² for an example natural image session with 32
Layer size # Layers
3
1
1000
0.91
0.89
0.89
500
0.91
0.92
0.89
250
0.92
0.91
0.91
Table Normalized R² for an example natural image session with 64
Layer size # Layers
3
1
1000
0.97
0,98
0.96
500
0.95
0.92
0.98
250
0.98
0.99
1.00
Due to computational constraints, we subset the full population of neurons into 1000 source and 500 target neurons. Varying these numbers up to 2000 neurons does not substantially affect the results shown in Figure 4 Tables 5 repoots the fraction of explainable variance explained (normalized R²) for representative example session of grating responses, fit with varying subset sizes. Increasing the number of cells beyond 2000 somewhat reduces performance and we suggest two potential explanations. First, due to animal health concerns our experiments are of limited duration, hence typical session produces only 1000 training examples. This may not enough to train a relatively large encoder network with 1000 input neurons. Second, only fraction of neurons recorded have reliable stimulus responses, and including more cells introduces a large amount of non-stimulus- related variance.
33

---

Table 5: Normalized R2t with an NCE model for an example grating session with d
Target # Source
100
500
1000
2000
100
0.98
0.96
0.96
0.97
500
0.98
0.97
0.94
0.93
1000
0.98
0.97
0.94
0.91
2000
0.99
0.97
0.97
0.92
Table 6: Normalized R with an RRR model for example grating session with d =
Target Source
100
500
1000
2000
100
0.92
0.88
0.86
0.86
500
0,91
0.85
0.84
0.84
1000
0.89
0.84
0.83
0.83
2000
0.90
0.85
0.87
0.84
For the remaining hyperparameters, we explored the hyperparameter space with the Optuna optimizer (Akiba al. KDD 2019). We found that the performance of the NCE is not sensitive modest variations in learning rate, batch size, and momentum.
1.4 Latent dynamics of spontaneous activity
To determine whether the extracted latent variables related to the behavioral state of the animal, compared them to the running speed of the animal. For each session of spontaneous activity, we trained an NCE model with pre-activation dimension of 4, described above. After training, took the activity of the source neurons across the entire session (removing the train/validationntest splits), and computed the inferred latents.
The running speed of the mouse is highly correlated with inferred latent variables. We computed correlation coefficient between each latent and the running speed, and found high positive or negative correlations for most variables. This supports the hypothesis that the latents inferred by the NCE are related to the behavioral state of the animal.
cuaning -11 MAL. Mo w LAMNNU
(0.44)
(0.00
camii
1068)
1036)
10.18)
Dunning
running
NUA
wWW.
(40.0)
1-7721
(0.08)
10:031
(0.66)
(024)
10.601
Time (nlintes)
Cninutes)
Figure 7: Extracted latents correlate with running speed. For four example sessions of sponta- neous activity, the dynamics the extracted latents 2; (purple) are correlated with the running speed (black). Correlation coefficient between each latent and running speed is reported in parentheses on the y-axis. All curves are smoothed with 1s Gaussian kernel for visualization only.
34

---

1,5 Compute Resources
Volumetric calcium movies were processed using custom workstation running Ubuntu 20.04.6 LTS with an Intel()) Xcon(R) w9-3475X CPU (36 cores, 2.20 GHz), NVIDIA RTX A4500 GPU (20GB VRAM), and 512 GB DDR5 RAM. Preprocessing with Suite3D took <20 hours per session.
All NCE experiments were performed on workstation running Windows 10 Pro, with an Intel (R) CT))(TT 17-11700k CPU (8 cores, 3.6 GHz), NVIDIA RTX 3060 GPU (12GB VRAM), and 128 GB DDR4 RAM. Total runtime of all NCE experiments was <96 hours. All preliminary experiments that are not reported in the paper were conducted on this workstation, with similar runtimes.
1.6 Software and Licenses
This work makes use of several publicly available software resources including: Python (PSF License), Matplotlib (PSF License), NumPy (BSD License), SciPy (BSD License), CuPy (MIT License), Py PyTorch (BSD-3 License), Suite3D (AGPL-3 License).
1.7 Code and Data Availability
All data and code used to train the NCE and generate the results, as well as code used for simula- tions in Fig- 2 are shared via Figshare: httss//ffgharr..com///aa222222 2e2.aa772222222997746443633 1829b7ec406 and Github:
35


Reviews:
[{"review_id": "hJyJp2yREp", "reviewer": "Reviewer_jahV", "summary": "The authors present a theoretical framework to resolve the apparent contradiction between experimentally observed high-dimensional neuronal activity and commonly hypothesized low-dimensional neuronal latents. Following their theoretical model of neuronal computations as a linear-nonlinear system, the authors propose a latent variable identification method, namely Neural Cross-Encoder, to identify the low-dimensional latents underlying high-dimensional neuronal activity. They demonstrate that this method reveals low-dimensional latent variable space underlying spontaneous neuronal responses from the visual cortex and responses to drifting gratings but not responses to naturalistic stimuli. Taken together, this work advances our understanding of neuronal computations in response to stimuli as well as build on methods to identify latent variables underlying these computations.", "questions": "1. While I enjoyed reading Section 2 (and the appendix), I am not quite sure what it adds to the paper. Section 3 and further sections do not require an RNN model assumption, and can be understood even without the results in Section 2. Could the authors clarify why they feel Section 2 is important in the context of the paper?\n2. The choice of the size of source and target neuronal populations is unclear. Do the results change if these choices are changed? Also, can you add an experiment to demonstrate the utility of NCE in simulated settings?\n3. Can you compare NCE to CEBRA or Rastermap or other latent variable methods? Given this is one of the core contributions of this work, it is imperative to show that NCE is indeed a performant latent variable method, at least comparable to existing methods.\n4. In the theory section, the latent variables were assumed to be drawn i.i.d. from a uniform distribution. Can you comment on this assumption, how reasonable is it and whether this can be relaxed? Does this affect NCE performance?", "limitations": "Yes.", "rating": 5, "confidence": 5, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 4, "clarity": 3, "significance": 3, "originality": 4, "strengths_and_weaknesses": "Strengths:\n1. The paper is well-written, reasonably easy to follow and deals with a very pertinent question in the field. Resolving this apparent contradiction in theoretical/modelling approaches and experimental observations will have a significant impact on unifying and advancing the field further. \n2. The theoretical connections between the covariance matrix eigenspectrum and the eigenspectrum of the integral operator is very interesting and well-described. Section 3, overall, is a nice read and sets up the reader for the NCE method. \n3. The NCE method description is very well written and shows promising results, especially on the neuronal recordings. The semantic relevance of some of the identified latent variables is also a very interesting result, and potentially highlights the utility of the method. \n\nWeaknesses:\n1. I think a core contribution of the paper is the NCE method, and it falls short on highlighting the correctness and utility of the method. Firstly, the method has certain hyperparameters (e.g. size of the source and target neuronal populations) that seemed to have been chosen without proper justification. Moreover, the authors do not comment (in the main text) on the robustness of the method to the choice of these hyperparameter values. Finally, the linear model baseline seems like a weak comparison and underscores the utility of the method. Ideally, the authors would also show the method's abilities to uncover latent variables in some simulation setting before demonstrating their findings in recorded neuronal data.\n2. While the theoretical contributions are very interesting, they rely on the assumption of infinite neurons. Real systems are finite dimensional, and it would be interesting to see how these correspondences hold in (high-dimensional) finite dimensional systems. At the least, I think the authors should add a note or comment on these assumptions. \n3. Conjecture 1 is interesting (and a very cool result), but it warrants further justification. In the current version, the reader lacks context as to how the authors came up with this conjecture. The authors could probably add a few statements in this section to improve the readability of this section.", "was_revised": true, "final_justification": "The authors' responses clarified my concerns around the method details and hyperparameter chocies. While I believe more work is required to concretely establish the relevance (and novel contribution) of Section 2 as well as the utility of NCE compared to other methods in literature, I believe this is a technically solid submission and makes a notable contribution towards resolving the apparent contradiction between empirically observed high-dimensional neural representations and the commonly hypothesized low-dimensional latent factors. Moreover, I think NCE is a neat contribution and will have a wide impact on the field of neural data analysis. Therefore, I recommend accepting this submission.", "rebuttal": "The reviewer appreciate the relevance of our results for systems neuroscience and highlight the quality of our text.\n### Weakness 1 \n#### (i) Choice of NCE hyperparameters\nThe reviewer is correct, the current manuscript does not provide sufficient justification of the selection of 500 source and 1000 target neurons and other parameters. Varying the number of source and target neurons up to 2000 does not substantially affect the results shown in Figure 3. We selected 500 and 1000 as intermediate values for computational reasons, since we fit 13 models with 5 shuffles each in 10 sessions for a total of 650 models (we have added 3 more experimental sessions since the original submission, confirming and strengthening results). To illustrate this stability, in the table below, we show the fraction of explainable variance explained for a model with $d=4$ on drifting grating responses, for NCE and RRR models, across combinations of source and target cells.\n| N_target \\ N_source | 100               | 500            | 1000           | 2000           |\n|---------------------|------------------|----------------|----------------|----------------|\n| 100                 | NCE: 0.98 / RRR: 0.92 | 0.96 / 0.88   | 0.96 / 0.86   | 0.97 / 0.86   |\n| 500                 | 0.98 / 0.91       | 0.97 / 0.85     | 0.94 / 0.84     | 0.93 / 0.84     |\n| 1000                | 0.98 / 0.89       | 0.97 / 0.84     | 0.94 / 0.83     | 0.91 / 0.83     |\n| 2000                | 0.99 / 0.90       | 0.97 / 0.85     | 0.97 / 0.87     | 0.92 / 0.84     |\n\nWe note that increasing the number of source cells to 2000 and beyond somewhat reduces performance and we suggest two potential explanations. First, due to experimental constaints,  a typical session produces ~1000 training examples. This may not be enough to train a large encoder network with >1000 input neurons. Second, only a fraction of neurons recorded have reliable stimulus responses, and including more cells introduces a large amount of non-stimulus-related variance.\n\nOf the remaining hyperparameters, we have found that the performance is not sensitive to variations in learning rate, batch size, and momentum by exploring the hyperparameter space with the Optuna optimizer (Akiba et al. KDD 2019). The results are mildly sensitive to the depth and width of the encoder network and weight decay; these hyperparameters were seleced via grid search on an example dataset.\n\nTo demonstrate the extent of the dependence on encoder architecture, we present the fraction of explainable variance explained for two example datasets as we vary the number of layers, and the size of the first layer. For multi-layer encoders, each successive layer is smaller by a factor of two. \n\n**Drifting gratings (d = 4):**\n\n| Layer size \\ # Layers | 3     | 2     | 1     |\n|-----------------------|-------|-------|-------|\n| 1000                  | 0.97  | 0.97  | 0.96  |\n| 500                   | 0.97  | 0.97  | 0.94  |\n| 250                   | 0.95  | 0.96  | 0.94  |\n\n**Natural images (d = 32):**\n\n| Layer size \\ # Layers | 3     | 2     | 1     |\n|-----------------------|-------|-------|-------|\n| 1000                  | 0.91  | 0.89  | 0.89  |\n| 500                   | 0.91  | 0.92  | 0.89  |\n| 250                   | 0.92  | 0.91  | 0.91  |\n\nWe will repeat the measurements shown in the tables above across multiple sessions and different values of $d$, and include them as supplementary figures in the appendix. We will also include in the appendix the details of hyperparameter selection, described above. \n#### (ii) Linear model baseline \nThe linear peer-prediction method we use as a baseline is closely related to linear models used in several influential works, e.g., for the estimation of dimensionality of communication subspaces (Semedo et al. Cell 2019) or spontaneous activity (Stringer et al. Science 2019). For a discussion of alternative latent variable models as a baseline, see our response to Question 3 below.\n#### (iii) Testing NCE on simulated data\nWe performed a test of the NCE method on simulated data generated by the solvable model presented in Section 2. This new analysis shows that the NCE can perfectly recover the two latent variables and it gives an accurate estimate of the pre-activation dimensionality (which is 2). In comparison, the linear model overestimates the dimensionality, estimating it at 4. These new results will be presented in an additional figure placed at the beginning of Section 4, making a natural link between the theoretical results (Sections 2 and 3) and the data analysis part of our work (Section 4).\n### Weakness 2 \nOur theory, which assumes an infinite population, gives a good prediction of the dynamics and covariance eigenspectrum of finite-size RNNs. To verify this claim, we simulated a network with 1000 neurons ($N=1000$) and showed that it produces a limit cycle similar to that shown in Fig. 1D and the 400 largest eigenvalues of its covariance eigenspectrum closely match our theoretical prediction. This new result will be mentioned in the main text. \n### Weakness 3\nWe agree with the reviewer. The sentence preceding the statement of the conjecture, “Known results and simulations suggest the following conjecture”, was too concise. We will replace this sentence by \n> Drawing intuition from Fourier analysis, the smoothness of a function (here the kernel) should be related to the decay rate of its Fourier transform (here the eigenspectrum)—the smoother the function, the faster the decay rate of its Fourier transform. Known results on the eigenvalues of random feature kernels for the cases $p=0$ and $p=1$, with $c=0$, confirm this intuition and show how it extends to general integers $d$ (Bach JMLR 2017). Extrapolating those results to any nonnegative $p$ and any real $c$, we get the following conjecture:\n### Question 1 \nThe paper starts with the solvable RNN model for two reasons. \n- First, the RNN model highlights the nontriviality of the notion of “pre-activation dimensionality” which this paper aims to introduce. The model shows that the autonomous dynamics (i.e. dynamics that are not input-driven) of an RNN can be both low-dimensional in the space of pre-activations and high-dimensional in the space of post-activations, which is a somewhat counterintuitive fact. This fact provides a mechanistic motivation for the notion of “pre-activation dimensionality”, as opposed to the more phenomenological notion of intrinsic dimensionality often considered in neuroscience. To make this point clearer, we will add the following paragraph in the introduction:\n> Here, we first construct a solvable RNN model that reconciles the low- and high-dimensional perspectives on population activity by carefully disambiguating the *linear* dimension of the system *before* and *after* the neurons' nonlinearity, which we refer to as the pre- and post-activation dimension, respectively. This dichotomy refines the usual distinction between linear and “intrinsic” dimension (Stringer et al. Nature 2019; Jazayeri and Ostojic CONB 2021; Humphries NBDT 2021), since the intrinsic dimension of a system is the same before and after any continuous, injective nonlinearity. Using the notions of pre- and post-activation linear dimensions, we show that our RNN can be exactly reduced to a low-dimensional dynamical system in the space of pre-activations, making the pre-activations low-dimensional. Then, we show that these latent dynamics generate high-dimensional post-activation activity that has a power-law covariance eigenspectrum. (In this work, dimension will always refer to linear dimension, unless stated otherwise.)\n\n- Second, the RNN model we present in Section 2 constitutes an original contribution to the theory of large neural networks, in and of itself. While, in this paper, we use the model to reconcile different views on dimensionality in systems neuroscience, this model should also be of interest to researchers in computational/theoretical neuroscience, machine learning, and the physics of complex systems, who are interested in the dynamics of RNNs.\n### Question 2 \nSee answer to Weakness 1(i) above.\n### Question 3\nNCE, CEBRA, and Rastermap are three different dimensionality reduction methods––the latter two being primarily for exploratory analyses––and are therefore hard to compare (e.g., it is difficult to compare CEBRA and Rastermap, hence (Stringer et al. Nat. Neuro. 2025) does not contain such comparisons). The purpose of NCE is to estimate the pre-activation dimensionality of neuronal activity and it tries to do so with as little inductive bias as possible. To our knowledge, there is no other latent variable model that is specifically designed to estimate pre-activation dimensionality. Methods such as LFADS and its derivatives could be used to estimate the pre-activation dimensionality but they come with complex inductive biases that are not necessary to the task of estimating pre-activation dimensionality. How NCE compares to methods such as LFADS (or RADICaL), which have more inductive bias, when the number of recorded neurons is large is an interesting question that lies beyond the scope of this work.\n### Question 4 \nIn Section 3, we assume that the inputs are uniformly distributed on the sphere. This is an assumption used for mathematical convenience only and it is probably not biologically realistic. Interestingly, a recent paper by (Li et al. JMLR 2024) (cited in the text) shows that results on the decay rate of the eigenspectrum of random feature kernels with uniformly distributed inputs on the sphere can be generalized to a broad class of distributions on more general domains. We will briefly mention this in the main text.\n\nIn G.3 of the Supplementary Material, we show that NCE can infer latents that are correlated with running. Since these latents (and running) are not uniformly distributed, this suggests that NCE performs well even when the latents are not uniformly distributed."}, {"review_id": "srEUQUkLwx", "reviewer": "Reviewer_zfNk", "summary": "The contributions of the paper are twofold: First, it presents an analytically solvable recurrent neural network model that combines an efficiently low-dimensional dynamics of the neuron activations in the pre-activation space (activations inside the output function) and a very high-dimensional dynamics with a power-law eigenspectrum for the post activations (post the output threshold). The eigenspectrum of the post activations is derived for different types of low-dimensional solutions for the pre-activations, and RKHS theory a relationship beween the form of the activation function (rectified power function) and the decay exponent of the spectrum is derived and numerically verified. \n\nThe second contribution is an analysis the the estimated pre-activation dimensionality of mouse visual cortex responses using a neural cross-encoder architecture. This analysis shows low dimensionality for grating and high dimensionality for natural image responses.", "questions": "What exactly of the model in (13) is fitted to the neural activities?\n\nApparently sloppy notation in L. 117 ff: For even m the eigenvalues should be always 0. Should not for the uneven ones sin (m \\theta) be the eigenfunction?\n\nL. 957 (appendix).  should second z_s be z_t ?\n\nIs the factor (2/ \\pi) in the equation below L. 921 correct?\n\n(17) =>  C_{ij} = 4 * [term in (18)] - 4 * (1/4)", "limitations": "Limitations are thoroughly discussed,", "rating": 5, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 3, "significance": 3, "originality": 3, "strengths_and_weaknesses": "The mathematical analysis part of the paper is very strong and combines methods from probabilistics, neural field theory, and functional analysis / RKHS theory to derive and analyze the solutions of recurrent networks with very specific assumptions for the lateral coupling structure. \n\nIt is not entirely clear in how far the results break down for more general forms of lateral connectivity (different width of the interaction kernel, deviation form cos shape). Generality in presence of such changes would be required for a validity of the model for larger classes of real neurons.", "was_revised": true, "final_justification": "Author's feedback is plausible and addresses as far as possible my concerns. I thus support the publication of this paper.", "rebuttal": "The reviewer highlights the strengths of the theory we present in the paper and they make detailed and accurate comments on the mathematical parts of the paper, which are very much appreciated.\n\n### Answers to Strengths and Weaknesses\nThe reviewer asks whether the solvable model can be generalized to lateral connectivity that is not of cosine form, Eq. (2). This is an excellent question. In fact, we can get a very similar solvable model with a more localized connectivity profile by replacing the cosine connectivity, Eq. (2), by a connectivity of the form: \n$$W_{ij}=J\\sum_{l=1}^4 \\beta_l \\cos(l(\\theta_i - \\theta_j) - \\Delta).$$\nThis extended model is also solvable and can produce a limit cycle on the circle in the space of pre-activations while exhibiting high-dimensional post-activations, as the original model. The analysis of this extended model will be presented in a new section in the Supplementary Material and briefly mentioned in the main text.\n\n\n### Answers to Questions\n**Q: What exactly of the model in (13) is fitted to the neural activities?** \n\nA: The parameters that are learned are the activation parameter $p$ (that is shared across target neurons), the vector of biases $\\mathbf{c}$, the matrix of readout weights $\\mathbf{u}$, and a vector of baseline firing rates $\\mathbf{r}$. This is stated in the sentence following Eq. (13). The latent variables $\\mathbf{z}$ are the output of the encoder and the input of the decoder; hence, they are inferred but not learned. To clarify the notation, we will update the text so that $\\tilde{b}$ denotes the recorded activity of target neurons, and $b$ denotes the NCE reconstruction of these neurons.\n\n**Q: Apparently sloppy notation in L. 117 ff: For even m the eigenvalues should be always 0. Should not for the uneven ones sin (m \\theta) be the eigenfunction?**\n\nA: There was indeed a minor mistake in Eq. (8). The expression in Eq. (8) will be replaced by $$\\lambda_n = \\frac{4}{\\pi^2}\\left(2 \\left\\lfloor \\frac{n-1}{2} \\right\\rfloor + 1\\right)^{-2}, \\forall n \\in \\mathbb{N}_*.$$ \n\n**Q: L. 957 (appendix). should the second z_s be z_t ?**\n\nA: The reviewer is right, we will correct this typo.\n\n**Q: Is the factor (2/ \\pi) in the equation below L. 921 correct?**\n\nA: The factor is correct but there were indeed some factors $1/(2\\pi)$ missing in ll. 915-920. We will correct these typos. \n\n**Q: (17) => C_{ij} = 4 * [term in (18)] - 4 * (1/4)**\n\nA: We will do the recommended modification, which helps readability.  \n\n### Final remarks\nWe hope that the strength of our theoretical results and their relevance to several disciplines related to the theory of neural systems, the fact that our solvable model can be extended to lateral connectivity that goes beyond the cosine form, and the mathematical clarification/corrections we have made (thanks to the reviewer’s detailed and accurate comments) will prompt them to reconsider their rating."}, {"review_id": "UbIo21PQ1I", "reviewer": "Reviewer_YAj9", "summary": "The authors provide theoretical proof using an analytically solvable recurrent neural network model that high-dimensional observations can be generated from low-dimensional latent variables. Further, they demonstrate that the generated high-dimensional observations can have a long-tailed covariance eigenspectrum, as is commonly seen in real neural data. Finally, the authors present a novel nonlinear latent variable model for calcium recordings, called Neural Cross-Encoder, and use it to demonstrate that high-dimensional neural activity in some regimes (i.e., spontaneous activity; response to natural images) can be described using low-dimensional latents while in other regimes they cannot (i.e., response to natural images).", "questions": "1. It would be helpful if the authors could more clearly communicate that the analyses provided in section 4 are not related to the theoretical work in section 2 (i.e., solvable RNN model), but rather only for supporting / validating the claims made in section 3 (i.e., the relationship between pre-activation dimensionality, activation function, and post-activation eigenspectrum). Also, it would be informative/interesting if the authors visualized the covariance eigenspectrum of the true vs. model-predicted neural activity. As the authors pointed out, the goal here isn’t to estimate (or predict/recreate) the tail of the eigenspectrum, but it would be interesting to see how the decay changes as a function of dimensionality + activation function in real data.\n2. How are the pre-activation dimensionality results in section 4 novel / surprising / noteworthy? First, natural images are statistically complex stimuli and so it is not entirely surprising that their latent representations are of higher dimensionality than the drifting gratings. Second, how do we know that the results are not a byproduct of the model class used in NCE? It might be that for more complex stimuli, which require more complex neural encoding, the current encoder architecture choice is not sufficiently descriptive to learn the complex neural encoding [1]. Third, comparing against only a linear activation function seems insufficient because calcium recordings are known to be nonlinear. It would be interesting / more compelling if the authors explored using different activation functions and looked at the impact on explained variance and the covariance eigenspectrum. The analyses provided in section 3 only considered general rectified power activation functions so it would be more compelling / demonstrative if the authors compared the performance of NCE within this family of functions (not just against linear, i.e., reduced rank regression).\n3. It would be helpful if the authors could provide additional information in the Appendix about which specific regions were recorded and roughly how many neurons in each of the regions were determined “responsive”. Were neuron splits between source and target region agnostic (the authors stated that they collected data from primary and higher visual cortices)? It would be helpful if authors could also provide information on how many neurons and from which regions were assigned to source and target. This information may be important towards interpreting the results for each of the three task conditions, especially gratings vs natural images.\n4. What was the encoder architecture used in NCE? Please add this to the Appendix.\n5. Is the analysis presented in section 2 only limited to neurons on the unit circle/sphere manifold? Would the authors expect that this analysis could generalize to other manifolds?\n\n\n**Minor**\n- Could the authors provide an explanation as to why specifically for the spontaneous condition multiple timepoints were used as input for estimating the latent variable as a single timepoint? Was something similar attempted for the stimulus conditions as well? Why/why not? If yes, what were the results?\n- Was there a particular reason why the authors first pretrained the NCE model using knowledge distillation rather than directly training on the target neurons?\n- Equation (2), how did the authors choose the values for $J$ and $\\Delta$?\n- Could the authors clarify the following: “...per experiment, across 7 experiments (3 spontaneous activity, 3 natural images, 1 drifting gratings, across three mice)...” Does this mean 7 experiments were conducted for each of the 3 mice or across the 3 mice aggregated (i.e., only 1 mouse had the drifting gratings stimulus or all 3).\n- Fig. 3F would benefit from a few more y-axis labels for the points not close to the extrema.\n- In line 101, minor typo that it should be a pair of neurons $i$ and $j$. Also minor typo line 905 in Appendix: “standard deviation [of]”. Also, the equations following lines 980 and 981 in Appendix E have a reversed bracket in the Appendix, e.g., $\\theta \\sim \\mathcal{U}(\\left[0, 2\\pi\\right])$\n- In Appendix line 1044 (last line of page 27), it would be useful to remind readers of what $\\mathcal{F}$, $a$, and $b$ denote, especially as $ \\mathcal{F} $ and $ a $ were not defined.\n- Clarification for the equations after line 897 of Appendix A: is that an inner product or elementwise product in the integral: $ \\begin{pmatrix}\\cos(\\theta) \\sin(\\theta)   \\end{pmatrix} \\cdot  \\begin{pmatrix}\\cos(\\mathcal{v}) \\sin(\\mathcal{v}) \\end{pmatrix} $?\n\n\n**References**\n\n[1] D.L.K. Yamins, H. Hong, C.F. Cadieu, E.A. Solomon, D. Seibert, & J.J. DiCarlo. Performance-optimized hierarchical models predict neural responses in higher visual cortex, Proc. Natl. Acad. Sci. U.S.A. 111 (23) 8619-8624, https://doi.org/10.1073/pnas.1403112111 (2014).", "limitations": "Limitations with respect to the experimental results in section 4 and the NCE model architecture are mostly sufficient. I recommend the following revision: “The fact that we were not able to accurately predict neuronal responses to natural images with a low-dimensional NCE model does not necessarily exclude the possibility for such a model to exist, and one could possibly find it with a larger training set **or a different model class/hyperparameter configuration**.” This is because the proposed NCE model may not have sufficient descriptive power/modeling capacity to robustly identify the shared variance in populations of neurons, especially in the context of more semantically complex stimuli wherein the encoding is less clear / more distributed [1,2].\n\nAlso, depending on the answer to my 5th question above, it may also be worth explicitly mentioning that the theoretical analysis provided in section 2 is limited to only specific manifold structures (e.g., unit ring/sphere).\n\n**References:**\n\n[1] Bolaños, F., Orlandi, J.G., Aoki, R. et al. Efficient coding of natural images in the mouse visual cortex. Nat Commun 15, 2466 (2024). https://doi.org/10.1038/s41467-024-45919-3\n\n[2] Michael Weliky, József Fiser, Ruskin H Hunt, David N Wagner. Coding of Natural Scenes in Primary Visual Cortex. Neuron 37(4), 703-718, 2003.", "rating": 5, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 3, "significance": 2, "originality": 3, "strengths_and_weaknesses": "**Strengths**\n- The theoretical/analytical work presented in sections 2 and 3 are the strongest part of the manuscript. Proofs and derivations (not closely checked) seem complete and well-supported by the cited prior work.\n\n- Proofs provided in the appendix were helpful for understanding the main results / claims in the manuscript.\n\n**Weaknesses**\n- The relevance of the experimental analyses to the theoretical work in sections 2 and 3 is not particularly clear. Moreover, the results/claims in section 4 (i) are not particularly novel / surprising, (ii) can be improved with additional experimental details and a few more/modified analyses, and (iii) do not necessarily help prove / support the theoretical claims made previously – as currently presented. Please see questions below for specific recommendations and points of confusion where additional clarification is needed.", "was_revised": true, "final_justification": "The authors present novel theoretical work that aims to address a very important question in neuroscience: the relationship between high-dimensional observable neural activity and the low-dimensional latent dynamics that drive it. In addition to complete derivations and proofs, the authors also provide relatively complete empirical demonstrations supporting their theoretical work. Moreover, the authors provided ample ablations, clarifications, and explanations during the rebuttal period that better demonstrated the robustness of their proposed model and its ability to interrogate the true pre-activation dimensionality in the measured system. Although I believe certain open questions remain (i.e., extending the theoretical work towards more general manifold definitions, demonstrating how well their model can reconstruct the covariance eigenspectrum of the observations), the work presented here will be an important contribution to the field. I recommend the paper for acceptance.", "rebuttal": "We thank the reviewer for the high quality of their review and for highlighting the strengths and clarity of our theoretical results. Below are our answers to their questions.\n### Question 1 \nWe understand what the reviewer means when they challenge the link between Sections 2 and 4 but we would like to argue that there is a clear line of reasoning connecting them. \n\nSpecifically, we use the experimental data to provide examples confirming the theoretical result that high post-activation dimension can arise from either low or high pre-activation dimension: drifting gratings and spontaneous activity (high post-activation dimension, low pre-activation dimension) vs. natural image responses (high post- and pre-activation dimension).\n\nTo clarify this line of reasoning, we have added the following piece of text in the introduction:\n>Here, we first construct a solvable RNN model that reconciles the low- and high-dimensional perspectives on population activity by carefully disambiguating the *linear* dimension of the system *before* and *after* the neurons' nonlinearity, which we refer to as the pre- and post-activation dimension, respectively. This dichotomy refines the usual distinction between linear and “intrinsic” dimension (Stringer et al. Nature 2019; Jazayeri and Ostojic CONB 2021; Humphries NBDT 2021), since the intrinsic dimension of a system is the same before and after any continuous, injective nonlinearity. Using the notions of pre- and post-activation linear dimensions, we show that our RNN can be exactly reduced to a low-dimensional dynamical system in the space of pre-activations, making the pre-activations low-dimensional. Then, we show that these latent dynamics generate high-dimensional post-activation activity that has a power-law covariance eigenspectrum. (In this work, dimension will always refer to linear dimension, unless stated otherwise.)\n\n>Before analyzing experimental recordings, we revisit the spectral theory of infinite-width neural networks (random feature kernels) [Refs.] to quantitatively relate the pre-activation dimension, the neuronal activation function, and the post-activation covariance eigenspectrum. This three-way relationship tells us that high-dimensional activity is consistent with both low- and high-dimensional pre-activations. [...]\n\nHence, we need latent variable models to infer the pre-activation dimensinality (Section 4).\n\nRegarding the second part of Question 1, we agree with the reviewer that it would be interesting to see how the eigenspectrum decay rate changes with the dimensionality of the latents and neuronal nonlinearities on NCE fitted to neural data. Unfortunately, we are not able to do this analysis at this stage because we do not have a method for estimating the decay rate of the eigenspectrum that is reliable and accurate enough, when the data is noisy and the number of neurons is relatively small (the number of target neurons is 1,000).\n### Question 2 \n#### (First) \nThe idea that neural activity is low-dimensional is dominant in the neuroscience literature (e.g. Gao & Ganguli CONB 2015, Gallego et al Neuron 2017 and refs. within).  Recent observations of high post-activation dimension have challenged this view, but leave open the possibility that pre-activation dimension could be high or low. \n\nFor example, whether spontaneous activity in visual cortex has high or low pre-activation dimensionality is not known: the study of (Stringer et al. Science 2019) and the recent work of (Manley et al. Neuron 2024) stress the high-dimensional aspect of spontaneous activity and do not suggest the pre-activations could be low-dimensional. Hence, the low-dimensional pre-activations we find during spontaneous activity are novel. We consider low pre-activation spontaneous dimensionality our main empirical finding. \n#### (Second)\nAn important difference between this work and previous studies (eg Yamins et al) is that we are predicting neurons from neurons, rather than from images.  Yamins and others use deep networks to process images, but a single layer to predict neurons from latents. \n \nWe will add to the Limitations the suggested text.  Also, to demonstrate the robustness of results to architecture changes, we present the fraction of explainable variance explained for two example datasets as we vary the number of layers, and the size of the first layer. For multi-layer encoders, each successive layer is smaller by a factor of two.\n\n**Drifting gratings (d = 4):**\n\n| Layer size \\ # Layers | 3     | 2     | 1     |\n|-----------------------|-------|-------|-------|\n| 1000                  | 0.97  | 0.97  | 0.96  |\n| 500                   | 0.97  | 0.97  | 0.94  |\n| 250                   | 0.95  | 0.96  | 0.94  |\n\n**Natural images (d = 32):**\n\n| Layer size \\ # Layers | 3     | 2     | 1     |\n|-----------------------|-------|-------|-------|\n| 1000                  | 0.91  | 0.89  | 0.89  |\n| 500                   | 0.91  | 0.92  | 0.89  |\n| 250                   | 0.92  | 0.91  | 0.91  |\n\nFurthermore, because we use the same model architecture across all experimental conditions (drifting gratings, natural images, and spontaneous activity), the comparison between these three conditions says something about the data that is not simply a byproduct of the model architecture. In addition, the fact that NCE always performs above the linear baseline excludes trivial cases of severe under- or over-fitting, which could have led to misleading results.\n#### (Third) \nWe followed the review’s suggestion and performed an additional analysis comparing NCE with varying fixed activation parameters $p$. We focused on the case of drifting gratings with $d=4$ and fitted NCE with $p$ ranging from $0.5$ to $4$. The results are summarized in the table below.\n\n| $p$   | R²    |\n|-------|--------|\n| 0.5   | 0.155  |\n| 0.75  | 0.161  |\n| 1.0   | 0.165  |\n| 1.25  | 0.167  |\n| 1.5   | 0.167  |\n| 1.75  | 0.169  |\n| 2.0   | 0.169  |\n| 3.0   | 0.164  |\n| 4.0   | 0.110  |\n\nThese results indicate that within the standard range $p\\in[1,2]$ (the fitted $p$ is 1.89 in the original analysis), the R2 is relatively stable. This numerical experiment suggests that the estimation of the pre-activation dimensionality should not be too sensitive to $p$, as long as $p$ is within the standard range $p\\in[1,2]$. \n### Question 3 \nWe will add to the appendix:\n>The recordings covered part of primary visual cortex, higher visual areas AM, PM, LM, and RL, with occasional inclusion of some somatosensory areas. We confirmed the recorded area via retinotopic mapping using sparse noise stimuli. To determine stimulus responsiveness, we used the signal-related variance metric to determine what fraction of a neuron’s variance could be explained by stimuli (described in Stringer et al. Nature 2019). In a typical recording with 22,000 neurons, the top 1000 cells had 49% stimulus-related variance, the next 1000 had 22%, 14%, 10% and so on, where the 10,000th cell had ~3%. The responsive neurons were distributed across primary and higher visual cortices, and the selection of source and target neurons was agnostic to location.\n\nWe agree that comparing the latent dimensionality of responses across different visual areas is an interesting question, though it is outside of the scope of this work.\n### Question 4 \nWe will add to the appendix:\n>The encoder is a 3-layer feedforward network with ReLU activation functions. The hidden layers contain (500, 250, 100) units.\n### Question 5 \nWe believe that the general idea of having a recurrent network with low-dimensional dynamics but high-dimensional activity is a general phenomenon that is not specific to the circle manifold we consider. However, finding a model whose latent dynamics and covariance eigenspectrum can be analytically solved, that is not a ring model, remains an open theoretical problem. \n### Minor comments \n1. The stimulus responses are the mean activity over a response window of 2.0s for gratings and 0.8s for natural images, and this averaging changes the statistics of the activity and helps reduce the shot/independent noise. To have a similar effect in spontaneous activity, we chose to use multiple timepoints in the source neurons to better estimate the latents. \n2. We found that the knowledge distillation method allowed the NCE to learn in fewer steps.\n3. We chose $J$ and $\\Delta$ such that, in Eq. (5), the expression for the latent dynamics is as simple as possible. \n4. The experiment numbers reported are aggregated across all mice. Since the submission, we have expanded the analysis to additional sessions of drifting gratings, such that we now have three mice with at least one session of spontaneous activity, images, and gratings, which has confirmed the results presented.\n5. We will add intermediate ticks to 3F.\n6. We will correct these typos.\n7. We will update the notation to include $\\tilde{b}$ to denote the recorded neural activity of the target neurons, and we will clarify the definitions of $a$ as the recorded activity of source neurons and $b$ as the predicted activity of source neurons. Hence the expression for loss in 1044 will be updated to $||b - \\tilde{b}||^2$ where $b = \\phi(u \\mathcal{E} (a) + c) + r$ and we will remove $\\mathcal{F}$. \n8. It is an inner product. We will clarify this in the Appendix.\n### Final remarks\nIn summary, the experimental results demonstrate the theoretically-predicted consistency of either low- or high- preactivation dimension with high post-activation dimension. The low spontaneous pre-activation dimensionality our method finds is a novel finding. We have also provided additional evidence showing that in the neuron-to-neuron prediction setting, the complexity of the encoder may not be a limiting factor since our results are robust to changes of encoder architectures. In light of the clarification we have provided and the additional analyses we have performed, we hope that the reviewer will reconsider their rating."}, {"review_id": "yhYw7WToj1", "reviewer": "Reviewer_Xtzt", "summary": "This paper introduces theoretical models and empirical results supporting the hypothesis that high-dimensional neural activity can arise from a low-dimensional latent space combined with nonlinear activations. The paper begins by first presenting a RNN model whose dynamics can be explicitly solved. This allows for direct analysis of the covariance eigenspectrum of the model’s post-activation activity, which the paper shows to be high-dimensional. The paper then relaxes assumptions on this solvable model and proposes a conjecture about the decay exponent of the covariance eigenspectrum in this relaxed model. Finally, the paper develops the Neural Cross-Encoder (NCE), a nonlinear autoencoder-inspired architecture, which they apply to neural recordings on mice during a visual task. They find that the NCE distinguishes between response types with low and high-dimensional pre-activations.", "questions": "- In Fig 3I, it seems that the linear model is overfitting for large $d$ while the NCE does not. Why do you think this is the case, especially since the NCE is more flexible than the linear model?\n- How does the NCE deal with noise in the observed neural recordings and still extract seemingly meaningful latent variables (at least for the drifting gratings and spontaneous response types)?", "limitations": "The authors address limitations about potential difficulties when fitting NCE, and in particular they point out that even if NCE is unable to find a well-fitting model, this does not mean such a model doesn't exist. I agree with this limitation and I think it could be better incorporated into the interpretation of results (section 4.3) as well. For example, on line 217 the authors write that the results on natural images suggest that natural images produce high-d representations in the space of pre-activations. However, according to the limitation, this is not necessarily the case -- it could also be due to getting stuck in local optima during training.", "rating": 5, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 2, "clarity": 3, "significance": 4, "originality": 3, "strengths_and_weaknesses": "Strengths: This paper tackles an interesting and fundamental question about how we should interpret and test for latent dimensionality in neural population activity. Overall, the paper is easy to read and the technical steps are clearly explained. The solvable RNN model presented in Section 2 is relatively simple yet illustrative of their hypothesis that post-activation activity may still be generated from a low-d pre-activation latent space. Their overarching question and empirical results in Section 4 are well-motivated by real findings from visual cortical population activity. \n\nWeaknesses: My main concern with this paper is about the cohesiveness between the theoretical results presented in Section 2-3 and the empirical results presented in Section 4. Right now, the paper reads as three separate contributions across these sections, and it is not clear to me how the theoretical models allow us to reinterpret the notion of dimensionality in neural recordings. I would have liked to see some experiments about how well Conjecture 1 holds for fitted NCE models of different dimensionalities $d$ in Section 4. Otherwise, I am not sure about the role of Section 3 in this paper: the conclusion that high-d activity might arise from a nonlinear transformation of low-d states is already apparent from the solvable RNN model. \n\nI also have some concerns about the Neural Cross-Encoder methodology. In 4.2, the authors state that the NCE has a 3-layer encoder to allow for flexible estimation of latent variables; however, the choice of 3 layers is a bit arbitrary and might also lead to overfitting in some cases. I think the authors should address how one would decide the number of encoder layers / flexibility of the encoder. It is also not clear to me why the NCE is better at discarding variability that is not shared across neurons compared to a standard autoencoder fit on the entire set of neurons, since that should also discard such variability. I would have liked to see a more detailed explanation to motivate NCE as well as empirical comparisons with the standard autoencoder. \n\nFinally, some suggestions on clarity:\n- On line 135 the authors introduce the notation $\\mathbf{z}$, and it would help to draw the connection to $\\mathbf{\\kappa}$ from the previous section.\n- Line 200 typo in “variability”", "was_revised": true, "final_justification": "The authors have sufficiently addressed my concerns and questions during the rebuttal period. I appreciate the additional text in the introduction which clarifies the link between the theoretical and empirical contributions of the paper. I also am pleased to see the additional experimental results showing that the NCE is quite robust to varying the layer size, number of layers, and numbers of cells. The authors also provided some helpful insight into the differences between the NCE and standard autoencoder, as well as the NCE's ability to extract latents given noisy observations. Thus, I chose to raise my score to a 5.", "rebuttal": "The reviewer gives an accurate and comprehensive summary of the paper, stating that we tackle “an interesting and fundamental question”. \n### Answers to the weaknesses\n#### 1. Lack of cohesiveness between the theoretical results presented in Section 2-3 and the empirical results presented in Section 4. \n\nWe use the experimental data (section 4) to provide examples confirming the theoretical result (sections 2,3) that high post-activation dimension can arise from either low or high pre-activation dimension: drifting gratings and spontaneous activity (high post-activation dimension, low pre-activation dimension) vs. natural image responses (high post- and pre-activation dimension).\n\n\nTo clarify the link between the three result sections, we have added the following text in the introduction\n\n>Here, we first construct a solvable RNN model that reconciles the low- and high-dimensional perspectives on population activity by carefully disambiguating the *linear* dimension of the system *before* and *after* the neurons' nonlinearity, which we refer to as the pre- and post-activation dimension, respectively. This dichotomy refines the usual distinction between linear and “intrinsic” dimension (Stringer et al. Nature 2019; Jazayeri and Ostojic CONB 2021; Humphries NBDT 2021), since the intrinsic dimension of a system is the same before and after any continuous, injective nonlinearity. Using the notions of pre- and post-activation linear dimensions, we show that our RNN can be exactly reduced to a low-dimensional dynamical system in the space of pre-activations, making the pre-activations low-dimensional. Then, we show that these latent dynamics generate high-dimensional post-activation activity that has a power-law covariance eigenspectrum. (In this work, dimension will always refer to linear dimension, unless stated otherwise.)\n\n>Before analyzing experimental recordings, we revisit the spectral theory of infinite-width neural networks (random feature kernels) [Refs.] to quantitatively relate the pre-activation dimension, the neuronal activation function, and the post-activation covariance eigenspectrum. This three-way relationship tells us that high-dimensional activity is consistent with both low- and high-dimensional pre-activations. [...]\n\n#### 1b. I am not sure about the role of Section 3 in this paper\n\nSection 2 provided one example of a phenomenon described more generally in section 3: determining the pre-activation dimensionality from post-activation eigenspectrum is an ill-posed problem.  We add to this a specific conjecture relating three concepts: the post-activation and pre-activation dimensionalities, and neuronal nonlinearities.  While this conjecture is not required for Section 4, we believe it is of sufficient general interest to warrant inclusion. \n#### 2. How well does Conjecture 1 hold for fitted NCE models of different dimensionalities in Section 4?\n\nThis is a natural question to ask. However, answering this question would require a reliable and accurate method for estimating the eigenspectrum decay rate of noisy (and finite) data and we have not found a satisfying method for that yet. This prevents us from answering this question. \n#### 3. 3 layers is a bit arbitrary and might also lead to overfitting in some cases. I think the authors should address how one would decide the number of encoder layers / flexibility of the encoder. \n\nWe agree that the current justification included in the manuscript for the selection of hyperparameters, i.e. the number/width of encoder layers and the number of source/target neurons is not sufficient. We present additional data to justify these decisions.\n\nThe results are only mildly sensitive to the depth and width of the encoder, which were selected via grid search on an example dataset. To demonstrate the minor effect of encoder architecture, we present the fraction of explainable variance explained for two example datasets with as we vary the number of layers, and the size of the first layer. For multi-layer encoders, each successive layer is smaller by a factor of two. \n\n**Drifting gratings (d = 4):**\n\n| Layer size \\ # Layers | 3     | 2     | 1     |\n|-----------------------|-------|-------|-------|\n| 1000                  | 0.97  | 0.97  | 0.96  |\n| 500                   | 0.97  | 0.97  | 0.94  |\n| 250                   | 0.95  | 0.96  | 0.94  |\n\n**Natural images (d = 32):**\n\n| Layer size \\ # Layers | 3     | 2     | 1     |\n|-----------------------|-------|-------|-------|\n| 1000                  | 0.91  | 0.89  | 0.89  |\n| 500                   | 0.91  | 0.92  | 0.89  |\n| 250                   | 0.92  | 0.91  | 0.91  |\n\nVarying the number of source and target neurons up to 2000 does not substantially affect the results shown in Figure 3. We selected 500 and 1000 as intermediate values for computational reasons, since we fit 13 models with 5 shuffles each in 10 sessions for a total of 650 models (we have added 3 more experimental sessions since the original submission, confirming and strengthening results). To illustrate this stability, in the table below, we show the fraction of explainable variance explained for a model with $d=4$ on drifting grating responses, for NCE and RRR models, across combinations of source and target cells.\n\n| N_target \\ N_source | 100               | 500            | 1000           | 2000           |\n|---------------------|------------------|----------------|----------------|----------------|\n| 100                 | NCE: 0.98 / RRR: 0.92 | 0.96 / 0.88   | 0.96 / 0.86   | 0.97 / 0.86   |\n| 500                 | 0.98 / 0.91       | 0.97 / 0.85     | 0.94 / 0.84     | 0.93 / 0.84     |\n| 1000                | 0.98 / 0.89       | 0.97 / 0.84     | 0.94 / 0.83     | 0.91 / 0.83     |\n| 2000                | 0.99 / 0.90       | 0.97 / 0.85     | 0.97 / 0.87     | 0.92 / 0.84     |\n\nIncreasing the number of cells beyond 2000 somewhat reduces performance and we suggest two potential explanations. First, due to ethical constraints our experiments are of limited duration, hence a typical session produces ~1000 training examples. This may not be enough to train a relatively large encoder network with >1000 input neurons. Second, only a fraction of neurons recorded have reliable stimulus responses, and including more cells introduces a large amount of non-stimulus-related variance.\n\nOf the remaining hyperparameters, we have found that the performance is not sensitive to variations in learning rate, batch size, and momentum by exploring the hyperparameter space with the Optuna optimizer (Akiba et al. KDD 2019). \n\nWe will repeat the measurements shown in the tables above across multiple sessions and different values of $d$, and include them as supplementary figures in the appendix. We will also include in the appendix the details of hyperparameter selection, described above.\n#### 4. Why is NCE better at discarding variability that is not shared across neurons compared to a standard autoencoder fit on the entire set of neurons, since that should also discard such variability?\n\nA standard autoencoder would discard non-shared variability only when the bottleneck $d$ is small. When $d$ is large (without suitable regularization), nothing prevents an autoencoder from accounting for non-shared variability. For example, if $d=N$, an autoencoder could learn a function that is close to the identity function (see Steck NeurIPS 2020 and references therein). Because we want to look at the effect of varying $d$, we need to exclude the tendency of a standard autoencoder to predict independent variability better with larger $d$, which precludes estimating the pre-activation dimensionality. \n\nBy design, NCE discards variability that is not shared across neurons, since it predicts the activity of target neurons from a disjoint set of source neurons, which prevents learning the identity function.\n#### Suggestions on clarity \non l. 135, we will replace “denoted $\\mathbf{z}$” by “henceforth denoted by $\\mathbf{z}$ (instead of $\\boldsymbol{\\kappa})$” and we will correct the typo on l. 200.\n### Answers to questions\n1. The iterative optimization procedure for the NCE described in the Appendix, where the model with the best validation loss is selected at each phase, has a similar effect to the early stopping method in preventing overfitting. We believe this prevents the NCE from overfitting to the noise in the data, whereas the L2 regularization is not sufficient in preventing overfitting in high-d linear models. Note that when training the linear model, at each d we train with 6 different values for the regularization parameter, and we have selected the range such that the selected models are not at the extrema, so increasing the L2 regularization does not prevent the overfitting. Note that all reported scores are on a held-out test set.\n2. This is a good question because the encoder in NCE is taking noisy inputs (the activity of the source neurons). Our best guess for how NCE can extract seemingly noiseless latents is the following: The use of cross-prediction during training strongly encourages NCE to learn encoder parameters leading to the inference of latents representing shared and distributed axes of variability across the population (in the space of pre-activations). To infer these latents, the weights in the encoder have to be distributed. In turn, such distributed weights help the encoder to produce outputs where input noise is averaged out.\n\n### Limitations\nRegarding the reviewer’s comment under Limitations, we will add, on l. 218, “(but see Limitations below)”.\n\n### Final remarks\nIn conclusion, we hope that the text we added to the introduction, which should improve the cohesiveness of the paper, the additional analyses and the explanation we provide on the choice of NCE hyperparameters, and the clarification of the difference between NCE and a standard autoencoder will prompt the reviewer to reconsider their rating."}]

Author rebuttal:

Meta-review: Since Stringer et al's 2019 paper, there wasn't a satisfactory explanation of how the relationship between high-dimensional observable neural activity and the low-dimensional latent dynamics that drive it in early visual system. This paper presents a very strong mathematical analysis that combines methods from various subdisciplines to derive and analyze the solutions of recurrent networks with very specific assumptions for the lateral coupling structure. The low-dimensional pre-activation and the power-law post-activation as a possible resolution represents a significant progress and opens future research directions. Authors also present a novel methodological contribution called the NCE to analyze the high-dimensional deconvolved neural responses.
Decision: Accept (oral)