High-dimensional neuronal activity from low-dimensional latent dynamics: a solvable model

The reviewer gives an accurate and comprehensive summary of the paper, stating that we tackle “an interesting and fundamental question”.

Answers to the weaknesses

1. Lack of cohesiveness between the theoretical results presented in Section 2-3 and the empirical results presented in Section 4.

We 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).

To clarify the link between the three result sections, we have added the following text in the introduction

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.)

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. [...]

1b. I am not sure about the role of Section 3 in this paper

Section 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.

2. How well does Conjecture 1 hold for fitted NCE models of different dimensionalities in Section 4?

This 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.

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.

We 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.

The 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.

Drifting gratings (d = 4):

Natural images (d = 32):

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.

Increasing 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.

Of 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).

We 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.

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?

A 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.

By 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.

Suggestions on clarity

on 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.

Answers to questions

1. 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.
2. 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.

Limitations

Regarding the reviewer’s comment under Limitations, we will add, on l. 218, “(but see Limitations below)”.

Final remarks

In 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.

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.

Question 1

We 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.

Specifically, 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).

To clarify this line of reasoning, we have added the following piece of text in the introduction:

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.)

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. [...]

Hence, we need latent variable models to infer the pre-activation dimensinality (Section 4).

Regarding 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).

Question 2

(First)

The 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.

For 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.

(Second)

An 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.

We 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.

Drifting gratings (d = 4):

Natural images (d = 32):

Furthermore, 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.

(Third)

We 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.

These 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]$.

Question 3

We will add to the appendix:

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.

We 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.

Question 4

We will add to the appendix:

The encoder is a 3-layer feedforward network with ReLU activation functions. The hidden layers contain (500, 250, 100) units.

Question 5

We 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.

Minor comments

1. 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.
2. We found that the knowledge distillation method allowed the NCE to learn in fewer steps.
3. We chose $J$ and $\Delta$ such that, in Eq. (5), the expression for the latent dynamics is as simple as possible.
4. 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.
5. We will add intermediate ticks to 3F.
6. We will correct these typos.
7. 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}$.
8. It is an inner product. We will clarify this in the Appendix.

Final remarks

In 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.

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.

Answers to Strengths and Weaknesses

The 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: $W_{ij}=J\sum_{l=1}^4 \beta_l \cos(l(\theta_i - \theta_j) - \Delta).$ This 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.

Answers to Questions

Q: What exactly of the model in (13) is fitted to the neural activities?

A: 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.

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?

A: 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}_*.$

Q: L. 957 (appendix). should the second z_s be z_t ?

A: The reviewer is right, we will correct this typo.

Q: Is the factor (2/ \pi) in the equation below L. 921 correct?

A: The factor is correct but there were indeed some factors $1/(2\pi)$ missing in ll. 915-920. We will correct these typos.

Q: (17) => C_{ij} = 4 * [term in (18)] - 4 * (1/4)

A: We will do the recommended modification, which helps readability.

Final remarks

We 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.

The reviewer appreciate the relevance of our results for systems neuroscience and highlight the quality of our text.

Weakness 1

(i) Choice of NCE hyperparameters

The 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.

We 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.

Of 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.

To 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.

Drifting gratings (d = 4):

Natural images (d = 32):

We 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.

(ii) Linear model baseline

The 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.

(iii) Testing NCE on simulated data

We 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).

Weakness 2

Our 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.

Weakness 3

We 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

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:

Question 1

The paper starts with the solvable RNN model for two reasons.

• 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:

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.)

• 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.

Question 2

See answer to Weakness 1(i) above.

Question 3

NCE, 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.

Question 4

In 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.

In 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.

We thank the reviewer for your thoughtful comments, we appreciate your interest in our work. The suggestions have improved the overall quality and impact of the paper. We will take your final thoughts into consideration for any future follow-up work.