Taking the neural sampling code very seriously: A data-driven approach for evaluating generative models of the visual system
Authors: Suhas Shrinivasan
Konstantin-Klemens Lurz
Kelli Restivo
George Denfield
Andreas S. Tolias
Edgar Y. Walker
Fabian H. Sinz
Venue: NeurIPS 2023 poster (2023)
Keywords: Neural Sampling Code
Probabilistic Inference
Bayesian Brain
Macaque V1
Natural Images
Population Recordings
Normalizing Flows
Probabilistic Models
Computational Neuroscience
Theoretical Neuroscience
TL;DR: This paper provides a formalization of neural sampling code that allows one to learn and quantitatively test the generative models of the brain.

Abstract: Prevailing theories of perception hypothesize that the brain implements perception via Bayesian inference in a generative model of the world.
One prominent theory, the Neural Sampling Code (NSC), posits that neuronal responses to a stimulus represent samples from the posterior distribution over latent world state variables that cause the stimulus.
Although theoretically elegant, NSC does not specify the exact form of the generative model or prescribe how to link the theory to recorded neuronal activity.
Previous works assume simple generative models and test their qualitative agreement with neurophysiological data.
Currently, there is no precise alignment of the normative theory with neuronal recordings, especially in response to natural stimuli, and a quantitative, experimental evaluation of models under NSC has been lacking.
Here, we propose a novel formalization of NSC, that (a) allows us to directly fit NSC generative models to recorded neuronal activity in response to natural images, (b) formulate richer and more flexible generative models, and (c) employ standard metrics to quantitatively evaluate different generative models under NSC.
Furthermore, we derive a stimulus-conditioned predictive model of neuronal responses from the trained generative model using our formalization that we compare to neural system identification models.
We demonstrate our approach by fitting and comparing classical- and flexible deep learning-based generative models on population recordings from the macaque primary visual cortex (V1) to natural images, and show that the flexible models outperform classical models in both their generative- and predictive-model performance.
Overall, our work is an important step towards a quantitative evaluation of NSC. 
It provides a framework that lets us \textit{learn} the generative model directly from neuronal population recordings, paving the way for an experimentally-informed understanding of probabilistic computational principles underlying perception and behavior.

Taking the neural sampling code very seriously: A data-driven approach for evaluating generative models of the visual system
Suhas Shrinivasan'.1, Konstantin- Lurz', Kelli Restivo?,
George H. Denfield', Andreas S. Tolias?. Edgar Y. Walker*, Fabian H. !.nz!
Institute for Computer Science and Campus Institute for Data Science, University of Göttingen, Göttingen, Germany
Center for Neuroscience and Artificial Intelligence, Department of Neuroscience, Baylor College of Medicine, Houston, USA
3 Department of Psychiatry, Columbia University, New York City, USA Department of Physiology and Biophysics, and Computational Neuroscience Center, University of Washington, Seattle, USA
5 Department of Electrical and Computer Engineering, Rice University, Houston, USA Equal contribution
"Correspondence: suhas. wewiiiiasannnnnngettttttee
Abstract
Prevailing theories of perception hypothesize that the brain implements perception via Bayesian inference in a generative model of the world. One prominent theory, the Neural Sampling Code (NSC), posits that neuronal responses to a stimulus represent samples from the posterior distribution over latent world state variables that cause the stimulus. Although theoretically elegant, NSC does not specify the exact form of the generative model or prescribe how to link the theory to recorded neuronal activity. Previous works assume simple generative models and test their qualitative agreement with neurophysiological data. Currently, there is no precise alignment of the normative theory with neuronal recordings, especially response to natural stimuli, and a quantitative, experimental evaluation of models under NSC has been lacking. Here, we propose novel formalization NSC, that (a) allows to directly fit NSC generative models to recorded neuronal activity in response to natural images, (b) formulate richer and more flexible generative models, and (c) employ standard metrics to quantitatively evaluate different generative models under NSC. Furthermore, we derive a stimulus-come connullu--ondiioone predictive model neuronal responses from the trained generative model using our formalization that we compare to neural system identification models. We demonstrate our approach by fitting and comparing classical- and flexible deep learning- based generative models on population recordings from the macaque primary visual cortex (V1) to natural images, and show that the flexible models outperform classical models in both their generative- and predictive-model performance. Overall, our work is an important step towards a quantitative evaluation of NSC. It provides framework that lets us learn the generative model directly from neuronal population recordings, paving the way for an ppperimentally-iffome understanding of probabilistic computational principles underlying perception and behavior.
Introduction
Our environment is riddled with sensory stimuli that are noisy, ambiguous, and often incomplete, necessitating organisms to handle uncertainty in their sensory observations. Bayesian models of
37th Conference on Neural Information Processing Systems (NcurIPS 2023).

---

perception and behavior have thus grown in prominence, successfully accounting for an extensive array of tasks across perception 2 cogmition [3]. sensory-motor learning 4 and decision making (588) These models posit that the brain maintains statistical generative model of the world, where sensory observations x are generated from unknown, world-state variable z. Upon encountering a stimulus x, perception in the brain is conceptualized as probabilistically inferring the world-state variable z that caused x. In other words, to perceive x, the brain would invert the generative model to compute the posterior distribution over z: p (z)x). One can express the posterior via Bayes' rule as: p (z|x) X p (x|z) p (z), where p (z) is the prior distribution of z and p (xlz) is the conditional distribution characterizing how well given describes x. While this has been an influential framework, the neuronal underpinnings of probabilistic inference remain challenging to conceptualize and test experimentally. To this end, the Neural Sampling Code (NSC) prominent theory that offers a unique link between neuronal responses and probabilistic inference. Specifically, NSC posits that neuronal responses, r, given stimulus, x, can be thought of as samples drawn from the posterior distribution: r ~ p (z) x) (Figure D
( [6 is
Background and related work Prevailing literature NSC uses simple and restrictive generative models and performs qualitative comparisons of model predictions with neurophysiological data to test the theory. Notably, existing NSC works use simple prior- and conditional distributions with pre-specified parameters. For example, a popular choice for the conditional distribution of images (stimuli) has been Gaussian with a likelihood function that linearly combines pre- pre-specified filters. Hoyer and Hyvärinen [9] learn these filters via independent component analysis on natural images, whereas Haefner, Berkes, and Fiser [12 use oriented Gabor filters instead. Similarly, popular choice for the prior is the exponential distribution with a pre-specified rate parameter (9]. These choices are inspired by (a) what is already known about sensory neurons, especially in the primary visual cortex (V1), and (b) the fact that renders posterior computation mathematically simpler. In the examples above, the choice of filters reflects well-known findings that the receptive fields of V neurons resemble (Gabor-like) orientation filters [1720], and the exponential prior is motivated by the principle of sparse coding 2
20 Importantly, these parameters and distributions and thereby, the generative models are not informed or learned explicitly from neurophysiological data. Rather, these works typically sample from the posterior of the assumed generative model in response to strongly parameterized stimuli (e.g., noisy oriented gratings). The models - and thereby the theory- - are then evaluated based on how well the samples qualitatively capture specific neurophysiological phenomena such as the variinar relationship 19 21 22 task-induced noise correlation structures
12] and contextual modulation in V1 neurons 16
In contrast, recent advances in deep learning-based neural system identification models have set new standards in providing expressive models that can faithfully predict neural population responses to natural stimuli (3/33), and offer experimentally verifiable insights at the single-neuron level (3T 40
42
Additionally, advances in generative modeling, especially of images, have clearly demonstrated
the effectiveness of deep, highly nonlinear, generative models such as auto-regressive models 43
44]. variational autoencoders (45)477), normalizing flows (48-5T), and diffusion models (52 54]
Given the complexity of high -igh-dimenional natural stimuli and real-world tasks, it is paramount that NSC be considered under a generative model that can match such complexity.
Our objective and contributions Here we ask: what exactly is the brain's 's generative model? More specifically, can we identify the brain's 's generative model from NSC population responses natural stimuli? Although simple generative models and qualitative evaluations in the NSC literature have offered us great insight into the potential generative models of the brain and engendered support for the theory, there remains a conspicuous gap in the quantitative evaluation of NSC, particularly in response to natural stimuli. In this work, we bridge this gap by proposing a formalization of NSC that allows us to directly fit NSC generative models to recorded neuronal activity in response to natural images, formulate richer and more flexible generative models, and employ standard metrics such as log- likeliheeooo and single trial correlation to quantitatively evaluate different generative models under NSC. As opposed to specifying generative model that ought to be maintained by the brain, our framework allows us to learn the generative model directly from neurophysiological data. Learning expressive generative models in a data-driven fashion additionally lets us take advantage of population recordings of large and ever-increasing scale in the field 57. Furthermore, our formalization lets us derive a stinuuuu-condittooeed predictive model of neuronal responses from the trained generative model, which can be directly compared to state-of of-toe--tthttt art system identification
2

---

A
B
11
22
(z|x
model
Neuron 2 respionsenn seppsnnsss Fz
r r~p(z|x0)
Z
P
Z
(xlx)
sampie
p(xiz)
(z|x
X
((1
Neuron responses.
Figure I: Conceptualizing NSC. A. Latent variable model of the world (stimulus): z is the world state variable (intensity of oriented Gabor filter here; figure inspired from Orbán et al. (13) and x is the observed sensory stimulus (c.g, an image brain encounters stimulus x, it inverts its generative model, combining the likelihood (x (z) and prior (z), to obtain the posterior (z) x), and are samples from the posterior. C. Neural response distribution under NSC: Each point corresponds to single response pair of two NSC neurons under three distinct stimuli depicted by distinct colors. The distribution of neurons matches the distribution over the corresponding latent variables
tiger). B. Responses r under NSC: As the
models. The predictive model has the ability to provide experimenntally-verffiaai neuron-specific predictions from the normative theory.
We demonstrate our approach by fitting classical generative models from NSC literature and flexible deep learning-basse generative models on macaque primary visual cortex (V1) population responses to natural images. We show that the flexible models outperform classical models in both their generative- and predictive-model performance. Overall, this work presents an important step towards a quantitative evaluation of NSC, paving the way for data-driven approach learning the generative model of the brain.
2 Fitting the Neural Sampling Code
2.1 Theory
An explicit formalization of NSC We begin by formalizing NSC as latent variable probabilistic model - X r. where z represents the world state variable that underlies the observable stimulus x (Figure TA). Subsequently, the stimulus x gives rise to the neuronal responses r via the posterior p (z|x) (Figure BB. NSC posits that the neuronal responses r elicited by stimulus x can be interpreted as stochastic samples from the posterior distribution y (z|x) (Figure IB, C). However, the exact relation between and r is often left unspecified. For instance, it is not clear what aspect neuronal response (e.g., firing rate, presence or absence of spikes, membrane potential) should be treated as a sample. In fact, most previous works do not make a distinction between r and z, and simply equate aspect of the neuronal response such firing rate with the latent sample. Here, make this assumption explicit and treat the neural response r as a random variable that matches the latent random variable in stimulus-conditiones distributions:
ssample p (z
(1)
r Zsample
(2)
Equation is slight abuse of notation state the equivalence in the stimulus-condiiooeed distribu- tions of r and z, more formally stated as:
(r|x) d p (z)x))
where _ denotes equality in distribution or density function. By marginalizing the stimulus, we find that the marginal distribution of must also match that of
p ( / p (x)x) d (
(4)
3

---

Explicitly formalizing NSC with distinct and z has two advantages Firstly, the resulting formulation provides the crucial link between the generative model - x and observed responses serving as the hasis for learning the generative model from the responses. This also provides the basis for a neuron-specific model comparison between different NSC models as well as the possibility to make predictions for specific neurons that can be experimentally tested. Secondly, the explicit link highlights the possibility to explore more flexible mappings between z and r. For instance, one could assume that the latent variable z is encoded in the membrane potential, but what we observe are spike counts, i.e. r = f(z) for some stochastic mapping f. This relation can, in turn, become part of the model, which can then be fitted to real data and compared to alternative versions of the model. In this work, we choose to learn the generative models from the data under the simplest mapping of r = z
(Section 3.2). Please see Section 4 for discussion alternative mappings between and z.
In NSC, we note that the latent variables are what underlie the stimulus (such as the intensity of an oriented filter in [A) and are not necessarily any task-relevant experimenter-defined variables (such as orientation in an iiiiratiionddscciiinatiion task). This is contrast to the alternative theory of probabilistic population codes [58] [5], where typically, the latent variable is explicitly defined to be
the task-relevant experimenter-define variables.
Learning the generative model under NSC One way to quantitatively test an NSC generative model p (z. x) testing how well the response distribution p(r(x) approximates the posterior p (z)x) of the generative model. i..,, testing equation 3 However in reality (z, x), and consequently p (z]x). unknown to us. However, our formalization allows to learn the generative model p (z, x) via learning the joint distribution x). The equivalence between p (z, x) and p (r, x) follows from equations and
P (z, = pz(x) (x) - (r(x) (x) = (r,
(5)
= px(r) (r,,
(6)
The joint distribution can then simply be fitted to recorded stimulus-repponee pairs fx(") (x))) J by maximizing the likelihood:
0" arg max
log
(x
(1)))
arg max
Σ logp ( (0) logp (p()0pp) (7)
0g.0p
Joint
Likelihood
Prior
where 0 are the parameters of the generative model, that we split into 9r and Op for parameters relevant to the likelihood and prior, respectively. Provided that 0r and Op do not overlap, the generative models can be learned by learning the likelihood and prior separately.
Evaluating NSC on data Fitting generative models under NSC on recorded data allows us to compare the generative models quantitatively by evaluating their performances as log- log-likelioood on held-out test set. Furthermore, once we have learned the generative model (r, x; 0*), we can invert it to arrive at the posterior p (r|x; 0*). This provides a neuronal encoding model p (r(x) under specific assumptions of the NSC model, allowing us to predict neural responses to arbitrary new stimuli, The performance of this predictive model can serve as yet another metric for quantitative model comparison under NSC. Additionally, the posterior allows us to compare the generative models to the cormative-theory-oe system identification models. is important to note that our quantification does not make any assumption about the kind of stimuli x. Existing works on NSC use parametric stimuli from classical neuroscience experiments and perform a qualitative comparison between model- and real-neuronal responses to the same stimuli. Our formulation, on the other hand, allows compare different NSC models on natural images, the type of stimuli the visual system has evolved to process.
2.2 Models
Following previous work in NSC (9-16), here we focus on vision and develop generative models under NSC for natural image stimuli x and spike counts r recorded from the visual cortex (Figure A). Developing generative models entails developing models for the prior p (r) and the likelihood ())
(Figure 2 B). Additionally, we also fit an approximate posterior y (r|x). We summarize our fitting
methodology in an algorithm towards the end of the section (Algorithm I

---

A
B
Pouston
depth probe
recoraligs
Salke
Op
2/r
-o
g (..1)
Pr (F
sorthng
Frow tranaform,
Quartizrr Dequrtieen
(xlp,
Graycoale natural Images
Ausika Siattig mosus maciqus
Spakna
Parametertoed
Prior
Likothood
Posterlor
Figure 2: A. Overview of the experimental setup for neuronal recordings (see Section 3.2 for
summary of the data description and Cadena et al. (29] for the complete details). B. Parameterized generative model for NSC: 0p: parameter of the prior p (r; 8p); Oz: parameter of the likelihood (x|z;0z) d: parameters of the (approximate) posterior (r)x; Φ). C. Flexible prior, likelihood, and posterior models. The prior follows our dequantization framework consisting of components: (1) continuous prior p( ((::w) (normalizing flow with Gaussian base distribution), (2) quantizer P(rlÇ) (floor function), and (3) variational dequantizer 9 (C(r: v/ conditional normalizing flow, Appendix B). The likelihood is an isotropic Gaussian distribution over where the parameters are functions (MLP) of r. The posterior is a Gamma distribution over whose parameters are functions of x, modeled as system Dentiifiattion-base convolutional neural-network model (Appendix C
Prior Spike counts r are discrete. Hence, we can neither directly fit standard literature models, such as an exponential ( Laplace distribution 20 the discrete variable r, nor can we straightforwardly fit more common flexible density models such as normalizing flows 44 60 as they are all continuous density models. A common approach remedy this is to employ uniform dequantization, where the discrete quantity is converted into continuous signal by adding uniform random noise [43 61 62]. We adopt the more general approach of variational dequantization where the noise distribution is learned instead of being fixed to be uniform (63 In this method, prior distribution over discrete is captured by positing a generative model involving continuous latent variable Ç linked to the discrete response via a deterministic quantizer function: P ((r)= P(r()))(C)dC where P(r(C) is the quantizer and p(C) is the continuous prior.
Since the integral is usually intractable, the whole model is fit by optimizing the evidence lower bound (ELBO):
Continuous prior
Dequantizer
log P(r) 2 (c~g((ir:))
log p (S:w)
+ H
(CIr;
(8)
where (C|r; u) is the approximate posterior distribution with parameters p (S: is the continuous prior with parameters w. and H (q (C]r;)) -Ec=g(((r::)) [log q (C(r; v)] the conditional entropy
of the dequantizing distribution. Note that Equation (8) only provides lower-bound to log P (r),
and a tighter bound via importance-weiihtee sampling 68] (Appendix A].
In our work, we only consider factorized prior distributions, i.e., we treat neurons to be priori independent log P (r) = -, log P (rx). This choice is informed by both the nature of the neural data that showed limited correlation across all stimuli and the simplicity of the fit it provides. The same independence assumption was applied for the continuous prior distribution over the dequantized responses Ç. Given this dequantizer framework, explore three different NSC priors by varying the distribution over the continuous latent p (C:
1. Exponential (Exp), 1 exp - (S), as found in the original NSC model by Hoyer & Hyväri- nen D
2. Half-normal (HN),
L2 oVT exp
2² 2002 (-
H (S), where H (C) is the heavyside function.
3. Normalizing flow (Flow): p(C; w) Pouse (T (T aT ar-
'((=)) where we choose Plase to be aç
a standard normal, and - represents the following series of invertible mappings with learnalee parameters w: Jaffine, tanh, affine, tanh, affine, tanh, affine, softplus where softplus log(e" - 1), affine ay + b with learnable parameters and b. softplus -1 ensures that the support of Ç is non-negative, since we are ultimately interested in modeling the distribution of (non-negative) spike counts (Figure Cc "Prior" sub-panel).
5

---

For the dequantizer distribution, / (C(F; v),, we utilize conditional normalizing flow-based flexible distribution, as [66] (Appendix
Likelihood
We model the likelihood as an isotropic Gaussian distribution:
p (xxx(()) = =N (xlxx()) 2(6) 1)
(9)
where the parameters mean, (1) E Rl and variance, o (2) E 15 are functions of response, r(1) and |x] is the number of dimensions of x. We consider (1) a linear function where μ war(() + by and exp">r")) (Lin) and (2) a nonlinear function μ Wy MLP((()) + by and a exp the ML where (.) :) is neural network (MLP) (Figure "Likelihood" sub-panel).
Posterior In many cases, the posterior distribution for desired generative model is not analytically tractable and must be approximated, commonly using variational inference or Markov Chain Monte Carlo sampling 69 71 Here, since we learn generative model p (x, r, 0*), we can approximate the true posterior [ x; 0*) by fitting a model posterior q (r|x; Φ) to samples from p p rx)) directly via maximum log- likeleikhooo::
arg Σ
(10)
where x', (x, r"),, are samples from the trained generative model. We model the posterior distribution of responses conditioned on images factorized Gamma distribution, following state-of-the-art (SOTA) work in system identification 72 p (r|x(0) II Pr (ryla("), (())
where ( is the ith image, r, is the jth neuron out of r = S total neurons, and the parameters concentration, (+) and rate, (3() are functions of the image, x(0), Since these functions map an image to response distribution parameters, we model them using a convolutional neural network
model (Figure 2 "Posterior" sub-panel), following SOTA system identification work 225 (Appendix
Prior
Likelihood
Name
Algorithm 1 Learning NSC models from data
Exponential (A = 1)
Linear
Exp1-Lin
Require: N pairs of stimuli and neuronal re-
sponses respectively, (x("), ( N
Exponential
Linear
Exp-Lin
Half-Normal
Linear
HN-Lin
Learning generative model p (x, r; 0p, (p,00,), 0r.)
Normalizing Flow
Linear
Flow-Lin
1: 8% + arg max >N log 0p)
Exponential
MLP
Exp-MLP
07 + arg max,
logp (x() (0)0))(00)))))))))))))
Half-Normal
MLP
HN-MLP
Learning approx posterior model (r|x; Φ)
Normalizing Flow
MLP
Flow-MLP
3: Sample (x'(i), r()) rx
Table 1: Generative models that we fit as be- ing composed of priors and likelihoods.
4: o² + arg max E; log (1100
Experiments
3.1 Synthetic data
We simulated 10,000 pairs of images and neuronal responses from the following three classical NSC models: Hoyer & Hyvärinen model (HNH) with exponential prior 191.
Olshausen
& Field (ONF) model where the prior is a Laplace distribution (20). and full Gaussian model (Gauss) where the prior is an isotropic Gaussian with mean O and variance o². All the three models share common linear, isotropic Gaussian likelihood (x|r) = N (x Ar, o²1). where A is the factor loading matrix learned via standard independent component analysis model (ICA) with a complete sampled image-response pairs from flexible model with Flow prior (described in Section 2.2 and MLP-based likelihood (Section 2.2). where all parameters were randomly initialized. For any given generative model, first sample neuronal responses from the prior via (5) p (r). and then sample corresponding images via
basis natural image patches 9 74 Additionally,

---

x(() (x/x()), where i E I...., 10, 000. We hold out set of 1,000 pairs as the test set. We fitted all models on the datasets simulated from the classical as well as the flexible models via Equation and computed joint lig-liiclihoods of the trained models on the held-out test set as log (x r(*)
logp (x(+) r(()) + log p ). For the classical models, maximum likelihood estimates of the param- eters were obtained analytically (Appendix F.. We trained the flexible model using gradient descent.
We find that (1) the flexible model fits responses and images simulated under other NSC models well, i.e., learns p (r) and p (x(r) and closely approximates the log- likelihood of the true models. Importantly, it outperforms the fit of other NSC models with mismatched generative distributions, consistently being the best model after the ground-trute truth model (first columns in Figure 3 Further- more, the flexible model is capable of generating complex image and response distributions that could not be easily captured by the classical generative models (column 4 in Figure 3) This demonstrates that our framework allows for NSC model fitting and that the flexible model has the ability to flexibly capture the data distribution across wide range of generative models. Critically, a flexible model could fit complex generative models that cannot be modeled well by other classical models.
Models
HNH
ONF
Gauss
Flex
HNH
-92.55
-49.81
Best
ONF
-22,77
-381.58
-360.19
-460991
Gauss
131.58
388.44
335.80
523.58
Worst
NA
Flex
-93.62
3386.00
338.554
217.94
Figure 3: log (x, r) of models on sim- ulated data (trial averaged, in bits). Col- umn denotes the model generating the samples (data) and rows the trained NSC model. Since the exponential prior in HNH has a non-negative support and does not match that of ONF and Gauss, the scores for HNH under ONF and Gauss data are unavailable.
3.2 Neurophysiological data
Data description Next, we demonstrate applying our ap- proach to real neuronal data. We used 32-channel laminar
NeuroNexus arrays (Figure ( to record population activity from the primary visual cortex (V1) of two awake male rhesus macaque monkeys (Macaca mulatta) (29] as they fixated on grayscale natural images sampled from the ImageNet dataset (75.. All the experiments concerned with the recordings adhered to the National Institutes of Health, United States guidelines, and received approval from the Institutional Animal Care and Use Committee. Each image was presented for 120 ms, and spike counts between 40 ms and 160 ms after the image onset were computed and used as the neuronal response r. The image stimulus x used for modeling is 41 x 41 px. For more details on the experiments and data collection, refer to Cadena et al. [29]] We collected data across 12 recording sessions, each having approximately 16,000 image-response pairs and at least 16 well-isolated single units. We split the dataset into approximately 10,000 pairs for training. 3,000 pairs for validation, and 3,000 for testing (for exact details on all sessions, see Appendix DD. We do not aggregate data across sessions and fit models separately for each session since the images can differ from session to session, and not all neurons have seen every image.
Fitting the generative model Given a dataset of images and responses, {x, r) we fit the likelihood p [x|r; 0z) on the umage-respons. pairs and the prior p 0p) on the responses, r as in Equation (7) We fit all of the generative models on each recording session as following procedure described in Table I Below, we describe our results for the session with the largest number of neurons (29 well-isolated single units) in detail and report summary results on all sessions.
For the prior models (Figure 4 A)), we report the test-set log- ligeliidl performance of all models (Exp, HN. Flow) relative to the Empl-model as the baseline. We find that our flexible normalizing flow model (Flow) achieves the best performance, improving the score from the exponential distribution (Exp) by 0.095 bits per neuron per trial, amounting to 2.755 bits across 29 neurons per trial. For the likelihood models (Figure 4B), we find that using an MLP likelihood function, the model improved by 0.052 bits per pixel per trial, amounting to 87.19 bits of improvement across all 41 x41 pixels per trial, relative to the model with linear likelihood function. For the joint distribution (Figure A C). we find that the flexible model (Flow-MLP) achieves the highest log- likelihood score, offering an improvement of 1.8452e-3 bits per pixel per neuron per trial, amounting to 89.951 bits across all 41 x41 pixels and 29 neurons per trial. We observed that in each of the cases, flexible models (Flow prior, MI. likelihood) offer much higher log-likelihood performance, with the same trend found across all sessions (Appendix D)
7

---

A
B
0.0030
Palative log- aop-hhllhhdd (bover bound)
34
-1 86
Ralative log /kalhood (lovr bound)
-1 87
0.0025
1.32
Hoyer & Hyvärinen model
Log Hig-ioo
-1 88
0.0020
30
0.0015
18
-1.90
-199
0.0010
-1.95
124
0,0005
EXP
HN
Fo
Lin
MLP
HN-Lin
Flow-Lin-
EXp-MLP
HN MLP
Fow-MLP
Exp
Prior (r)
Likelihood p (x|r)
Joint p (x, r)
Figure 4: Log-likelihood scores (in bits) of generative models on population recordings (test set) as fit on recording session with the highest number of neurons (n = 29). Error-bars denote the standard error of mean across trials. A. Prior models, (r): log- ligliiilihood (lower bound) relative to the baseline Expl prior model. The score is averaged across neurons and trials. Note that for the prior models discrete spike counts, bound (Equation
only obtain lower bound p (r). Here we show with 1000 samples. B. Likelihood models p rx(r):
the apportance-samplingg absolute log- likelihood of likelihood functions, averaged across image pixels and trials. C. Joint models p (x, r): log- ligllikeihoods relative to the baseline Expl-Lin generative model. The score is averaged across pixels, neurons and trials.
Figure Posterior performance of NSC models, along with system identification model ("Sysident").. Error bars represent the standard eтror of mean over neuronal responses. All metrics are averaged over number of neurons and trials. A. The lower bound of bog-likellhhod bits, for session with 29 neurons. We compute the lower-bound since we are evaluating the Gamma posterior (discrete) spike counts, and full-likelihood is intractable (Appendix G). B. Same as A but across all 12 sessions (purple: average, gray: single session). C. Single-trial correla- tion, for session with 29 neurons. D. Same C across all 12 sessions
-3.5
Log igglikehhodd
Nower bound)
4.0
55
-6
-5.0
-55
-8
Conelation
0.3
0.3
0.2
0.2
Exp1 Lin
Exp-Lin
HM HN--in
Flow-
Exp MLP
HNV-MLD
flow MLO
Sysident
Exp-Lin
HN-Lin
Flow Lin
MLP
HN--MLP
Flow-MuP. Sysident
NSC genarative models
NBC generalve models
Learning the posterior distribution For each of the trained generative models, (xx 8*) we
approximated the model's posterior distribution p (r|x; 8*) using an approximate posterior g (r x; o) trained on samples drawn from the trained generative model (Equation (T)) Algorithm Tb. We
evaluated the posterior distribution for each generative model by computing their mean log-likelihood on real neuronal responses conditioned on real images from the test set (Figure 5 B). We also computed a single-trial correlation between the mean of the learned posterior distributions and neuronal responses (Figure C, D). Finally, compared the posterior distribution to deep system identification model.
We find that, in general, more flexible trained generative model tends to yield a higher posterior predictive performance. Based on the log-likelihood evaluation, our flexible generative model (Flex- MLP) gained as much as 1.39 bits per neuron per trial compared to the baseline Expl-Lin model and 0.61 bits per neuron per trial compared to Exp-Lin model (the Hoyer & Hyvärinen model (9). ). In terms of single-trial correlation performance, our flexible generative model (Flex-MLP) achieved 10% higher correlation compared to the Expl-Lin baseline and 8% higher correlation compared to Exp-Lin. When averaged across all sessions, we find that Flow-Lin performs best, almost

---

par with the Flow-MLP, which achieves 0.019 bits per neuron per trial less. Furthermore, a system identification model trained on the dataset of real neuronal responses performed better than the best NSC generative model, gaining 0.17 bits per neuron per trial and 17% higher correlation per neuron per trial compared to Flex-ML.P.
All model training was performed using backpropagation and gradient descent and we provide training, compute and infrastructure details in Appendix B
4 Discussion
The main focus of this work to develop way answer the question, "How well does NSC explain neurophysiological data quantitativel??.. While NSC is a prominent normative theory for probablistic computation in the brain, and the literature has provided much qualitative insight, our work is the first to offer a quantitative paradigm for empirically testing it using brain responses to ecological, natural stimuli. Our framework additionally lets us formulate more flexible generative models- - which can be better informed by the data - and employ standard metrics such as log- likelihood to quantitatively evaluate alternative generative models under NSC. Furthermore, inverting the learned generative model has allowed us to obtain the posterior distribution, which is equivalently a neuronal response predictive model. Importantly, this let us compare NSC models models outside of NSC's theoretical framework, such as system identification models, allowing us to benchmark the predictive performance of NSC models. Our results demonstrated that the flexible generative models outperformed classical models in terms of both generative and predictive model performance, yet system identification models achieve superior response- response-predictive performance compared even our best generative models. We now discuss some limitations of our current study and discuss a number of open questions and imminent future directions.
Limitations I: Assumption of strict 1:1 neuron-latent mapping: One limitation in our current study is that we only use 1:1 identity mapping between the activity of neurons and latent variables in our formulation of NSC. Abiding by this restriction could limit the capacity of the NSC models, especially considering some existing work in NSC that have qualitatively explored more flexible mappings. For example, Orbán et al. [13] model membrane potential values (responses) as a nonlinear function of posterior samples. Furthermore. Savin and Denève (6) map responses of N neurons to D latent variables where N > D. Many more ways of how r and relate are conceivable. However, our formulation with a separate r and allows us to, in principle, incorporate different mappings and learn the corresponding generative models. Since the focus of this study was on the aspect of fitting NSC models to data, we chose the simplest (original) interpretation of NSC where r = z.
Limitations I:: Definition of a "sample" as total spike counts: We defined a "sample" as the total spike count of neuronal activity within a specific time window following the stimulus, which is not necessarily what literature works do. However, to our knowledge, there is no generally agreed upon or rigorous definition of a "sample" in NSC. While NSC was originally motivated with firing rate/spike counts over a 500ms window as the sample 9 2T). many alternative definitions such as membrane potential over 10ms (13) have been employed. is unclear on what generally applicable metric other than goodness of fit to data- such a definition could be evaluated. This in fact served to us as another motivating factor for striving towards a fata-driven evaluation of sampling models that would allow one to compare such choices in an informed manner. In this work, we chose the total spike count as the working definition.
Limitations III: Better generative models are needed: Advances in deep learning architectures, latent variable models, and transfer learning have greatly enhanced the capabilities of generative models in machine learning. We believe the models we chose, although more expressive than classical models, are still limiting, especially considering that our likelihood p (x/r) uses linear or MLP decoding from neurons to images, with simple Gaussian noise model. To capture the rich and complex nature of neuronal representations of natural images, we believe it is necessary to consider more sophisticated generative models, that even incorporate natural image prior, that would eventually close the gap predictive performance between system identification performance and NSC generative models. Furthermore, an important avenue of research is identifying biological mechanisms that underlie NSC (i.e. sampling from the posterior) [1 77 82] It worth noting that our deep learning- based generative models are not meant to be mechanistic models of NSC neurons.

---

Rather, we believe that our approach lays the foundation for alternative biologically plausible models to be quantitatively evaluated and compared.
Why do system identification models perform better than NSC generative models? System identi- fication models directly trained discriminatively, i.e., ming 1xx [PKL (Puse (r|x) Pnodel (r|x; 0))] to predict neuronal responses to natural images and deep-learning based ones are currently SOTA. There is still much room to build better generative models that would better explain the data (see Limitations III.. However, for given dataset of responses to stimuli from a fixed stimulus distribution, we do not expect the posterior of even the ideal generative model to surpass the performance of the ideal system identification model because the generative model training, i.c., ming [ DKL (Pune (r|x) Pmodel (I x 0)) + DKL (Puue (x) || Pmodel (x: B))}. ;, does not provide any advantage over system identification in response prediction unless some specific inductive biases are introduced in the generative model.
Why bother fitting NSC models if they fail to quantitatively compete with system identification? If we change the stimulus distribution p (x) to Phew (x) with markedly different stimulus statisties and let the sensory neurons adapt to Pnew (x), we would expect the system identification model's s performance to drop on Pnew (x). The system identification model might have to be retrained on new dataset of responses under Pnew (x). This is the case where we would expect the NSC's learned generative model be beneficial. Specifically, change in p (x) to Pnew (x) may entirely derive from the change in prior p (z) to Pnew (z), while p (x (z) remains fixed. Hypothetically. this is since (x/z) represents the invariant "physical" process by which the latents (c.g., the identity of an animal) give rise to observations (e.g.. the appearance of the animal). Consequently. if NSC accurately describes a neural population, i.e, r ~ p (z) x) X p (x(2) p (z). the neuronal adaptation can be accounted for by simply learning Pnew (z), i.e., Fnew ~ Pnew (z|x) X p (x ()) Pnew (z), keeping p (x z) fixed. We believe such out-of-context generalization is a theoretical strength of NSC, and is a consequence of its normative nature (responses being "samples" from the posterior distribution). Such normative hypotheses are neither present in the purely phenomenological system identification models and nor is it straightforward to equip them with normative assumptions.
The above insight thus helps us identify potential future experiments to test NSC models utilizing our framework since it lets us learn the generative model p (x|z) (z) via p (x/r) p (r) (NSC assumption). Namely, one could perform experiments in which we let the neural population adapt to different sensory contexts with expected shifts in (z). Using our NSC framework, we would expect to be able to predict how neuronal responses should change (as reflected in updated Fnew ~ Pnew (zlx)) under new contexts.
Why have previous works not fit NSC models to data? We attribute the lack of such attempts to (1) limitations in data availability, (2) complexities involved in training flexible machine learning and inference algorithms on recorded data, and (3) the philosophical approach behind normative theories. Normative theories describe how biological system ought to function in order to tackle fundamental tasks. They propose models with parameters that are optimized for those tasks, without relying on actual experimental data 83 84]. Typically normative theories are evaluated using qualitative agreements between the proposed models and data. NSC is itself a normative theory. In contrast, phenomenological approaches such as system identification propose models whose parameters are directly learned from experimental data. Normative and phenomenological approaches have historically been developed independently of each other. Similar to Młynarski et al. [83), who interpolate between phenomenological and normative models via maximum entropy priors, our approach allows us to get the best of both worlds: state-of-the-art deep learning-based system identification models from phenomological approaches and the thcoretical underpinnings of the normative NSC. System identification provides with expressive models that faithfully model and predict the activity of thousands of neurons to rich natural stimuli. NSC, on the other hand, goes beyond what experimental data alone could offer by letting us hypothesize how neurons encode uncertainty about the stimulus, reflecting the posterior distribution over latent variables in a generative model of the world, thus allowing us make novel predictions such as about generalizability across stimulus contexts and design relevant experiments.
10

---

Acknowledgements
We thank all the reviewers for their valuable and constructive feedback. We additionally thank Jakob Macke, Xaq Pitkow, Ralf Haefner, Gergö Orbán, members of Sinz-, Walker-, Tolias-lab for helpful and stimulating discussions. SS and FHS are supported by the German Research Foundation (DFG): SFB 1233, Robust Vision: Inference Principles and Neural Mechanisms, TP 06, project number: 276693517. KKL is supported by German Federal Ministry Education and Research through the Tübingen AI Center (FKZ: 011S 18039A), KR and AST are supported by the National Eye Institute, National Institutes Health (NTH), USA with award numbers RO1 EY026927, and Core Grant for Vision Research with grant number 752-EY-002520-37, AST, KR and FHS supported by the National Science Foundation Collaborative Research in Computational Neuroscience, USA with grant number (IS-211311H Germany with FKZ: 01GQ2107. GHD is supported by The National Institute of Mental Health, NIH, USA with grant number T32MH015144 1444. EYW supported by the National Institute of Neurological Disorders and Stroke, NIH, USA with grant number 1 1U19N$107609-01.
References
[1] Yair Weiss, Eero P Simoncelli, and Edward H Adelson. "Motion illusions optimal percepts". In: Nature neuroscience 5.6 (2002), pp. 598-604.
[2] David Knill and Whitman Richards. Perception as Bayesian inference. Cambridge University Press, 1996.
[3] Charles Kemp et al. "A probabilistic model of theory formation". In: Cognition 114.2 (2010),
pp. 165-196.
[4] Konrad P Körding and Daniel M Wolpert. "Bayesian integration in sensorimotor learning". In: Nature 427.6971 (2004), pp. 244-247.
[5] Marc O Ernst and Martin S Banks. "Humans integrate visual and haptic information in statistically optimal fashion". In: Nature 415.6870 (2002), pp. 429-433.
[6] Aaron C Courville, Nathaniel D Daw, and David S Touretzky. "Bayesian theories of condition- ing in a changing world". In: Trends in cognitive sciences 10.7 (2006), pp. 294-300.
[7] Konrad Kording. "Decision theory: what" should" the nervous system do?" In: Science
318.5850 (2007), pp. 606-610.
[8] Nathaniel Daw al. "Cortical substrates for exploratory decisions in humans". In: Nature 441.7095 (2006), pp. 876-879.
[9] Patrik Hoyer and Aapo Hyvärinen. "Interpreting neural response variability as Monte Carlo sampling of the posterior". In: Advances in neural information processing systems 15 (2002).
[10] József Fiser et al. "Statistically optimal perception and learning: from behavior to neural representations".. In: Trends in cognitive sciences 14.3 (2010), pp. 119-130.
[11] Pietro Berkes et "Spontancous cortical activity reveals hallmarks of an optimal internal model of the environment". In: Science 231.6013 (2011), I). pp- 83-87.
[12] Ralf M Haefner, Pietro Berkes, and József Fiser. "Perceptual decision-making as probabilistic inference by neural sampling". In: Neuron 90.3 (2016), pp. 649-660.
[13] Gergö Orbán et al. "Neural variability and sampling-based probabilistic representations in the visual cortex". In: Neuron 92.2 (2016), pp. 530-543.
[14] Rodrigo Echeveste et "Cortical like dynamies recurrent circuits optimized for sampling- based probabilistic inference". In: Nature neuroscience 23.9 (2020), pp- 1138-1149.
[15] Camille Rullán Buxó and Cristina Savin. "A sampling-bard based circuit for optimal decision making". In: Advances in Neural Information Processing Systems 34 (2021), pp. 14163- 14175.
[16] Dylan Festa et al. "Neuronal variability reflects probabilistic inference tuned to natural image statistics". In: Nature communications 12.1 (2021), p. 3635.
[17] David H Hubel and Torsten N Wiesel. "Receptive fields of single neurones in the cat's 's striate cortex". In: The Journal of physiology 148.3 (1959), p. 574.

---

[18] S Marčelja. "Mathematical description of the responses of simple cortical cells". In: JOSA 70.11 (1980), pp. 1297-1300.
[19] John G Daugman. "Uncertainty relation for resolution in space, spatial frequency, and orienta- tion optimized by two-dimensional visual cortical filters". In: JOSA A 2.7 (1985), pp. 1160- 1169.
[20] Bruno A Olshausen and David J Field. "Emergence of simple-cell receptive field properties by learning a sparse code for natural images". In: Nature 381.6583 (1996), pp. 607-609.
[21] AF Dean. "The variability of discharge simple cells the cat striate cortex". In: Experimen- tal Brain Research 44.4 (1981), pp. 437-440.
[22] David Tolhurst, J Anthony Movshon, and Andrew F Dean. "The statistical reliability of signals in single neurons in and monkey visual cortex". In: Vision research 23.8 (1983),
pp. 775-785.
[23]
Daniel LK Yamins et al. Perforraacce-ppiiiizze hierarchical models predict neural responses in higher visual cortex". In: Proceedings of the national academy of sciences 111.23 (2014),
pp. 8619-8624.
[24] Alexander S Ecker et al. "A -uuiiian-quivaaiiant convolutional neural network model of primary visual cortex". In: arXiv preprint arXiv: /809. 10504 (2018).
[25] Santiago A Cadena et al. "Deep convolutional models improve predictions of macaque VI responses natural images". In: PLoS computational biology 15.4 (2019), e1006897.
[26] Benjamin Cowley and Jonathan W Pillow. "High-contrast "gaudy" images improve the train- ing of deep neural network models of visual cortex". In: Advances in Neural Information Processing Systems 33 (2020), pp. 21591-21603.
[27]
Menttantin-Klemenn Lurz et al. "Generalization in data-driven models of primary visual cortex". In: BioRxiv (2020), pp. 2020-10.
[28] Mohammad Bashiri et al. "A flow-based latent state generative model of neural population responses to natural images". In: Advances in Neural Information Processing Systems 34 (2021), pp. 15801-15815.
[29] Santiago A Cadena al. "Diverse task-driven modeling of macaque V4 reveals functional specialization towards semantic tasks". In: bioRxiv (2022). pp. 2022-05.
[30] Konstantin F Willeke et al. "Deep learning-drive driven characterization of single cell tuning in primate visual area V4 unveils topological organization", In: bioRxiv (2023), pp. 2023-05.
[31] Pawel A Pierzchlewicz et al. "Energy Guided Diffusion for Generating Neurally Exciting Images". In: bioRxiv (2023), pp. 2023-05.
[32] Charles F Cadieu et al. "Deep neural networks rival the representation of primate IT cortex for core visual object recognition". In: PLOS computational biology 10.12 (2014), 0039633.
[33] Elcanor Batty et al. "Multilayer recurrent network models of primate retinal ganglion cell responses". In: International Conference on Learning Representationss 2016.
[34] Ján Antolık et al. "Model constrained by visual hierarchy improves prediction of neural responses to natural scenes", In: PLoS computational biology 12.6 (2016), e1004927.
[35] Lane McIntosh al. "Deep learning models of the retinal response to natural scenes". In: Advances in neural information processing systems 29 (2016).
[36] William F Kindel, Elijah D Christensen, and Joel Zylberberg. "Using deep learning to reveal the neural code for images in primary visual cortex". In: arXiv preprint arXiv:1706.0623 (2017).
[37] David Klindt et al. "Neural system identification for large populations separating "what" and "where""" In: Advances in Neural Information Processing Systems 30 (2017).
[38] Martin Schrimpf et al. "Brain-score: Which artificial neural network for object recognition is
most brain-like?" In: BioRxiv (2018), p. 407007.
[39] Fabian Sinz et al. "Stimulus domain transfer in recurrent models for large scale cortical population prediction on vidco". In: Advances in neural information processing systems 31 (2018).
12

---

[40] Pouya Bashivan, Kohitij Kar, and James J DiCarlo. "Neural population control via deep image synthesis". In: Science 364.6439 (2019), eaav9436.
[41] Carlos Ponce al. "Evolving images for visual neurons using a deep generative network
reveals coding principles and neuronal preferences". In: Cell 177.4 (2019), pp. 999-1009.
[42] Edgar Y Walker et al. "Inception loops discover what excites neurons most using deep predictive models". In: Nature neuroscience 22. 12 (2019), pp. 2060-2065.
[43] Benigno Uria, lain Murray, and Hugo Larochelle. "RNADE: The real-valued neural autoregres- In: Advances in Neural Information Processing Systems 26 (2013). [44] Aaron Van den Oord al. "Conditional image generation with pixelcnn decoders". In: Ad- vances in neural information processing systems 29 (2016).
sive density
[45] Diederik P Kingma and Max Welling. "Auto-encoding variational bayes". In: arXiv preprint arXiv:1312.611 (2013).
[46] Danilo Jimenez Rezende, Shakir Mohamed, and Daan Wierstra. "Stochastic backpropagation and approximate inference in deep generative models", In: International conference on machine learning. PMLR. 2014, pp.
[47] Rewon Child. "Very deep vaes generalize autoregressive models and can outperform them on images". In: arXiv preprint arXiv:2011.10650 (2020).
[48] Laurent Dinh, David Krueger, and Yoshua Bengio. "Nice: Non- Mon-linear independent components estimation". In: arXiv preprint arXiv:1410.8516 (2014).
[49] Laurent Dinh, Jascha Sohl-Dickstein, and Samy Bengio. "Density estimation using real nvp". (2016).
In: arXiv preprini 160Xvv6605..80
[50] Mathieu Germain et al. "Made: Masked autoencoder for distribution estimation". In: Interna-
tional conference on machine learning. PMLR. 2015, pp. 881-889.
[51] Durk P Kingma and Prafulla Dhariwal. "Glow: Generative flow with invertible 1x1 convolu- tions". In: Advances in neural information processing systems 31 (2018).
[52] Jonathan Ho, Ajay Jain, and Pieter Abbeel. "Denoising diffusion probabilistic models". In: Advances in neural information processing systems 33 (2020), pp. 6840-6851.
[53] Jiaming Song, Chenlin Meng, and Stefano Ermon. "Denoising diffusion implicit models". In: arXiv preprint arXiv:2010.02502 (2020).
[54] Robin Rombach et al. "High-resolution image synthesis with latent diffusion models". In: Proceedings of the IEEE/CVF conference on computer vision and pattern recognition. 2022,
pp. 10684-10695.
[55] Ian H Stevenson and Konrad P Kording. "How advances in neural recording affect data analysis". In: Nature neuroscience 14.2 (201 pp. 139-142.
[56] Peiran Gao and Surya Ganguli. "On simplicity and complexity in the brave new world of large-scale neuroscience". In: Current opinion in neurobiology 32 (2015), pp. 148-155.
[57] Cole Hurwitz et al. "Building population models for large-scale neural recordings: Opportuni- ties and pitfalls". In: Current opinion in neurobiology 70 (2021), pp. 64-73.
[58] Wei Ji Ma et al. "Bayesian inference with probabilistic population codes". In: Nature neuro- science 9.11 (2006), pp. 1432-1438.
[59] Jeffrey M Beck et al. "Probabilistic population codes for Bayesian decision making". In: Neuron 60.6 (2008), pp. 1142-1152.
[60] Danilo Rezende and Shakir Mohamed. "Variational inference with normalizing flows". In: International conference on machine learning. PMLR. 2015, pp. 1530-1538.
[61] L Theis, A van den Oord, and M Bethge. "A note on the evaluation of generative models". In: International Conference on Learning Representations (ICLR 2016). 2016, pp. 1-10.
[62] Aaron Van den Oord and Benjamin Schrauwen. "Factoring variations in natural images with deep gaussian mixture models". In: Advances in neural information processing systems 27 (2014).
13

---

[63] Jonathan Ho et al. "Flow+:: Improving flow-based generative models with variational dequan- tization and architecture design". In: International Conference on Machine Learning. PMLR. 2019, pp. 2722-2730.
[64] Christina Winkler et al. "Learning likelihoods with conditional normalizing flows". In: arXiv preprint 19 10191120022222 (2019).
[65] Emiel Hoogeboom et al. "Integer discrete flows and lossless compression". In: Advances in Neural Information Processing Systems 32 (2019).
[66] Emiel Hoogeboom, Taco Cohen, and Jakub M Tomezak. "Learning discrete distributions by dequantization". In: arXiv preprini arXiv:2001. /2355 (2020).
[67] Yuri Burda, Roger Grosse, and Ruslan Salakhutdinov. "Importance weighted autoencoders". In: arXiv preprini arXiv: 1509.005 (2015).
[68] Justin Domke and Daniel R Sheldon. "Importance weighting and variational inference". In: Advances in neural information processing systems 31 (2018).
[69] David Blei, Alp Kucukelbir, and Jon McAuliffe. "Variational inference: A review for statisticians". . In: Journal of the American statistical Association 112.518 (2017), pp. 859-877.
[70] Christopher M Bishop and Nasser M Nasrabadi. Pattern recognition and machine learning.
Vol. 4. 4. Springer. 2006.
[71] David JC MacKay. Information theory, inference and learning algorithms. Cambridge univer- sity press, 2003.
[72] Konstantin-Klemens Lurz et al. "Bayesian Oracle for bounding information gain in neural encoding models". In: The Eleventh International Conference on Learning Representations. 2023. URL: pettps://openreviee
[73] Luca Baroni al. "Learning invariance manifolds of visual sensory neurons". In: NeurIPS Workshop on Symmetry and Geometry in Neural Representations. PMLR. 2023, pp. 301-326.
[74] Aapo Hyvärinen. "Independent component analysis: recent advances". In: Philosophical Trans- actions of the Royal Society A: Mathematical, Physical and Engineering Sciences 371.1984 (2013), p. 201 10534.
[75] Jia Deng al. "Imagenet: large-scale hierarchical image database". In: 2009 IEEE confer- ence on computer vision and pattern recognition. leee. 2009, pp. 248-255.
[76] Cristina Savin and Sophie Denève. temporl Representations of Uncertainty in Spik- ing Neural Networks", In: Advances in Neural Information Processing Systems. Ed. by Ghahramani et al. Vol. 27. Curran Associates, Inc., 2014, pp. 2024-2032.
[77] Lars Buesing et al. "Neural dynamics as sampling: a model for stochastic computation in re- current networks of spiking neurons". In: PLoS computational biology 7.11 (201 1), e1002211.
[78] Yanping Huang and Rajesh P Rao. "Neurons Monte Carlo Samplers: Bayesian Inference and Learning in Spiking Networks". In: Advances in neural information processing systems 27 (2014).
[79] Guillaume Hennequin, Laurence Aitchison, and Máté Lengyel. "Fast sampling-basse based inference in balanced neuronal networks". In: Advances neural information processing systems 27 (2014).
[80] Xingsi Dong et al. "Adaptation Accelerating Sampling-basee Bayesian Inference Attrac- tor Neural Networks". In: Advances in Neural Information Processing Systems 35 (2022), pp. 21534-21547.
[81] Paul Masset et al. "Natural gradient enables fast sampling in spiking neural networks". In: Advances in Neural Information Processing Systems 35 (2022), pp. 22018-201 822344.
[82] Shirui Chen et al. "Expressive probabilistic sampling in recurrent neural networks". In: arXiv preprint arXiv:2308.. 809 (2023).
[83] Wiktor Młynarski et al. "Statistical analysis and optimality of neural systems". en. In: Neuron 109.7 (Apr. 2021), 1227-1241.e5.
[84] Daniel Levenstein et al. "On the role of theory and modeling in neuroscience". In: Journal of Neuroscience 43.7 (2023), pp. 1088-1088.
14

---

[85] Adam Paszke et al. "Pytorch: An imperative style, bigh-performance deep learning library". In: Advances in neural information processing systems 32 (2019).
[86] Diederik P Kingma and Jimmy Ba. "Adam: A method for stochastic optimization", In: arXiv preprint arXiv: 14 62.69980 (2014).
15


Reviews:
[{"review_id": "zWwBR4gcJR", "reviewer": "Reviewer_AuQ9", "summary": "This paper proposes a statistical framework to infer the generative models hold in the brain from population activity. This study first analytically proves that it is feasible to infer the generative model from population activity. Then, the study compares various combinations of prior and likelihood to fit synthetic and empirical data. They found that the flow-MLP model outperformed other models. In particular, the flow-MLP model allows to flexibly capture the true generative models.", "questions": "1. The logic of neural sampling codes is very clear. I am not doing this line of research but I have know  the difference between neural sampling codes (NSC) and probablistic population codes (PPC) for a long time. More insights are needed to compare NSC and PPC. Is it also feasible to infer generative models under the framework of PPC. More discussions are needed.\n2. I am really confused by several mathematical notations here.\n   * In Eq1, I understand the meaning of r, x, z. but the footnote $r|x$ and $z|x$ are unclear to me. Given this is the key equation in this paper. \n   * Eq2, \"This subtle and seemingly trivial match implies that there exists a one-to-one mapping between the  space of neuronal responses and the latent variable underlying the stimulus\". This sentence is very confusing to me.\n   * Lines85-87, \"However, intriguingly, our formalization in equation 1 allows us to directly learn the generative model pz,x (z, x) by instead learning the probabilistic mapping between the neural responses and the stimuli\". I think this is they key logic. However, given the mathematical notations in Eq1&2 are confusing. This logic is really hard to follow.\n   * In Eqs.6&7, what are the $w$ and $v$? I assume they are parameters of the continous latent distribution? Better to clarify.\n3. In what situations is this approach useful. Can the authors give some cases where experimentalists can apply this approach to their data, or can design a better experiment to test relevant theories.", "limitations": "This study overall is well done. Some notations should be better explained and more discussions are needed in terms of relevance to empirical experiments.", "rating": 6, "confidence": 4, "soundness": 4, "presentation": 2, "contribution": 4, "strengths": "Overall I am positive about this study. This study is well motivated, interesting, and technically sound. The discussion part is also inspiring. \n\n1. Neural sampling codes (NSC) is an important theory of population activity in the brain. This study definitely provides a practical framework that neural sampling codes can be useful in analyzing high-dimensional neural data.\n2. Using population activity to infer generative models is an excellent idea. This study is well motivated.\n3. The experiment is complete and neat. The study includes three priors x two likelihood functions, as well as a system identification model for comparison.", "weaknesses": "1. The writing can be improved. The overall texts are clear but some mathematical notations lack explicit explanations.\n2. The normalizing flow model needs more detail.\n3. More explanations are needed about in what situation this approach is useful in future empricial experiments. Or how experimentalistscan  design and perform better experiments to promote this line of research.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank you for your positive and supportive review. We are positive that we can address all your comments and questions below.\n\n**Re clarification of notation and meaning behind Eq 1, 2 and 3 (key logic)**: In the manuscript, we let the subscript of the probability distribution designate distinct functions (probability density function for continuous and probability mass function for discrete random variables). Namely, $p_{r|x}$: probability distribution of $r$ conditioned on $x$, and $p_{z|x}$: probability distribution of $z$ conditioned on $x$. In Eq 1, $p_{r|x}(r|x)$ denotes the probability density function $p_{r|x}$ as evaluated on values of $r$ and $x$, and, $p_{z|x}(r|x)$ denotes the probability density function $p_{z|x}$ as evaluated on the values $r$ and $x$.  \nIn its simplest form, NSC posits that _a neuronal response to a stimulus is a \"sample\" from the posterior distribution of latent variables given the stimulus_, which we mathematically write as: $ z_{sample} \\sim p_{z|x}(z|x)$ and $r = I(z_{sample})$, where $I$ is the identity function.  What this implies is that there exists a one-to-one correspondence between the $r$ and $z$. In other words, $r$ and $z$ are two random variables that share the same (1) stimulus-conditioned distribution (Eq 1) and (2) marginal distribution (Eq 2). Consequently, the equivalence allows for the generative model of the stimulus $p_{x, z}(x, z) = p_{x|z}(x|z) p_z(z)$ to be learned directly using stimuli $x$ and responses $r$, where $r$ can be used in place of $z$ (Eq 3).\n\n**Re details on normalizing flow**: For modeling the prior distribution over spike counts $P(r)$, we used a normalizing flow distribution to model the continuous prior $p(\\zeta)$ over the dequantized neural responses $\\zeta$ (see Equation 5).\nWe assumed that prior distributions over each neuron response are independent and thus can be factorized: $\\log P(r) = \\sum_j^n \\log P(r_j)$, where $j$ is the neuron index over $n$ neurons.\nThe same independence assumption was applied for the continuous prior for the dequantized neuronal responses  $\\zeta$. \nWe will now describe the model for a single dimension $\\zeta_j$, and for simplicity, we drop the subscript $j$.\n\nWe rephrase our normalizing flow distribution in two steps:    \n$p({\\zeta};\\omega)=p_{\\text{base}}(T^{-1}({\\zeta};\\omega))\\cdot|\\frac{\\partial T^{-1}({\\zeta};\\omega)}{\\partial {\\zeta}}| = \\mathcal{N}(T^{-1}({\\zeta}; \\omega)|0,1)\\cdot |\\frac{\\partial T^{-1}({\\zeta}; \\omega)}{\\partial {\\zeta}}|$\nwhere we choose $p_{\\text{base}}$ to be $\\mathcal{N}(\\cdot|0, 1)$, a standard normal distribution, and $T^{-1}$ represents the following (ordered) series of invertible mappings with learnable parameters $\\omega$: $[\\text{affine}, \\tanh, \\text{affine}, \\tanh, \\text{affine}, \\tanh, \\text{affine}, \\text{softplus}^{-1}]$, where $\\text{softplus}^{-1}(y) = \\log(e^y - 1)$, $\\text{affine}(y) = ay + b$ with learnable parameters $a$ and $b$. $\\text{softplus}^{-1}$ ensures that the support of ${\\zeta}$ is non-negative, reflecting (non-negative) spike counts.\n\nWe identified two errors in our normalizing flow description in the manuscript: (1) we wrote the $\\text{softplus}^{-1}$ layer at the wrong end of the transformation (Line 144) and (2) we used $\\xi$ as the parameters of the flow whereas it should have been $\\omega$ (Lines 125, 144). We will correct them both in the main text.\n\n**Re symbols in Eqs 6 & 7**: We thank you for bringing this to our attention. We missed the description of these two parameters, and will add them to our main text. $\\omega$ refers to the parameters of the continuous latent distribution modeled by the normalizing flow $p(\\zeta; \\omega)$, and $\\nu$ refers to the parameters of the dequantization distribution $q(\\zeta|r; \\nu)$.\n\n**Re feasibility of inferring generative models under PPC**: We thank you for prompting a discussion on PPC and we will add the following discussion to our main text. In its original formulation (Ma, et al 2006; Beck et al 2008), PPC focuses on capturing the likelihood function $\\mathcal{L}(x) \\equiv p(r|x)$, relating the population response to some specific experimenter-defined quantities such as orientation. Hence, fitting PPC models to neural responses amounts to capturing only this likelihood function and contains no notion of the prior over the stimulus $p(x)$. We hence believe that PPC does not readily offer a means for fitting a generative model, complete with the prior. That being said, there exist extensions to PPC where the responses are seen as encoding a posterior distribution of an experimenter-defined variable (e.g. orientation) via a specific distributional family (i.e. exponential family) (Shivkumar et al. 2018). However, in such work, the models are typically limited in that it only captures the relationship between the population response and some experimenter-defined variables such as stimulus orientation. This is in contrast to NSC in which latent variables are seen as what underlies the (natural) stimulus but not necessarily any experimenter-defined variable. To our knowledge, there currently exists no method to identify/learn the target latent variable z and the generative model over (z, x) based on the observed population responses r within the framework of PPC. While this is certainly an interesting avenue of research, deriving such a method would require dedicated work, perhaps similar in extent to our current work on the NSC. \n\n**Re relevance to experiment**: Please see overall response (“**Re 9, 10**”)."}, {"review_id": "2O3A1uyL1D", "reviewer": "Reviewer_1S3S", "summary": "The authors employed a generative model approach to predict V1 neural responses, aiming to evaluate the Neural Sampling Code (NSC) using a quantitative metric. The primary potential contribution lies in the theoretical formulation and the idea that neural response prediction can be utilized to assess the NSC hypothesis, especially on a per-trial basis.\n\n", "questions": "\nHow does the current work speak for NSC hypothesis? Does the result in this paper confirm or refute it? \n", "limitations": "No discussion of societal impact.", "rating": 6, "confidence": 4, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "The authors employed a generative model approach to predict V1 neural responses. This is a novel contribution in the area of single neuron response prediction, even though generative models like VAE, beta-VAE, and hierarchical VAE have been commonly used in analyzing fMRI data, particularly for reconstructing the input image. \n\nThe primary potential contribution lies in the theoretical formulation and the idea that neural response prediction can be utilized to assess the NSC hypothesis, especially on a per-trial basis.  The proposal that the performance of a  generative model for predicting neural response can provide a quantitative metric for evaluating the Neural Sampling Code (NSC) is new and interesting, even though I have not been completely convinced by this work presented (see below).  Nevertheless,  I put this paper as a \"marginal accept\" for the idea and the theoretical contribution.  \n\nNovel and substantial macaque V1 data.", "weaknesses": "The authors employed a generative model approach to predict V1 neural responses, aiming to evaluate the Neural Sampling Code (NSC) using a quantitative metric.   However, the response prediction performance, as measured by correlation, was  poor compared to the standard state-of-the-art feedforward neural network system identification model. Hence, the virtue and the effectiveness of utilizing a generative model in this context remain unclear.  \n\nMore importantly,  it is unclear how this work confirms or refutes the NSC hypothesis. Moreover,  the earlier conceptual framework and theoretical formulation of the NSC propose that samples of population activities within a 10-20 msec window are used to approximate the posterior via Markov chain Monte Carlo (MCMC). In contrast, this work uses responses within a 120 msec duration throughout the entire stimulus presentation as a sample. Consequently, I find it hard to relate the current results to the original NSC idea. Additional arguments and explanations are necessary to establish the connections more convincingly. \n\n", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank you for your encouraging review.\n\nWe are positive that we can address all your questions below. \n\n**Re the correlation performance gap between the posterior and system identification**: Please see our overall response (specifically “**Re 9, 10**”)\n\n**Re confirmation or refutation of NSC**: Our work can be used to (1) learn generative models directly from neurophysiological population recordings and natural images and to (2) subsequently quantitatively compare different generative models proposed for NSC. These comparisons can then form a solid basis to compare different normative assumptions, realized via specific choices in the generative model specification such as sparsity realized via exponential distribution, but all under NSC. Hence, our work in itself neither refutes nor confirms NSC. Rather, (1) we envision future work to use our framework for quantitatively comparing NSC models to models derived from other theories, but most crucially (2) if NSC were true, i.e responses are samples from the posterior distribution $r \\sim p(z|x) \\propto p(x|z)p(z)$, then changes in $p(z)$ would be reflected in the responses. We believe that designing experiments that aim to test this idea would be a strong way to arrive at confirmation or refutation of NSC, and our framework that lets us learn the generative model would serve as the necessary prerequisite. For a more detailed discussion, please refer to our overall response (**Re 9, 10**).\n\n**Re choosing “sample” as total spike count over 120 ms different from literature**: We would like to point out that, as far as we know, there is no generally agreed upon or a rigorous definition of a “sample” in NSC. While NSC was originally motivated with firing rate/spike count over a 500ms window as the sample (Hoyer and Hyvärinen, 2002 referring to A. F. Dean. 1981), some works even employed membrane potential over 10ms (Orban et al. 2016) as the definition of a sample. It is unclear on what generally applicable metric such a definition could be evaluated, _other than goodness of fit to data_. In this work, we chose the total spike count as the working definition. However, identifying the right sample definition falls within the scope of modeling and fitting NSC models to the data. Now that we can quantitatively fit neuronal responses under NSC, we could evaluate different definitions of samples to test which leads to the highest likelihood under  NSC. To us, this was another motivating factor for striving towards a data-driven evaluation of sampling models that allow us to compare these choices in an informed and quantitative way.\n"}, {"review_id": "8h7hgbzHiI", "reviewer": "Reviewer_D2be", "summary": "The authors propose a novel formalization of the Neural Sampling Code (NSC) in this paper. The authors aim to bridge the gap in existing NSC research that needs more quantitative alignment between normative theory and neuronal recordings. \n\nThe authors propose formalization enables fitting NSC generative models directly to neuronal activity, creating richer, more flexible models and allowing standard metrics to evaluate different generative models.\n\nThe authors propose to compare classical and deep learning-based generative models on population recordings in V1 in response to natural images. Their framework outperforms classical models in both generative- and predictive-model performance. They also propose a method for obtaining stimulus-conditioned predictive models of neuronal responses that can be compared to existing neural system identification models.\n", "questions": "How does the model perform in other brain regions? ", "limitations": "Clearly laid out by the authors", "rating": 6, "confidence": 4, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "The authors propose a novel formalization of the Neural Sampling Code (NSC), allowing for the direct fitting of NSC generative models to recorded neuronal activity. Also, the paper presents a new way to evaluate different generative models under NSC using standard metrics quantitatively. \n\nThe study compares classical and deep learning-based generative models, providing valuable insights into the capabilities of different generative models in capturing the structure of real neuronal responses.\n\nThe research is grounded in the practical application of the theory, using neuronal activity data from a macaque primary visual cortex (V1) in response to natural images. \n\nFinally, the authors spend a significant amount of space identifying the limitations and suggesting areas for future research, promoting ongoing dialogue and exploration in the field. \n", "weaknesses": "I believe that the main weaknesses of the current paper have been presented by the authors in the discussion section but do not undermine the validity of the current study. \n\nHowever, one important limitation lies in the experimental section. The paper provides limiting experimental results when it comes to different neural recordings. In order to clearly evaluate the validity of the claim, one would need to explore other brain areas, probably still in the early sensory processing regions. ", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank you very much for your supportive and motivating review. We hope that our work promotes the ongoing dialogue and exploration in the field.\n\n**Regarding different brain regions**: We share your interest and curiosity in exploring and testing our approach comprehensively in different brain regions and also across different animal species. Given that we have our current framework in place, we are excited that it gives us the opportunity to do so in an upcoming study. \nHowever, the focus of our current work was to establish a novel formulation of NSC that facilitates fitting them to recorded neurophysiological data. Testing it on different brain regions, which likely have no established normative models yet, was out of the scope for the current study. We chose the macaque primary visual cortex (area V1) because most of existing NSC literature has focused on modeling the V1 population responses to visual stimuli. We would like to emphasize that our approach is indeed not limited to V1, but serves as a broader framework that can be used to develop and fit models in other brain areas. \n\nThe consistency of results across all tested datasets gives us confidence in the reliability of our approach and its potential for generalizability of findings in other brain regions.\n"}, {"review_id": "FiuIoKAg0l", "reviewer": "Reviewer_GaXp", "summary": "The authors used several generative models to fit the image data and the corresponding neural responses. And for each generative model, they used the variational inference method to get their posterior.", "questions": "1. What do the two conditional probabilities in Eq.(1) mean? Could you provide a detailed explanation?\n When talking on NSC, I believe the latent variables z should be abandoned. Please see [1].\n\n2. What does the term normative theory/assumption mean?\n\n3. What insights can we gain from the the generative model fitted by neuronal population recordings?\n\nReference:\n[1] Shivkumar, Sabyasachi, et al. \"A probabilistic population code based on neural samples.\" Advances in neural information processing systems 31 (2018).", "limitations": "None", "rating": 2, "confidence": 4, "soundness": 2, "presentation": 3, "contribution": 1, "strengths": "The article is written in a clear manner.", "weaknesses": "The authors used several generative models to fit the image data and the corresponding neural responses. Neither the generative model nor the fitting method is biologically plausible. This cannot explain the working mechanisms of the neural system; it can only assess how well the model fits the data. And it is evident that more complex models tend to fit better. I fail to see the significance of this work.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank you for your critical feedback. We would like to address all of your points that you raise.\n\n**We would like to point out a misunderstanding that appears in your summary**. We do _not_ use variational inference in order to get the posterior, but rather we fit a model-posterior to the samples of the learned generative model (please see paragraph under “Posterior” in Section 2.2 “Models'' and Equation 9). We did so to ensure that we are not biased by the choice of approximation (i.e., ELBO). We realize this is a bit unconventional and can be misleading if not emphasized. We created an algorithm that summarizes our fitting methodology concisely and we state the posterior-learning part for your reference:\n\nSample $\\\\{x'^{(i)}, r'^{(i)}\\\\}\\_i^S \\\\sim p\\_{x, z}(x, r; \\\\theta_P^*, \\\\theta_L^*)$\\\n$\\phi^* \\gets argmax_\\phi \\sum_{i}^S \\log q_{z|x}(r'^{(i)}|x'^{(i)}; \\phi)$,\\\nwhere $p\\_{x, z}(x, r; \\\\theta_P^*, \\\\theta_L^*)$ is the learned generative model.\n\n\n**Re the meaning of the conditional probabilities in Eq 1**: In the manuscript, we let the subscript of the probability distribution designate distinct functions (probability density function for continuous and probability mass function for discrete random variables). Namely, $p_{r|x}$: probability distribution of $r$ conditioned on $x$, and $p_{z|x}$: probability distribution of $z$ conditioned on $x$. In Eq 1, $p_{r|x}(r|x)$ denotes the probability density function $p_{r|x}$ as evaluated on values of $r$ and $x$, and, $p_{z|x}(r|x)$ denotes the probability density function $p_{z|x}$ as evaluated on the values $r$ and $x$. Please also see our response to Rev **AuQ9**  for a reiterated version of our theoretical formulation.\n\n**Re latent variable $z$ in NSC**: We see your point about abandoning $z$ and replacing it with $r$ for NSC as is often done in the literature. We would like to highlight that such an approach entails an _implicit_ canonical assumption about the relationship between $r$ and $z$, specifically that there exists an equality in distribution between $r$ and $z$ (please also see our point under “Limitations I” under the “Discussion” section in the submitted manuscript). We deliberately use $z$ and $r$ as separate variables for two reasons: Firstly, our goal was to mathematically spell out this implicit equality assumption in NSC, and it was hence necessary to start with the idea that there exists a generative model of the world, and link neuronal responses to the latent variables in the generative model. Secondly, making the link explicit also allows one to formulate richer links between $r$ and $z$, which may be seen as extending NSC. For example, one could formulate that neuronal responses represent a subset of latent variables, involving nonlinear mappings. Ultimately, we believe that neuronal responses are a _reflection_ of the posterior samples. Being able to make these assumptions explicit and having a quantitative metric such as data log-likelihood performance allows one to see what assumptions describe the data best and aid follow-up experiments.\n\n**Re what normative theory/assumptions mean**: Normative theories aim to describe how a system _ought_ to function based on first principles, and involve assumptions with biophysical constraints (please see references [49] and [45] for a detailed exposition of normative theories). These principles often come in the form of optimality principles. A prominent example in visual neuroscience is sparse coding by Olshausen and Field (2004), where the underlying optimality principle is the maximization of energy efficiency and robustness to noise, among others. NSC itself is also a normative theory, postulating that the brain’s responses are samples from the posterior distribution over latent variables reflecting a generative model of the world and deriving neuronal coding properties from that. Commonly in NSC models, assumptions such as the equality of distributions between $r$ and $z$ (see point above) and the response/latent variables being sparse (hence modeled by an exponential distribution) are normative assumptions. In sharp contrast to system identification where the parameters of models are learned purely from data, normative theories often propose models that are either parameter-free or have parameters that are independent of data. Furthermore, normative theories are derived from broader computational frameworks, such as Bayesian inference, and have the ability to make predictions that are qualitatively different from system identification (such as how neuronal responses could encode uncertainty about the stimulus and adapt to changing stimulus distributions; please see our point in overall response, specifically “**Re 9, 10**”). Our work can be seen as an interpolation between system identification and normative modeling with the goal to make normative theory of NSC better testable on neurophysiological data.\n\n**Re biological implausibility of the fitting and models**: It was not our aim to specifically provide a biologically plausible model or algorithm. In the tradition of David Marr’s three levels, many prominent normative theories, such as the work by Olshausen and Field (2004), target the computational or representational level of what the computational goal of a system is and how stimuli are represented, but not how it is implemented in a neuronal system (implementational level). Our work targets models for these levels. Asking how a neuronal system would implement such a model and how it would learn it, is an interesting and scientifically important question, but adds a lot more complexity which would require a dedicated treatment. For that reason, we did not focus on that.\n\n**Re what insights can be gained by fitting generative models to data**: Please refer to our overall response (specifically “**Re 9, 10**”).\n"}]

Author rebuttal:
We thank the reviewers for the valuable feedback on our manuscript. We are happy that they found our work to be “well-motivated, interesting and technically sound” (Rev **AuQ9**), with “novel contribution in single response prediction” (Rev **1S3S**), “grounded in practical application of the theory” (Rev **D2be**) and providing “a practical framework that NSC can be useful in analyzing high-dimensional neural data” (Rev **AuQ9**). They remarked on our experiments as “complete and neat” (Rev **AuQ9**) and that we test on “substantial macaque V1 data” (Rev **1S3S**), and found our discussion to be “inspiring” (Rev **AuQ9**).

The reviewers also raised several points to be addressed:

**1**. To clarify our notation, some probabilistic objects (Rev **AuQ9**, Rev **GaXp**) and the logic behind our formulation (Rev **AuQ9**).  
**2**. To add details for the normalizing flow (Rev **AuQ9**).  
**3**. What normative theory/assumption means (Rev **GaXp**).  
**4**. To abandon latent variable $z$ (Rev **GaXp**).  
**5**. On biological implausibility of our methodology (Rev **GaXp**).  
**6**. To discuss our “sample” definition as spike count over 120ms different from literature (Rev **1S3S**).  
**7**. On feasibility of inferring generative models under PPC and make connections (Rev **AuQ9**).  
**8**. How models would perform in other brain regions (Rev **D2be**).  
**9**. What the benefit is from learning a generative model on responses and stimuli (Rev **GaXp**), especially when the posterior correlation is lower compared to system identification (Rev **1S3S**).  
**10**. What new experiments could be designed from using this approach on data, in order to test relevant theories (Rev **AuQ9**), and  
**11**. Whether our study confirms or refutes NSC (Rev **1S3S**).   

We are confident that we can address all of the points in our responses. In this overall response, we first summarize the **main intended contribution** of the paper, since it is relevant for all the questions. Then we  respond to points **9** and **10** since we believe it carries significance relevant to all reviewers. Other points (**1 - 8, 11**) are answered in the responses to the individual reviewers. 

**Summary**: This work is a contribution to the idea of the Neural Sampling Code (NSC) in computational neuroscience. While it is a prominent theory, there exists no work that formulates it such that NSC models can be fit to recorded neurophysiological data. The central goal of our current work is to (1) show how to mathematically express NSC which allows us to (2) fit generative models under NSC to neural responses to natural images, and (3) quantitatively compare different models with standard metrics such as data-likelihood. This allowed us to develop and fit deep learning-based generative models that are more flexible than existing literature models, and they outperform them when applied to recorded V1 responses to natural images. We also show (4) that we can derive a response-predictive model from NSC using the learned generative model’s posterior, which can be directly compared to system identificatio. Equipped with this framework, we (5) show that NSC models are not powerful enough to account for V1 responses to natural images compared to system identification. This motivates future work that can close the gap between normative NSC theory and data driven system identification.

**Re 9, 10**: Currently, system identification models are the state-of-the-art in predicting responses to natural images. To our knowledge, our work is the first to show that generative models under NSC fail to predict neuronal responses well compared to system identification models. There is still much room to build better generative models which would improve the posterior performance and our work builds the quantitative foundation for it. 

But, for a _given dataset_ of responses and stimuli, we do _not_ expect the posterior from even the _true_ generative model to surpass the performance of the _best_ system identification model, because the generative model does not provide any advantage over system identification in response prediction, unless some specific inductive biases were to be introduced in the generative model. 

However, if we change the stimulus distribution $p(x)$ to $p\_{new}(x)$ with markedly different stimulus statistics and let the sensory neurons adapt to $p\_{new}(x)$, we would expect the system identification performance to drop on $p\_{new}(x)$ and the system identification model might have to be retrained on a new dataset of responses under $p\_{new}(x)$.

This is the case where we could expect the generative model to be beneficial. Specifically, change in $p(x)$ to $p_{new}(x)$ may entirely derive from the change in prior $p(z)$ to $p_{new}(z)$, while $p(x|z)$ remains fixed, since $p(x|z)$ hypothetically represents the invariant “physical” process by which the latents (e.g., identity of an animal) give rise to observations (e.g., appearance of the animal). Consequently, if NSC were to be true i.e, $r \sim p(z|x) \propto p(x|z)p(z)$, the neuronal adaptation can be accounted for via simply learning $p_{new}(z)$, i.e., $r_{new} \sim p_{new}(z|x) \propto p(x|z)p_{new}(z)$, keeping $p(x|z)$ fixed. Importantly, $p_{new}(z)$ can be learned from the new stimulus samples alone from $p_{new}(x)$ _without the need of a new dataset of responses to the new stimulus samples_. 

This insight helps us identify potential future experiments to rigorously test NSC models utilizing our framework since our it lets us learn the generative model $p(x|z)p(z)$ via $p(x|r)p(r)$ (NSC assumption). Namely, we would perform experiments in which we let the neural population adapt to different sensory contexts with expected shifts in $p(z)$. Using our NSC framework, we would expect to be able to predict how these neuronal responses should change (as reflected in updated $r_{new} \sim p_{new}(z|x)$) under new contexts.

Meta-review: This paper describes a new data-driven approach for quantitatively assessing the neural sampling hypothesis. This question is timely and important, but the paper elicited a major difference in opinion from the four reviewers, with three giving a "weak accept" and one giving a "strong reject" score. After reading the reviews, author rebuttal, and discussion, I am inclined to side with the three reviewers who voted for acceptance. Reviewer GaXp raises a worthwhile point in criticizing the lack of biologically plausible implementations of neural sampling codes, which I feel the authors should certainly mention in their discussion as a shortcoming of the theory and an important avenue for future work. However, I ultimately agree with the authors that the lack of such an implementation is not a fatal flaw, and that a framework for quantitatively assessing the (normative) sampling hypothesis based codes is a major contribution in its own right. (Indeed, if we go back to the original sparse coding papers from Olshausen & Field, the normative theory made an important contribution despite the lack of plausible implementation; in my reading, it was only later that papers such as Rozzell et al 2008 took the additional steps toward proposing biologically realistic implementations of sparse coding). So ultimately I would like to recommend this paper for acceptance. However, all reviewers raised concerns about clarity (e.g., about the relationship between latents z and r), and so I would urge the authors to revise the paper to address all reviewer comments carefully and thoroughly.
Decision: Accept (poster)