Neural Sampling in Hierarchical Exponential-family Energy-based Models
Authors: Xingsi Dong
Si Wu
Venue: NeurIPS 2023 poster (2023)
Keywords: Baysian brain
sampling-based inference
energy-based models
local learning
exponential-family
TL;DR: We propose a hierarchical EBM whose learning process is localized both in time and space.

Abstract: Bayesian brain theory suggests that the brain employs generative models to understand the external world. The sampling-based perspective posits that the brain infers the posterior distribution through samples of stochastic neuronal responses. Additionally, the brain continually updates its generative model to approach the true distribution of the external world. In this study, we introduce the Hierarchical Exponential-family Energy-based (HEE) model, which captures the dynamics of inference and learning. In the HEE model, we decompose the partition function into individual layers and leverage a group of neurons with shorter time constants to sample the gradient of the decomposed normalization term. This allows our model to estimate the partition function and perform inference simultaneously, circumventing the negative phase encountered in conventional energy-based models (EBMs). As a result, the learning process is localized both in time and space, and the model is easy to converge. To match the brain's rapid computation, we demonstrate that neural adaptation can serve as a momentum term, significantly accelerating the inference process. On natural image datasets, our model exhibits representations akin to those observed in the biological visual system. Furthermore, for the machine learning community, our model can generate observations through joint or marginal generation. We show that marginal generation outperforms joint generation and achieves performance on par with other EBMs.

Neural Sampling in Hierarchical Exponential-family Energy-based Models
Xingsi Dong .2.3
Si Wu 1..33 siwu@pku eeweppuuuuuuuu.....
dxs19980605@pku.com edw.:e000@puuuuuuccr..... edu. cn
edu. cn
1. PKU-Tsinghua Center for Life Sciences, Academy for Advanced Interdisciplinary Studies.
2. School of Psychological and Cognitive Sciences,
Beijing Key Laboratory of Behavior and Mental Health, Peking University.
IDG/MeGovern Institute for Brain Research. Center of Quantitative Biology, Peking University.
Abstract
Bayesian brain theory suggests that the brain employs generative models to un- derstand the external world. The campling-basete perspective posits that the brain infers the posterior distribution through samples of stochastic neuronal responses. Additionally, the brain continually updates its generative model to approach the true distribution of the external world. In this study, we introduce the Hierarchical eaponential-famile Energy- Energy-based (HEE) model, which captures the dynamics of inference and learning. In the HEE model, we decompose the partition function into individual layers and leverage group of neurons with shorter time constants to sample the gradient of the decomposed normalization term. This allows our model estimate the partition function and perform inference simultaneously, circumventing the negative phase encountered in conventional energy-based mod- els (EBMs). As a result, the learning process localized both in time and space, and the model is easy to converge. To match the brain's rapid computation, we demonstrate that neural adaptation can serve accelerating the inference process. On natural image datasets, our model exhibits representations akin those observed in the biological visual system. Furthermore, for the machine learning community, our model can generate observations through joint or marginal generation. We show that marginal generation outperforms joint generation and achieves performance on par with other EBMs.
momentum term, significantly
Introduction
Human behavioral studies 1 2
and animal neurophysiological studies 4 have suggested that the brain performs statistically optimal Bayesian inference to interpret the external world [6] 7 g 9
One promising theory for implementing Bayesian inference in the brain is to interpret the variability of neural responses as Monte Carlo sampling of the posterior distribution [10]. This perspective
naturally accounts for the irregular firing patterns and other response properties observed sensory
cortex neurons [ 12 13]. Numerous sampling-basse models [10
14 [3] 16 12 13 [8] have
been proposed to elucidate neural dynamics and the underlying mechanisms of Bayesian inference.
To facilitate the brain's ability to derive meaningful representations from sensory input, it must continually update its generative model to approximate the true distribution of the external world [19]. Previous approaches often neglect this critical learning process. They either maintain fixed parameters for their generative model [15] [ [6), or they employ biologically implausible methods like the variational approach with backpropagation (BP) for training [9]. Is there a generative model that integrates sampling-based inference with the capability learn locally both time and space?
37th Conference on Neural Information Processing Systems (NcurIPS 2023).

---

Energy- based models (EBMs) (0 learning with spatially localized rules. The reason for non-locality in time is the need to estimate the partition function. To estimate the global partition function, the network has to perform a top- down pass to obtain negative sample. This process is referred to as the negative phase. During this pass, it disrupts the neural network's s stored inference results, which is essentially the same reason for the temporal non-locality as in the case of BP. Recently, predictive coding networks (PCNs) [21 22 23 24 use the Gaussian distribution to avoid the negative phase since the partition function of Gaussian is constant. But further study (5 reveals that the Gaussian assumption is restrictive when dealing with complex probability distributions. Moreover, setting aside biological constraints, estimating the partition function itself poses a significant challenge. In the field of machine learning, various sampling methods, such as amortized generation 26 and implicit generation 27) have been proposed to tackle this issue. However, both of these methods involve the use of the negative phase, which is known for being challenging converge. Generative adversarial networks (GANs) address this issue by utilizing a discriminative model, but GANs are notoriously difficult to train in practice.
provide a framework for inference with sampling method and
Besides, the inference dynamic of Gibbs sampling 28 20 or Langevin sampling 22 in EBMs essentially perform random walks in local regions rather than the whole posterior space, which is too slow to be compatible with brain functions [16). Therefore, it is crucial to investigate whether neural basedd circuits in the brain have the capacity of realizing sampling- based inference rapidly.
Summary of the work. In See we propose that our brain holds an EBM as the intrinsic generative model to interpret the external world. The neural dynamics employ a sampling baaeed method for Bayesian inference. Simultaneously. the learning dynamic aims to minimize the discrepancy between the observed distribution of intrinsic model and the real world. In Section 3 we introduce the Hierarchical Exponential-family Energy-based (HEE) model, which allocates the partition function across each layer. This allocation shifts the estimation of the total sample space required for calculating the partition function from a product of individual layer spaces to a sum of layer spaces. This approach enhances the convergence of our model. Furthermore, we efficiently sample the normalization term of the exponential-fane fammly in each layer using a group of neurons with fast dynamics. This localizes the learning process in both time and space. In See we find that incorporating noisy adaptation, a generic feature of neuronal responses, into the inference dynamics effectively yields second-order Langevin dynamic. In Sec 5 we validate the capabilities of the HEE model using 2D synthetic datasets and FashionMNIST [30] Then, we incorporate receptive field as the biological constrains to the HEE model training CIFARIO (3. We show that the HEE model can achieve performance comparable with previous EBMs (27. We also investigate the neural representation of semantic information, including orientation, color and category, which exhibit similarities to biological visual systems. And the neural adaptation can trigger neural phenomena including oscillations and transient which are widely observed in biological systems. In Sec 6 we discuss several related theories and models.
Main Contributions. We propose a hierarchical EBM whose learning process is localized both in time and space, which could potentially serve as a mechanism for the brain to utilize changes in synaptic strength for learning. And our brain-inspired EBM also presents technique for estimating the partition function, which is a challenging problem within the machine learning community.
2 The intrinsic generative model
In this section, we propose that our brain holds energy-based models (EBMs) as the intrinsic generative model consisting of two components: inference and learning. Inference is believed to be carried out through neural sampling [10 20). while learning is accomplished through long-temn term synaptic plasticity 32)
Let be the observation received by brain, and let be the latent variable represented by neurons. The joint distribution of the EBMs is written as Po (xx z) = po [x] (/z))o(z)) where 0 are stored in the connection weights of neurons (IA) The EBMs aims to minimize the difference between the
intrinsic marginal distribution pe (x) and the true distribution of the external world Prrue (x), which is described by the Kullback-Leibler divergence,
min DKI. |Pimee || (poox)) 0
(I)
2

---

(A) Probabilistic model
(B)
Inference learning
(C)
Joint generation
(D)
Marginal generation
x
x
x
pe (x(2)
4
z
z
(z(
Figure 1: (A) The directed graphical model of the energy-based model (EBM). (B) The latent variable receives the likelihood information po (=(x)z from the observation x and combines with the prior knowledge po (z) to perform the inference dynamic. And the connected weights changes following the learning dynamic (dashed line). (C) The latent variable performs inference dynamic and the observation x performs generation dynamic, which leads to the distribution of x, 2 Pe 2). And the connection weights
fixed (solid line). (D) The latent variable doesn't 't receive likelihood
information from observation and (z(z), x pe (x) which is called marginal generation.
The neural system can adopts the gradient based learning method such as gradient decent. And the gradients is calculated (see SI for detailed proof),
V DKL [Pmne (x) pe (x)] Ee (x) (x)
Ve in (x. z)]
(2)
The second expectation of the above equation requires the posterior of a given observation x, which is calculated as pe (z|xx) = (x (x)))(x). Practically, the posterior is intractable for the denominator (x) requires complex integral. Variational inference is usually used solve this problem (such as VAE (33) And the EBMs adopts the sampling method to avoid complex calculations. By leveraging the relationship V. In po (z/x) 2 In (x, z). ,the samples of posterior can be offered by the neural dynamic (Langevin sampling),
dz
dt = V. In z) + V
where Ç is Gaussian white noise and Ta is the time constant. The stationary distribution of the above dynamic is our target distribution [z(2)x] x)). The sampling algorithm, along with the joint distribution, determines the connections of neurons and their inference dynamic (B).
For the learning dynamic, the brain receives the observations from the real world continuously (Ex-Praeexx)) and the neural dynamic mentioned above can produce samples of their posterior simul- taneously (Fzz)p(())))) The parameters can be updated according to the gradient, which often implemented by Hebbian learning [
d0 dt
DKL [Pave (x) || pa (x) = Ve In (x,,
(4)
where т is the time constant of synapses.
For neuroscience, this is the end of the story. Nevertheless, EBMs can also generate observations which the machine learning society follow with interest. In the generation process, the observation is fixed but undergoes the Langevin sampling,
dx
Vx Inpo(xxx V
(5)
The latent variable can follow the inference dynamic Eq. (3.. In this case, the observation x and
latent variable z together follow the joint distribution (x, 2). Thus, the observation x can produce order to get marginal distribution po (x), the latent variable z just needs to follow the prior distribution pa (z) rather than reaching the posterior. this case, the generation dynamic of latent variable written as,
samples following (x), which is called joint generation (C.. However,
dz V, (2) + V2T
To dt
(6)
3

---

And the stationary distribution of Eg.(5) performed by observation x still equals to pe (x), which called marginal generation (Fig DiigDp)...o S we will see that the marginal generation performs better than joint generation. Here is an informal understanding. The process of sampling (x, 2) can be understood as searching for a specifie pair. We assume that the size of the x-space is ((m) and the size of the z-Space is (n).. In joint generation, the search is conducted simultaneously in both the x-space and z-space, resulting in a required search space of O(n * m.). In marginal generation. the process involves initially searching in the z-space according to pa (z). Once z is found, it is fixed. This step's s search space size is (m)) Then, x is searched based on po (x|z)) ()) in the x-space. This step's 's search space size is (m). leading to combined required search space size of O( + m).
3 Exponential-family energy-based model
In this section, we provide a neural implementation of the HEE model and outline the specific dynamics involved in inference, learning, and generation, as discussed in Section [ Approximating the target distribution (x) which is diverse and complex requires a good representation ability of the model. Exponential families include many of the most common distributions (such normal, Poisson, gamma distribution and so on). Morcover, exponential families can be easily parameterized, allowing for generalization and flexibility in modeling various types of distributions.
Let Xg E Ro be the observation (such as an image) received by our brain. And there are layers of neurons representing the latent variables X1:1 = (x1, X2, ..... ) E R The joint distribution is a Markov chain (Fig 2 starting with p(xr.) = exp |n. (xx,) + g(x))) - A(nr)].
L-1
(x) = p(xi)
II
po x(+1), 1=0
(x(x| = exp [yi o(xi) + g(xi) - A(m)] (7)
The natural parameter np E R is function of X1+1 with parameters 0: R which is written as 17; = (x where f(-) is the activation function. And 17 is constant. The sufficient statistic o(x)) H and the base measure g(xi) IR is the function of x1. A(T) E E is the normalize term (log-partition function) to make sure the sum of the probability equals to 1. In order to get the inference dynamic, we substitute the joint distribution Eq.47) into the Langevin dynamic 33 obtaining,
dx: dt
Ts
((x)))[[[((xxxoo - A (m) + o' (x))n + (xx) + 2T,
(8)
where fxxx), o' (xi) R are diagonal matrices. The derivative of pog-parttiion A' A'n (m is
intractable for it needs complex integral. Here, we use group of interneurons EI- E " to rep-
resent the term ((x(-1) - 11) It can proved that A' 4'(7-11- = Ex--ipp (x= |xx) (o(x)-1) 1)] (See SI for detailed proof). Thus, in order to calculate A'( the interneurons need to produce samples X following the distribution - |xx) in short time compared with the inference dynamic. Therefore, the dynamic of E can be written as,
du: = d' ( - + g'(u-1) + V2T
E1 = <b (x <5 U:
Tu
dt
By setting the time constant T. << Ta we can ensure that uz-1 converges much faster than x1, and the stationary distribution of u corresponds to po (u:-1/x)) This leads to E ((x(-1) A' A'( (1)-
Now, we can rewrite the inference dynamic Eq. into,
Ts dx = f'xx + o' ('(x)) + g + 2
(10)
dt The first term the right side of the dynamic equation indicates that neurons x; receive feedback from interneurons 1;- which provides likelihood information po (x/-1|x)), The second term shows
'There are various types of interneurons that target on pyramidal cells, comprising approximately 10-20% of the overall neuron population in the cerebral cortex [34]. The interneurons in the HEE model bear the closest resemblance the Large Basket Cell or Nest Basket Cell 35) which collectively constitute around 50% of interneurons. Their clectrophysiological characteristics include fast spiking, non- accommodating, and non- adapting behaviors. These interneurons have also been identified in the visual cortex of ferrets (36). displaying short-duratios action potentials (approximately 0.5 ms at half height). This suggests that these neurons have shorter time constants compared to pyramidal cells.

---

(A)
Probabilistic model
BB)
Inference & learning
(C)
g(-)
(xx-1)
X1-1
X/-1
El-1
Mp-11
Pa (x-- (x1)
((x))
X1
X
13
m
pa (x) x(+1)
1+1
11
X1+1
X(+1
Figure 2: (A) The directed graphical model of the hierarchical exponential-family energy-based (HEE) model.. (B) The inference and learning dynamic of HEE model. The red arrows represent the likelihood information and the black arrows represent the prior information. Neurons X receive the likelihood information from E and receive prior information
from E₁. The interneurons Et-1 receive the natural parameter 7/-1 from neurons X; and perform the Langevin sampling dynamic to approximate A' (1) Then, the interneurons compare it with the the sufficient statistic é(x--1) received from X;-1 to calculate the value of et-1- The dashed line shows that the connection weights 0 will perform gradient decent in the learning dynamic. (C) The neural connection diagram between layer and layer +
that neurons X₁ receive prior knowledge (x(1)x14) from interneurons E, through a feedforward loop. The term (xx) controls the self-connections within layer (Fig 2B). E1-1 is also called the error term in PCNs. Here, we show that the predictions in PCNs essentially estimate the log peg-pariiioon.
Then, we can obtain the learning dynamic by substituting the joint distribution Eq. into the gradient decent dynamic Eq. 4).
dA [x(x)) A'(n)) T8 dt
=x((x)+1)*
(11)
The derivatives of the log peg-partiiion function are stored the interneurons E₁. As result, the synaptic changes are determined solely by local neurons, adhering to Hebbian rules.
After inference and learning, our model can also generate observations. During joint generation, the dynamics of neurons X1:1. follow the same principles as the inference dynamics described by (T). By substituting the joint distribution (q.(7) into the Langevin dynamic Eq.(5). samples from the model by.
generate new
dxo ø' (xo)To + g (xxx + 22
To dt
(12)
And for marginal generation, we can substitute the prior distribution described Eq. into Eq.(6) obtaining the dynamic of neurons XI:L-
dx
o'(xi m + gxi) + V2TzE.
(13)
dt
Neural adaptation accelerate the sampling process
Langevin sampling essentially performs random walks in local regions rather than the whole posterior space 37, because the drift term V. InInp lo z) in (3. will vanish near the local minima and only noise term remains (Ai Sampling the entire posterior space is time-consuming process and does not align with the brain't 's ability to perform tasks quickly. Therefore, it is important to implement faster sampling algorithms for our model.
In this section, we show that by including noisy adaptation, the network able to speed up the inference dynamic significantly. Adaptation is common phenomenon observed in neural systems,

---

(A)
(B)
second-order
(C)
(D)
Langevin sampling
Langevin sampling
Tesc
A
Trec
Time M(/r))
Energy function
0.5
In (z(x)
Layers L
Figure 3: (A)BB The sampling trace of Langevin dynamic and second -odder Langevin dynamic of mixture of 2D Gaussian distribution. (C) Illustration of the sampling process. States with lower energy indicates higher probability, which needs to be stay longer. And the network needs to cross the energy barrier to reach new local minima for the non-convex of the energy function. (D) X and det H)) increases and decreases, respectively, with the changes of layers L and the total number of neurons fixed (2 Th 10000).
where negative feedback mechanisms are employed suppress neuronal responses when they reach high levels. Here, we show that the neural adaptation can introduce an auxiliary variable to implement faster sampling algorithm, which is called second-order Langevin dynamics (SLD),
dz
Tz di. V2 (n v V2T2E.
(14)
dv
Tv dt
mV2T.
(15)
where >> 0 controls the adaptation strength. The auxiliary variable can keep the network moving
while reaching the local minima, which essentially serves as the momentum term 3B. When there is adaptation (m 0), the above dynamic will degenerate to the Langevin sampling (LS) described in Eq. And the stationary distribution of the above dynamic remains to be po (z|x) (See SI for proof). By substituting the joint distribution of exponential family Eq. into the above dynamic, we can obtain the SLDD inference dynamic,
dx:
T= ('(x)) 1 d' d'(x)) g'(xi) vi V2T.EI
(16)
dt
mv. To dt 2 m
dv:
(17)
Here we exemplify the neural adaptation with spike frequency adaptation (SFA) [38]. In the neural system, the adaptation current vi accumulates the noise term & coming from the concentrations, release of neural transmitters, activation/inactiiation of ion channels and so on, which is described by Eq. 17) The adaptation current induces suppression on neurons X;, acting as momentum variables to accelerate the sampling process (16).
We conduct further investigation elucidate the precise mechanism by which noisy adaptation facilitates the acceleration of the sampling process in the HEE model. Considering that the energy function In pa po(x), (1. (x) non-convex, the sampling process
be divided into parts. In the first part, the network needs to find local minima and samples near the local minima, which takes a certain amount of time called recurrence time Tree (Ci In the second part, the network needs leave the local minima and find a new one, which takes a certain amount of time called the escape time Teac- Typically, we have Tree & Teac- The total time to get stationary distribution can be approximated by T Texe + Tree- It can be proved that (39] the recurrence time is bounded by Tree (1))( (H,)
the smallest eigenvalue of the matrix H.,.
H H/Ts 0
1 1/T m2/(2T))
(18)
where H is the Hessien matrix of the energy function - In po (x.1xx.//X) The smallest eigenvalue of H., calculated as A1 (H,) min(Ae (H)/Tx m2/(2r.)) Thus, in the case m > 2\ (H))// the SLD can reduce the recurrence time Trec to accelerate the sampling process. And the escape

---

Preaf x
HEE-L Joint
Marginal
HEE-LL
HEE-NL Marginal
HEE-NL-A HEE-L-A
Figure 4: Evaluation on 2D synthetic datasets and FashionMNIST: a mixture of four Gaussian distribution (first line), a mixture of four Sanana-shaped distribution (second line), pinwheel-shaped distribution (third line).
time is bounded by Tese O (VT/dee(HH)) The determinant of H. is calculated det(H)) det(H) (mn) ( /2r,)_ ni Thus, det(H)) is monotonically increaseing with m, indicating that the larger m is, the shorter the time it takes to escape from the local minima.
The analysis presented above demonstrates that the convergence speed of the inference dynamic is determined by the values of (H) and det H). This valuable insight can be leveraged to guide the design of the network architecture, enabling the creation of more efficient and effective models. Specifically, with a fixed total number of neurons Σ, n4, increasing the number layers L results in deeper network, while decreasing the number of layers results wider network. Practically, show that del((H) will decrease with layers while A will increase (Fig(D), which indicates that there is trade-off between the recurrence time Tree and the escaping time Tesc with different layers (See SI for detailed setting and analysis).
5 Experiment
In this section, firstly validate the capability of HEE model for approaching complex distribution by examining the quality of generation. Then, we show that our model demonstrates similarity in the representation of natural images the biological visual system. And adaptation induce oscillatory behavior and transient overshoots in neurons during the inference phase.
5.1 Generation
Firstly, we conducted experiments using three variations of the HEE model, each with different Φ(x) functions and sampling methods (Tabel []) evaluate their capabilities. The experiments were performed on both 2D synthetic datasets and the FashionMINST.. We use the fully connected architecture, i.e., has
elements. The results (second and third column) show that HEE with linear statistic (HEE-L. struggles to capture the complex distribution. Moreover, we theoretically prove that HEE-L can only approach unimodal distributions (See SI for detailed proof). For HEE-NL, some modes are missing while using the joint generation. And when the spacing between modes large, there is an issue of non-uniformity among different modes. And the marginal generation converge much faster than the joint generation. In FashionMNIST, it takes less time for marginal method to get the generation high-quality images.

---

Then, we employ the HEE-NL-A with layers 10 on the CIFARIO unconditional. The sparse connection is employed as method to mimic the receptive field behavior found in biological systems. We quantitatively evaluate image quality of HEE-NL-A with Inception score [40] and FID score 4I in
Model
<(x)
Sampling method
HEE-L
I
LS
HEE-NL
sigmoid((s) sigmoid(x)
L.S
HEE-NL-A
SLD
Table 1: Table of different HEE models.
T Overall, we achieve a performance comparable with the previous EBMs. And the generation quality of the marginal method is better than joint method, which agrees with the previous results 42]
5.2 Inference
We further use the HEE-NL-A trained on CIFARIO to explore the relationship between the latent features and semantic information, including orientation, color and category.
Orientation: Simple cells 43 and complex cells 4 are the most prominent and widely observed neurons in the biological visual system that exhibit tuning to orientation. They are found in the primary visual cortex (V1) of numerous animal species 45 46 We present the model with gabor images of different orientations commonly used in experiments (A) and compute the mean and variance of the neural responses for each neuron. Then, we use Gaussian curve with bandwidth limited from 20" to 90" and two-modes Gaussian curve to fit the simple cell and complex cell, respectively (Fig B).. We find that the proportion of simple cells and complex cells remains relatively consistent across each layer and both decrease with increasing layers in our model (Fig@C)
End stopping: In the HEE model, interneurons essentially represent the error term in the PCNs. We have observed the phenomenon of end stopping' in interneurons (Fig 6D), which aligns with the end-stopping behavior observed in error neurons in the PCNs (21]
Color: Recent study shows that there a hierarchical representation for chromatic processing across the ventral pathway of macaque [47. We present model using reshaped natural images (48) and employ principal component analysis to demonstrate that the middle layer's neural representation's most informative dimensions carry chromatic information (Fig
Category: Visual object recognition believed to be solved by the brain hierarchically [49] A recent study [50] demonstrate that the inferotemporal cortex, situated in the deeper layer of the visual pathway, is capable of constructing a lincar map of the object space. For each layer in our model, employed a linear support vector machine (SVM) to classify the ten labels of the CIFARIO. The SVM was trained using the average neural responses as features. Fig(6F illustrates the projection of the features in the last layer onto the SVM weights corresponding to the "cat" and "dog" labels. Furthermore, we observe that the classification accuracy improves as we move up the layers of our model (Fig (G) which is consistent with findings in both biological visual systems (50] and the artificial neural networks [51]
Phenomena:
Oscillations [52] and transients [53 are two kinds of spatial-temporal dynamic features
in neural systems, which play a crucial role during the sampling process [3] Here, we show that by adjusting the adaptation strength m, the oscillation frequency of the HEE model can span within the
Model
IS
FID
HEE-NL-A (Joint)
5.95
43.21
EBM (single) 27
6.02
40.58
HEE-NL-A (Marginal)
6.47
37.05
MEG (Generator)
42
6.49
35,02
EBM (10 ensemble)
6.78
38.20
MEG (MCMC)
7.31
33.18
Table 2: Table of Inception and FID scores.
Figure 5: Marginal generation of HEE-NL-A on CIFAR10.
8

---

(A)
Simple cell
(B)
Complex cell
(C)
Orientation luning
0.3
simgle cell
0.5
0.5
complex
Firing rate
Propartion
0.2
0.1
0.0
0.0
4
B
Layer
(D)
Firing rate
10
20
Bar length x (pixels)
(G)
Decoding accuracy
(H)
Oscillation peak
Transient
0.8
Frequency (Hz))
10"
Average step size
Accuracy
0.6-
10"
0.4
10
0.2
2
2
Layer
ml(T./T))
mlTs/T))
Figure 6: (A)(B) The red dots show the average firing rate of two neurons in x with different gabor-like stimulus. Different tuning curves are used to fit simple cells and complex cells. (C)
The proportion of simple cells and complex cells in each layer. (D) We show horizontal bars of varying lengths to the HEE, The red dots show the average firing rate of neuron in E₁ whose corresponding neuron in X simple cell preferring 0 degree. (E) The images are plotted the location corresponding to the projection of their average neural response in layer 5 onto the first two principal components, Red images are located in the upper left corner, while blue-green images are located in the lower right corner. (F) The true label of cat, dog, and other categories should be respectively located the second, fourth, and third quadrants. (G) The classification accuracy SVM in each layer. (H) We sampled and statistically analyzed the distribution of the highest firing neurons during the inference phase in the first 100T, for different values of The gamma band is centered around the dashed line. (I) We sampled and statistically assessed the maximum change in firing rates of all neurons during the inference phase in the first 100r, for different values of We refer to the mean of the maximum change values the 'average step size'. We utilize the average step size of the ncurons during the sampling process as an indicator of transients.
rate frequencies
range of 20-80 Hz (gamma band), which is widely observed visual systems 54 stimulus-onset transients of the firing rate can also be enhanced by the adaptation (Fig
Fig 6
(H). And
Discussion
The present study investigates the sampling-based inference and learning dynamic within the frame- work of intrinsic generative model. We introduce the HEE model neural implementation that utilizes neural dynamics and Hebbian learning. Additionally, demonstrate that the inclusion of

---

neural adaptation can significantly accelerate the sampling process and give rise to various dynamic phenomena throughout the network. this section, we will discuss several related theories and models.
Probabilistic Population Code (PPC) 55 56 is another theory that explains how the brain perform Bayesian inference, in which neural responses are interpreted as the parameters of the probability distributions. We adopted the idea [57 that PPC theory incorporates two generative models. In the framework of PPC, the experimenter presents the subject with observation x (gabor image) based on the semantic information s (orientation), which actually defines an external generative model from to x. And the subjects holds an intrinsic generative model with latent variable z represented by neural response to interpret the observation x. PPC theory integrates two generative models into a single generative model, in which the neural response z is regarded as the observation generated from the semantic information s. We propose that the learning dynamic occurs exclusively within the intrinsic generative model, as the brain is not aware of the external generative model.
Energy-based models (EBMs) 20 58 When EBMs were initially proposed [20), they had latent variables corresponding to neurons. Later, to enhance the model's expressive power, hierarchical structures were introduced [59 60 Such EBMs could typically ensure local learning in space. As
artificial neural networks have become increasingly powerful, it has been observed that for generative tasks. there's no need to explicitly introduce neurons as latent variables within EBMs. Instead, one can directly employ a neural network to represent the energy (27). Training such EBMs often involves utilizing BP. The distinction between these EBMs and traditional EBMs is akin to the difference between dynamic systems and recurrent neural networks
Regardless of whether it's the traditional EBMs or the new type of EBMs, both involve the challenge of estimating the partition function. This difficulty arises from the fact that as the depth of the energy function increases, the total sample space required multiplies the space for each layer. In the case of HEE, we allocate the partition function across each layer. As a result, the total sample space required the sum of the spaces for each layer. This significantly reduces the required sample space.
Predictive coding networks (PCNs) 22 22 The interneurons in the HEE model serve a similar role to the prediction error in PCNs. And our theoretical analysis shows that the predictions in PCNs essentially represent the decomposed log- pegrarttiionn function. And PCNs don't stress the sampling- based inference, which requires them approximate the energy function using variation inference by delta function. Sampling- inference can assist the network in exploring the posterior probability space, leading to more accurate estimation of the energy function. Additionally, it can account for the observed neural variability in experiments.
Diffusion models (DDPMs) 66 62 The HEE model and DDPMs share the same joint distribution and both exhibit a hierarchical Markov structure, which may contribute to the HEE model's ' strong expressive potential. The marginal generation is also called latent space MCMC 663 42), which is similar to the generation process of DDPMs. While DDPMs unfolds the Markov chain over time, the HEE model unfold it between layers of neurons. However, in order to reach better performance, DDPMs use a fixed diffusion process as the inference dynamic, which may not be adopted by our brain since the latent variables in our brain carry semantic information (such as simple cells 43)
Speed up sampling 17 38 previous work, inhibition neurons were used to serve as momentum terms to accelerate sampling [17. However, this approach required a one-to-one correspondence between inhibition neurons and excitatory neurons. In our approach, we consider the adaptive properties inherent in each neuron itself to serve momentum, naturally resolving the one-to-one correspondence issue. Furthermore, our consideration extends to sampling in a non-convex energy space, which differs from the prior focus solely on convex space convergence properties (38).
Acknowledgement
I'd like to express my gratitude to Tianqiu Zhang and Chaoming Wang for their assistance in configuring the experimental environment. I also want to thank BrainPy [64] for their support throughout this work. Special thanks to Yumeng Cao for reviewing the article's grammar and expression. This work was supported by Science and Technology Innovation 2030-Brain Science and Brain-inspired Intelligence Project (No. 2021ZD0200204).
10

---

References
[i] Marc Ernst and Martin S Banks. Humans integrate visual and haptic information statistically optimal fashion. Nature, 415(6870):4294444 2002.
[2] Konrad P Körding and Daniel M Wolpert. Bayesian integration in sensorimotor learning. Nature, 427(6971):244-2444 2004.
[3] Wei Ji Ma, Vidhya Navalpakkam, Jeffrey M Beck, Ronald Van Den Berg, and Alexandre Pouget. Behavior and neural basis of near-optimal visual search. Nature neuroscience, 14(6):783, 2011.
[4] Yong Gu, Dora E. Angelaki, and Gregory C DeAngelis. Neural correlates of multisensory cue integration in macaque mstd. Nature Neuroscience, 11110)120111210, 2008.
[5] Christopher R Fetsch, Gregory C DeAngelis, and Dora E Angelaki. Bridging the gap between theories of sensory cue integration and the physiology of multisensory neurons. Nature Reviews Neuroscience, 14(6))429-442, 2013.
[6] Tai Sing Lee and David Mumford. Hierarchical bayesian inference in the visual cortex. JOSA
A, 20(7):1434 1448, 2003.
[7] David Knill and Alexandre Pouget. The bayesian brain: the role of uncertainty in neural coding and computation. TRENDS in Neurosciences, 27(12):112-11, 2004.
[8] Alan Yuille and Daniel Kersten. Vision bayesian inference: analysis by synthesis? Trends in cognitive sciences, 10(7):301 -308, 2006.
[9] James CR Whittington, Timothy H Muller, Shirley Mark, Guifen Chen, Caswell Barry, Neil Burgess, and Timothy EJ Behrens. The tolman--ichenben machine: unifying space and relational memory through generalization in the hippocampal formation. Cell, 183(5):1249- 1263, 2020.
[10] Patrik Hoyer and Aapo Hyvärinen. Interpreting neural response variability monte carlo sampling of the posterior. Advances in neural information processing systems, 15, 2002.
[11] Ralf M Haefner, Pietro Berkes, and József Fiser. Perceptual decision-making as probabilistic inference by neural sampling. Neuron, 90(3):649-660, 2016.
[12] Gergö Orbán, Pietro Berkes, József Fiser, and Máté Lengyel. Neural variability and sampling- based probabilistic representations in the visual cortex. Neuron, 92(2):530-543, 2016.
[13] Rodrigo Echeveste, Laurence Aitchison, Guillaume Hennequin, and Máté Lengyel. Cortical-like dynamics in recurrent circuits optimized for sampling-bas based probabilistic inference. Nature neuroscience, 23(9):1 138- 149, 2020.
[14] Lars Buesing, Johannes Bill, Bernhard Nessler, and Wolfgang Maass. Neural dynamics as sampling: a model for stochastic computation recurrent networks of spiking neurons. PLoS computational biology, 7(1 7(11)= 2011.
[15] Agnieszka Grabska-Barwinsk.. Barwinska. Jeff Beck, Alexandre Pouget, and Peter Latham. Demixing oasts-fast inference olfaction. Advances in Neural Information Processing Systems, 26, 2013.
[16] Guillaume Hennequin, Laurence Aitchison, and Máté Lengyel. Fast sampling- based inference in balanced neuronal networks. In Advances in neural information processing systems, pages 2240-220 2248, 2014.
[17] Laurence Aitchison and Máté Lengyel. The hamiltonian brain: efficient probabilistic inference with excitatory-inhibioryy neural circuit dynamics. PLoS computational biology. 12(12), 2016.
[18] Yang Qi and Pulin Gong. Fractional neural sampling theory of spatiotemporal probabilistic computations in neural circuits. Nature communications, 1-11), 2022.
[19] Karl J Friston and Cathy J Price. Dynamic representations and generative models of brain function. Brain research bulletin, (4(3):275-288) 2001.
[20] David H Ackley, Geoffrey E Hinton, and Terrence Sejnowski. learning algorithm for boltzmann machines. Cognitive science, 91):147-169, 1985.
[21] Rajesh PN Rao and Dana H Ballard. Predictive coding in the visual cortex: functional interpretation of some extra-classical receptive-field effects. Nature neuroscience, 2(1):79-87, 1999.
11

---

[22] James CR Whittington and Rafal Bogacz. An approximation of the error backpropagation algo- rithm in a predictive coding network with local hebbian synaptic plasticity. Neural computation. 29(5):1229-1200) 2017.
[23] Tommaso Salvatori, Yuhang Song, Yujian Hong, Lei Sha, Simon Frieder, Zhenghua Xu, Rafal Bogacz, and Thomas Lukasiewicz. Associative memories via predictive coding. Advances in neural information processing systems, 34:3874-3886 2021.
[24] Alexander Ororbia and Daniel Kifer. The neural coding framework for learning generative models. Nature communications, 13(1):2064, 2022.
(25) Luca Pinchetti, Tommaso Salvatori, Yordan Yordanov, Beren Millidge, Yuhang Song, and Thomas Lukasiewicz. Predictive coding beyond gaussian distributions. arXiv preprint arXiv:2211.03481, 2022.
[26] Taesup Kim and Yoshua Bengio. Deep directed generative models with energy-based probability estimation. arXiv preprint arXiv:1606.03453 2016.
[27] Yilun Du and Igor Mordatch. Implicit generation and modeling with energy based models. Advances in Neural Information Processing Systems, 32, 2019.
[28] Stuart Geman and Donald Geman. Stochastic relaxation, gibbs distributions, and the bayesian restoration of images. IEEE Transactions on pattern analysis and machine intelligence, (6):721- 741, 1984.
[29] Max Welling and Yee W Teh. Bayesian learning via stochastic gradient langevin dynamics. In Proceedings of the 28th international conference on machine learning (ICML-I1), pages 681-688 2011.
[30] Han Xiao, Kashif Rasul, and Roland Vollgraf. Fashion-mnist: a novel image dataset for benchmarking machine learning algorithms. arXiv preprint arXiv: Ttt::170 778800777, 2017.
[31] Alex Krizhevsky, Geoffrey Hinton, et al. Learning multiple layers of features from tiny images. 2009.
[32] Tim VP Bliss and Terje Lømo. Long-lasting potentiation of synaptic transmission in the dentate area of the anaesthetized rabbit following stimulation of the perforant path. The Journal of physiology, 232(():331-33
1973.
[33] Diederik P Kingma and Max Welling, Auto-encoding variational bayes. arXiv preprint arXiv:1312.6114, 2013.
[34] Therese Riedemann. Diversity and function of somatostatin expppesstttn--erreessin interneurons in the cerebral cortex. International journal of molecular sciences, 20( 12):2952, 2019.
[35] Henry Markram, Filif Muller, Srikanth Ramaswamy, Michael W Reimann, Marwan Abdellah, Carlos Aguado Sanchez, Anastasia Ailamaki, Lidia Alonso-Nanclares, Nicolas Antille, Selim Arsever, et al. Reconstruction and simulation of neocortical microcireuitry. Cell, 163(())456- 492, 2015.
[36] Vanessa F Descalzo, Lionel G Nowak, Joshua C Brumberg, David A McCormick, and Maria V Sanchez- Slow adaptation in fast-spiking neurons of visual cortex. Journal of neurophysi- ology, 1312):1111-11188 2005.
[37] Xiang Cheng and Peter Bartlett. Convergence of langevin meme in kl-divergence. In Algorithmic Learning Theory, pages 186-211. PMLR, 2018.
[38] Xingsi Dong, Zilong Ji, Tianhao Chu, Tiejun Huang, Wenhao Zhang, and Si Wu. Adaptation accelerating sampling- based bayesian inference in attractor neural networks. Advances in Neural Information Processing Systems, 15:215344215444 2022.
[39] Xuefeng Gao, Mert Gurbuzbalaban, and Lingjiong Zhu. Breaking reversibility accelerates langevin dynamics for global non-convex optimization. arXiv preprint arXiv: 2018.
[40] Tim Salimans, Ian Goodfellow, Wojciech Zaremba, Vicki Cheung, Alec Radford, and Xi Chen. Improved techniques for training gans, Advances in neural information processing systems, 29, 2016.
[41] Martin Heusel, Hubert Ramsauer, Thomas Unterthiner, Bernhard Nessler, and Sepp Hochreiter. Gans trained by a two time-scale update rule converge to a local nash equilibrium. Advances in neural information processing systems, 30, 2017.
12

---

[42] Rithesh Kumar, Sherjil Ozair, Anirudh Goyal, Aaron Courville, and Yoshua Bengio. Maximum entropy generators for energy- based models. arXiv preprint arXiv: 085008 2019.
[43] David H Hubel and Torsten Wiesel. Receptive fields of single neurones the cat's striate cortex. The Journal of physiology, 148(3):574, 1959.
[44] David H Hubel and Torsten N Wiesel. Receptive fields, binocular interaction and functional architecture in the cat's visual cortex. The Journal of physiology, 160(1))106, 1962.
[45] PO Bishop, J So Coombs, and GH Henry. Receptive fields of simple cells in the cat striate cortex. The Journal of physiology, 23 231(1):3,-60, 1973.
[46] David H Hubel and Torsten N Wiesel. Brain and visual perception: the story of 25-year collaboration. Oxford University Press, 2004.
[47] Ye Liu, Ming Li, Xian Zhang, Yiliang Lu, Hongliang Gong, Jiapeng Yin, Zheyuan Chen, Liling Qian, Yupeng Yang, Ian Max Andolina, et al. Hierarchical representation for chromatic processing across macaque v1, v2, and v4. Neuron, 108(3):538-55) -500, 2020.
[48] Jia Deng, Wei Dong, Richard Socher, Li-Jia Li, Kai Li, and Li Fei-Fei. Imagenet: A large- scale hierarchical image database. In 2009 IEEE conference on computer vision and pattern recognition, pages 248-255. leee, 2009.
[49] James DiCarlo, Davide Zoccolan, and Nicole C Rust. How does the brain solve visual object recognition? Neuron, 73(3):415-434, 2012
[50] Pinglei Bao, Liang She, Mason McGill, and Doris Y Tsao. A map of object space in primate 2020.
inferotemporal cortex. Nature,
[51] Sue Yeon Chung, Daniel D Lee, and Haim Sompolinsky, Classification and geometry of general perceptual manifolds. Physical Review X, 8(3):031033, 2018.
[52] Supratim Ray and John HR Maunsell. Differences in gamma frequencies across visual cortex restrict their possible use in computation. Newron, 65:558888 2010.
[53] Bilal Haider, Michael Häusser, and Matteo Carandini. Inhibition dominates sensory responses in the awake cortex. Nature, 493(7430);97-100, 2013.
[54] Mark Roberts, Eric Lowet, Nicolas M Brunet, Marije Ter Wal, Paul Tiesinga, Pascal Fries, and Peter De Weerd. Robust gamma coherence between macaque v1 and v2 by dynamic frequency matching. Neuron, 7833)552--53,, 2013.
[55] Jeffrey M Beck, Wei Ji Ma, Roozbeh Kiani, Tim Hanks, Anne K Churchland, Jamie Roitman, Michael N Shadlen, Peter E Latham, and Alexandre Pouget. Probabilistic population codes for bayesian decision making. Neuron, 60(6): 1400::142-111
2008
[56] Jeff Beck, Alexandre Pouget, and Katherine A Heller. Complex inference in neural circuits with probabilistic population codes and topic models. Advances in neural information processing systems, 25, 2012.
[57] Sabyasachi Shivkumar, Richard Lange, Ankani Chattoraj, and Ralf Haefner. A probabilistic population code based on neural samples. In Advances in Neural Information Processing Systems, pages 7070-7079, 2018.
(58) Ian Goodfellow, Jean Pouget-Abadie, Mehdi Mirza, Bing Xu, David Warde-Farley, Sherjil Ozair, Aaron Courville, and Yoshua Bengio. Generative adversarial networks. Communications of the ACM, 63(1 1):1391144, 2020.
[59] Ruslan Salakhutdinov and Geoffrey Hinton. Deep boltzmann machines. In Artificial intelligence and statistics, pages 448 448.455. PMLR, 2009.
[60] Benjamin Scellier and Yoshua Bengio. Equilibrium propagation: Bridging the gap between energy-based based models and backpropagation. Frontiers in computational neuroscience, 11:24, 2017.
[61] Jonathan Ho, Ajay Jain, and Pieter Abbeel. Denoising diffusion probabilistic models. Advances in Neural Information Processing Systems, 33:6840-6851, 2020.
[62] Yang Song, Jascha Sohl-Dickstein, Diederik P Kingma, Abhishek Kumar, Stefano Ermon, and Ben Poole. Score-based generative modeling through stochastic differential equations. arXiv preprint arXiv:2011.13456 2020.
13

---

[63] Yoshua Bengio, Grégoire Mesnil, Yann Dauphin, and Salah Rifai. Better mixing via deep representations. In International conference on machine learning, pages 552-560. PMLR, 2013.
[64] Chaoming Wang, Xiaoyu Chen, Tianqiu Zhang, and Si Wu. Brainpy: a flexible, integrative, efficient, and extensible framework towards general-purpon brain dynamics programming. bioRxiv, pages 2022-10, 2022.
14


Reviews:
[{"review_id": "pXdVL3dmjs", "reviewer": "Reviewer_zTHc", "summary": "In this submission, the authors introduce a biologically plausible method to train a hierarchical energy-based model as well as perform inference over it, via Langevin sampling. This accomplishes a goal that has long been sought in computational neuroscience, namely, neurally plausible methods for inference over a hierarchical generative model of the world.\n\nPreviously, the main problem facing inference in hierarchical energy-based models is the difficulty of evaluating the partition function. The authors address this issue by modeling the partition function with Langevin sampling, but with a faster timescale than the Langevin sampling used for inference more generally. This appears to work well empirically. Also underlying their model's impressive generation performance is the use of a momentum term in the Langevin sampling, similar to Hamiltonian Monte Carlo, which they associate with short-term adaptation in the brain. Learning in this system is Hebbian and biologically plausible, and additionally the authors show some rough parallels with biological phenomena, such that the inferred activations of selected neurons have orientation- and hue-specific tuning curves.", "questions": " - How many long does convergence take to produce the samples shown in the manuscript?\n - What is the connectivity of the interneurons? What cell types might these be? What evidence is there for their faster time constants?\n - I don't fully understand the circuit diagram for biological neurons. Could a schematic be drawn with individual neurons, complete with a table of the key properties of those neurons in this theory?", "limitations": "There was not much mention of any limitations of this method. What are the least biologically plausible aspects? What aspects need to be experimentally confirmed?", "rating": 5, "confidence": 3, "soundness": 3, "presentation": 1, "contribution": 3, "strengths": "I was impressed with the quality of samples generated by this hierarchical EBM, and all the more so because the learning process is so simple and Hebbian in nature. The model aligns well with a wide literature in computational neuroscience, including hierarchical predictive coding networks, and it extends previous EBM models in an important way by enabling a pure sampling-based inference procedure. The model combines multiple previous ideas into a single framework, and the approach feels natural. ", "weaknesses": "Although generation was impressive, this is also a proposal for Bayesian inference, and there were no serious benchmarks quantifying inference. Some thought should be put into evaluating known benchmarks for inference over hierarchical energy-based models. Many examples are possible, e.g. hierarchical gaussian mixture models, etc. \n\nClarity was a major issue. This was difficult to read. If this is going to be widely read, it should to be rewritten. I’ll describe specific problems below.\n\nFirst, there was not sufficient context in the surrounding literature. The introduction is rather short about EBMs. Otherwise it is difficult to know about this work’s specific contribution as many of these ideas have separately been introduced before (H-EBMs, Hamiltonian MC, etc.).  Also, for example, section 2 should summarize previous hierarchal EBMs, as currently it reads like the authors’ original proposal (which it is not). The related work section should be greatly expanded.\n\nThere were not also enough details about the experimental methods. I would certainly not be able to reproduce these figures from the manuscript. If there is not enough room, at least please put more details in the supplementary figures. \n\nFor clarity, I would recommend adding more detail about the biological circuit proposal. It is difficult for me to imagine the connectivity of the inhibitory neurons. (Do they need to project from one area to another, for example?) A figure with biology in mind would help. More details about biology should be added, as right now it is rather vague as to how this would be implemented in the brain. \n\nIn addition, many sentences state a controversial hypothesis as a known fact. These should be weakened considerably. For example,\nLine 46: “In this study, we show that our brain holds an intrinsic energy-based model”. ‘Show’ implies a certain conclusion. Better would be ‘hypothesize’ or ‘propose’.\nLine 70: “We demonstrate that our brain holds an intrinsic energy-based model”. Again, ‘demonstrate’ is conclusive. You are demonstrating that the brain *might* use such a model.\nLine 1: “The brain engages in probabilistic inference” This is just a hypothesis, although a popular one. Neural networks can appear to be Bayes-optimal even though they are not at all probabilistic internally. See Orhan, A. Emin, and Wei Ji Ma. \"Efficient probabilistic inference in generic neural networks trained with non-probabilistic feedback.\" Nature communications 8.1 (2017): 138.\n\n\n**Minor comments:**\n\nThe word ‘society’ is used when I think the authors want ‘community’.\n\nIn general the sentences are often awkwardly phrased. I recommend asking for editing help, perhaps using a large language model for grammar.\n\nSection 4. This seems to be Hamiltonian Monte Carlo. Is that correct? If so, this should be cited (e.g. Neal (2010)) as well as previous HMC in the brain proposals, like the Aitchison and Lengyel (2016) paper cited elsewhere. If it is incorrect, at least please cite SLD.\n\nLine 31: Many have argued that variational methods are biological plausible. See for example:\nD. Rezende and W. Gerstner, “Stochastic variational learning in recurrent spiking networks,” Frontiers in Computational Neuroscience, vol. 8, p. 38, 2014, doi: 10.3389/fncom.2014.00038.\n\nLine 84: “avoid complex calculations”. Could you be more specific?\nLine 121-124. This sentence is unclear and needs to be unpacked. What is amortized generation method, etc.\n", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We acknowledge the encouraging and valuable comments of the reviewer, and would like to address the concerns of the reviewer in details below.\n\n**Weaknesses:**\n\na) Quantifying Inference:\nThank for the suggestion of incorporating benchmarks into the inference section. This aligns with our future plan. Besides the benchmarks you mentioned, we intend to compare our results with \"neural prediction accuracy\" benchmarks [1]. This benchmark entails using our model's representations to predict real brain neuron representations.\n\n[1] Zhuang, Chengxu, et al. \"Unsupervised neural network models of the ventral visual stream.\" Proceedings of the National Academy of Sciences 118.3 (2021): e2014196118.\n\nb) Clarity:\nThank for the suggestions, in the revised version, we will\n\n1. provide a more comprehensive introduction to EBMs to ensure a clear understanding of our work.\n\n2. largely re-write the figure 6 caption to provide detailed descriptions for each subfigure.\n\n3. merge figures 1 and 2 in the revision and present the neural connection diagram.\n\n4. improve the writing and we will diligently address the parts highlighted by the reviewer, alongside other instances.\n\nc) Minor Comments:\n\nIn Section 4, the second-order Langevin sampling we used indeed bears similarity to Hamiltonian Monte Carlo, although they are not entirely the same. Both methods fall within the broader family of sampling techniques involving auxiliary variables to accelerate the sampling process [1]. Neal (2010) first introduced the concept of augmenting sampling speed using auxiliary variables, and Aitchison and Lengyel (2016) extended this to neuroscience by employing a group of inhibitory neurons as extra variables. In our work, we propose using adaptation current as an auxiliary variable. We will appropriately cite the references in the revised version.\n\nAs pointed out by the reviewer, the potential application of variational methods in the brain have been explored before. Some studies have also compared variational methods with sampling-based methods to determine their relative merits [2]. Sampling-based methods face convergence challenge, while variational methods involve both inference and generation stages during training, posing difficulties in multi-layered structure. In our current work, we primarily discussed the potential of sampling-based methods. Determining the brain's exact strategy requires further investigation, and we will duly address variational methods in the revised version.\n\nThe inference process of a generative model involves calculating the posterior distribution $p(z|x)$, where $p(z|x)=p(x,z)/p(x)$. Computing the denominator term, $p(x) = \\int p(x,z) dz$, becomes challenging when the dimensionality of $z$ is high. This integral can lead to complex calculations. Variational methods circumvent this by approximating $p(z|x)$ using a tractable distribution $q(z|x)$, thus avoiding the need to calculate $p(x)$. In sampling methods, like ours, we only require knowledge of $\\nabla_z \\ln p(z|x)$ to obtain samples that follow the distribution $p(z|x)$. Considering $\\nabla_z \\ln p(z|x)= \\nabla_z \\ln p(x,z)$, we also avoid the need to calculate $p(x)$.\n\n[1] Ma, Yi-An, Tianqi Chen, and Emily Fox. \"A complete recipe for stochastic gradient MCMC.\" Advances in neural information processing systems 28 (2015).\n\n[2] Grabska-Barwinska, Agnieszka, et al. \"Demixing odors-fast inference in olfaction.\" Advances in Neural Information Processing Systems 26 (2013).\n\n**Questions:**\n\na) Convergence speed depends on network size and depth. For our CIFAR10 dataset, with a network of 10 layers and a total of 150k neurons, typically 50-100 $\\tau_z$ steps are needed for high-quality image generation.\n\nb) Regarding Figure 2B, interneurons $\\varepsilon_l$ and principal neurons $x_l$, $x_{l+1}$ are connected. The $\\varepsilon_l-x_l$ connections follow one-to-one pairing, and $\\varepsilon_l-x_{l+1}$ connections are established via the weight matrix $\\theta_l$. The internal connections of $\\varepsilon_l$ are determined by $g'(u_l)$ as described in Eq.(9). The discussion of interneurons is addressed in the **rebuttal to all reviewers (part 2)**.\n\nc) In the current manuscript, it is not possible to discern the specific neuronal connections from the figures alone. One would need to refer to the corresponding dynamic equations Eq.(9-10) to understand the precise connections. We will address this shortcoming in the revision by providing an illustration of the neuronal connectivity diagram.\n\n\n**Limitations:**\n\n*“What are the least biologically plausible aspects? ”*\n\nAs you and many other reviewers have pointed out, the assumption $\\tau_u \\ll \\tau_z$ (line 128) is indeed a strong one. The rationale behind making this assumption is that during the inference process described by Eq. (8), we need to have real-time knowledge of the value of $A'(\\eta_l)$. However, we compute $A'(\\eta_l)$ through sampling (as in Eq. (9)), and sampling inherently introduces some temporal delay. To minimize this delay, we need to ensure that the timescale of sampling $A'(\\eta_l)$, denoted as $\\tau_u$, is much smaller than the timescale of inference dynamics $\\tau_z$. Nevertheless, we have also identified some empirical evidence that supports the validity of this assumption in the reply to all reviewers.\n\n*“What aspects need to be experimentally confirmed?”*\n\n1.Our model assumes SFA as the adaptation term to accelerate the sampling process (as discussed in Section 4). Suppose one can manipulate the SFA strength in the experiment, then we can test this assumption. \n\n2.Additionally, our model offers substantial flexibility, including the choices of $\\phi(x)$, $\\eta_l$, and $g(x)$. Currently, our selections primarily consider their impact on the model's learning capability (as discussed in Section 5.1), without necessarily aligning with real neural systems. Perhaps in the future, we could conduct experiments to validate the appropriateness of these choices.\n\n"}, {"review_id": "SQWXzhrURe", "reviewer": "Reviewer_nSfE", "summary": "In this article the authors introduce a new model, called the ‘’Hierarchical Exponential family Energy’ (HEE) that is biologically plausible and captures the dynamic of inference and learning in the brain. The model introduces multiple layers to decompose the EBM normalizing constant in a bio-plausible fashion, and leverages a neural adaption mechanism to make the sampling process compatible with ‘biological time’. In a series of experiments on synthetic dataset, Fashion MNIST and CIFAR10, the authors demonstrate the generation abilities of the HEE model. In addition, the authors demonstrate that the representation elicited by the HEE are similar to those observed in biological systems.", "questions": "* You did not include any comparison with other EBM for Fashion MNIST or the synthetic dataset. Is there a reason for that ? The synthetic dataset would have offered the perfect opportunity to conduct the comparison described previously .\n\n* One of the major difficulty in the IEBM [1], is the lack of stability of the training procedure (partly due to the mode traversing issue in the negative phase). Have you faced the same issue with the HEE ?\n\n* Line 200 : I don’t understand the sentence : \"The dropout technique is employed as a method to mimic the receptive field behavior found in biological systems. » Could you further explain ?\n\n* Have you tried some convolutional architecture ? Don’t you think it would improve the results ?\n\n* The difference between HEE-NL-A and HEE-NL is only the number of layers ? This should be clarified.\n\n* In Table 1, could you define LS and SLD ?\n\n* In Fig5, the samples seems be visually lower quality than those of the IEBM [1], but the FID and IS are better. Do you have an explanation ? In particular, in the samples you present, the textures seem to be rather uniform…\n\n* What is the impact of the number of layers in the quality of the generation ??\n\n\nTypo : \nLine 65 : EMB —> EBM \nLine 135: PCN is used but not defined yet.\n\n[1] Du, Yilun, and Igor Mordatch. \"Implicit generation and generalization in energy-based models.\" arXiv preprint arXiv:1903.08689 (2019).", "limitations": "Limitations have been properly addressed by the authors", "rating": 6, "confidence": 3, "soundness": 3, "presentation": 2, "contribution": 4, "strengths": "This article offers an interesting and valuable perspective on Energy Based Models (EBM) taking inspiration from neuroscience. The proposed solution to the intractable normalizing constant is original and elegant, and its link to predictive coding is interesting. Overall, there are plenty of good ideas in these article. \n", "weaknesses": "My biggest concerns are related to the experimental part:  (1) ablation studies that would validate the mathematical choice of the authors are missing , and (2) in general the experiments lack careful description. In addition, (3) the last part concerning the similarities between biological and HEE representation lack better comparison with other models. See more detailed comment below:\n\n\n1 ) One of the innovation of this article lies in the decomposition of the intractable partition function using multiple layers. Intuitively I do see the advantage of the proposed method compared to the bufferization of the amortization method used in standard EBM. But I was expected a comparative experiment between the 3 different methods to confirm the advantage of the proposed method. The comparison with the IEBM of [1] is not carefully enough conducted to conclude to the superiority of the proposed method : the parametrization seems to be different (CNN versus FullyConnected + number of parameters) and other improvement are introduced (e.g. neural adaptation) that might biased the comparison. In general I would suggest more ablation experiments that would test the benefits of all the HEE components separately.\n\n2) The networks are poorly described in the experimental part. What is the exact parametrization of the network (dimension of the fully connected network) ? What is the exact sampling procedure you use in the HEE : Do you first wait for a stable point in the sampling of the latent variable x_l before sampling the next one ? Or do you do one step of sampling for all the variables until you reach a global stable point ? Such implementation details might lead to very different equilibrium state and should be at least discussed. Also, how many sampling steps are you using ? Or are you forcing a stoping criterion (based on the prediction error ?) to decide when to stop the sampling process ?\n\n3) The experimental results in section 5.2 lacks comparison with other similar model. There is nothing showing the results you show are not  general to all networks (and not specific to HEE). Could you compare the HEE representations with those obtained for other EBM (or even other generative model) ? The idea would be to demonstrate that an EBM compatible (in terms of implementation) with biology would produce representation more aligned with biology.\n", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We acknowledge the encouraging and valuable comments of the reviewer, and would like to address the concerns of the reviewer in details below.\n\n**Weaknesses:**\n1. Please see **rebuttal to all reviewers (part 3)**.\n\n2. We recognize the current paper lacks detailed experiment descriptions, and we commit to add them in the revised manuscript. We will also incorporate a diagram illustrating the neuronal connections and network parameters.\n\n    The network parameters, denoted as $\\theta_l$, represent the connection weights between $x_{l+1}$ and $\\varepsilon_l$ (dashed arrows in Figure 2B). In the experiment on CIFAR-10 dataset , 80% of $\\theta_l$ elements are randomly set to zero permanently, effectively disconnecting 80% of connections between $x_{l+1}$ and $\\varepsilon_l$. We refer to this operation as \"dropout\" following the terminology in machine learning, but perhaps “sparse connectivity” is a better word.\n\n    Our sampling approach involves performing every step of sampling for all variables till a global stable point is reached. This approach aligns better with the biological plausibility, and we have yet explored hierarchical sampling. We have presented the details of the training and generation processes in Supplementary Material (section 3). During training, we stop sampling after every 300 $\\tau_x$. time steps, using the Euler method with a time step of $dt = 0.01\\tau_x$.\n\n3. We appreciate the suggestion and plan to incorporate these experiments and results in the revision. We are considering to add data from IEBM and PixelCNN [1]. Additionally, in future work, we intend to compare our model's results with the neural prediction accuracy benchmarks [2], by using our model's representations to predict the real brain neuron representations.\n\n    [1] Aaron Van Oord, Nal Kalchbrenner, and Koray Kavukcuoglu. Pixel recurrent neural networks. In ICML, 2016.\n\n    [2] Zhuang, Chengxu, et al. \"Unsupervised neural network models of the ventral visual stream.\" Proceedings of the National Academy of Sciences 118.3 (2021): e2014196118.\n\n**Questions:**\n\n1. Because our primary focus is on theoretical neuroscience, our work aims to present a biologically plausible learning method. The purpose of the experimental section is to demonstrate that our proposed model (a generative model) could potentially underlie the mechanism behind perception in the brain. The experiments conducted on the 2D synthetic and FashionMNIST datasets in Section 5.1 were primarily intended to showcase the learning capabilities of our proposed 'local learning' approach. We aimed to compare different sampling methods and identify models with favorable performance. Subsequently, the study of using the CIFAR10 dataset was motivated by its inclusion of natural images that closely resemble visual inputs received by the brain. This evaluation aimed to provide evidence that our model indeed learns the underlying distribution of these natural images. This serves as a foundation to support the discussion in Section 5.2, where we explored the model's representations and their alignment with real brain neuron representations. In the revised manuscript, we will include performance comparisons with IEBM on 2D synthetic datasets and the FashionMNIST dataset.\n\n2. While we have not trained IEBM, we can't pinpoint the exact reason for the issue it encountered. However, we did face instances of sudden gradient increasing and model instability during the training. After careful analysis, we identified that the uneven data distribution during learning led to inflated eigenvalues of the Jacobian matrix in the sampling dynamic Eq. (10). To mitigate this, we introduced a regularization term to the function $g(x)$, as described in Supplementary Material (section 3 and $g(x)$ in Table 1).\n\n3. We randomly set 80% of $\\theta_l$ elements to zero, effectively disconnecting 80% of connections between $x_{l+1}$ and $\\varepsilon_l$. This implies that during inference, each neuron in the second neural network layer only receives 20% of information from the first layer directly, and this setting extends to subsequent layers (each neuron in the third layer can receive information from the first layer equivalent to 1 - 0.8 * 0.8.). This concept closely resembles the receptive field concept in the visual system.\n\n4. We have not attempted this approach, as it would compromise biological plausibility (weight sharing is not biologically plausible). To better compare with IEBM, we plan to explore the CNN structure in future work.\n\n5. HEE-NL-A and HEE-NL differ in their sampling methods. HEE-NL employs Langevin sampling without introducing adaptation, while HEE-NL-A incorporates second-order Langevin dynamics for sampling.\n\n6. LS stands for Langevin sampling, and SLD for second-order Langevin dynamics. We will provide detailed explanations in Table 1 in the revised manuscript.\n\n7. Our IS (6.47) is indeed less than the 10-ensemble approach of IEBM (6.78), while utilizing FID for evaluating our results yield better. Simply speaking, IS measures the quality of generated images, while FID assesses the similarity between generated images and those within the dataset. Therefore, in term of scoring, our generated image quality might not be as high as IEBM's, but the similarity of our generated images to CIFAR10 is closer compared to IEBM. Visually, the perceived lower quality and uniform texture might be attributed to the small size of Figure 5. Additionally, the images not being vector graphics could contribute to this. In the revised version, we will address this concern by including larger generated images in the Supplementary Material. This will provide a better opportunity to observe the finer details of generated images. \n\n8. For FashionMNIST, we tried L=5, 10, 15, and 30. L=5 yielded the worst result, L=10 and L=15 exhibited similar performance, L=30 faced gradient explosion during training. Hence, we employed L=10 on the CIFAR10 dataset.\n"}, {"review_id": "6P6yh7WJb8", "reviewer": "Reviewer_tHDe", "summary": "This paper proposes a new model for sampling-based neural inference. The model posits that the brain attempts to match the marginal distribution $p(\\mathbf{x})$ of an internal generative model $p(\\mathbf{x}|\\mathbf{z})$ (based on a neural representation $\\mathbf{z}$) to the observed distribution of sensory inputs. In the model, sampling occurs via Langevin dynamics. The contributions of the paper are 1) use of a hierarchical exponential family model, 2) introduction of a second latent represenation $\\mathbf{u}$ (putatively identified with interneuron dynamics) responsible for estimating the demeaned residual $\\phi(\\mathbf{x}) - A'(\\boldsymbol{\\eta})$, and 3) use of a second-order Langevin sampling scheme. Experiments show that the resulting generative model can perform on toy data sets and there are further claims about relation to neural data in the experiments.\n\nThis is a potentially interesting idea that I found to be muddled by lack of justification of several modeling choices and unclear exposition, particularly in discussing the rationale and results of the experiments.\n", "questions": "- ll. 152-156: It's clear from (17) that $\\mathbf{v}$ represents a a noise term with temporal autocorrelation. In typical neural models, this is assumed semi-empirically as the result of temporal correlations among a large number of inputs to a given system. Is there an intuition the authors can give in this text for why it is better to identify (17) with spike frequency adaption as opposed to the existence of autocorrelated noise in synaptic inputs?\n- It's not very clear from the results presented what benefits hierarchy in processing brings. Can the authors show in some more convincing way which posterior features hierarchy specifically helps to capture?", "limitations": "- The model in the paper assumes a feedforward generative model rather than a recurrent model, as in most brain circuits. \n- The assumption $\\tau_u \\ll \\tau_z$ (line 128) is a strong one, and one presumes is justified by the relatively higher firing rates of cortical interneurons compared with principal neurons. But can this be justified from firing rates alone? Is there an argument that needs to be made about relative equilibration times of neuronal voltage dynamics for these two types of cells to substantiate this?\n- It's clear from the material in the supplement that some versions of the model cannot deal with multimodal posteriors, and even the second-order Langevin dynamics proposed will have difficulty mixing well in high-dimensional latent spaces.", "rating": 3, "confidence": 4, "soundness": 2, "presentation": 1, "contribution": 2, "strengths": "- The neural sampling framework has been an interesting line of inquiry over the previous decade, particularly in thinking about how the brain might implement Bayesian inference, and this work furthers that approach.\n- I found the use of a hierarchical model and the second latent representation to be interesting technical twists that are potentially powerful.\n- The paper attempts to align both the additional latents and the second-order Langevin dynamics with known features of neural physiology.", "weaknesses": "- The paper is written as if its main contribution will be to theoretical neuroscience, but the experiments present evaluations of the generative model on toy data sets, and the experiments underlying Figure 6 have an unclear rationale.\n- Overall, the paper has some deficits in presentation. For instance, section 2 is somewhat confusingly written. There is some parallel material that might be merged and better motivation given to the proposed choices. Why Langevin dynamics? Why not some other sampling method? Similarly, I didn't find the diagrams in Figure 1 particularly helpful in giving an intuition for the math. \n- Likewise, the rationale behind most of the experiments is not clearly explained in the text, particularly, those underlying Figure 6. Captions are too sparse to explicate what the results mean, and the discussion in text only barely explains the relevant neuroscience background. I have no idea what analysis was done for Figures 6A and 6B, nor why. Figures 5 and 6E are impossibly small.\n- ll. 120-23: This is a pretty compressed discussion of previously proposed approaches and perhaps difficult to follow for a non-expert reader.\n- More generally, identifications of model constructs with neuroscience findings ($\\mathbf{u}$ with interneuron activity, $\\mathbf{v}$ with spike-rate adaptation) are simply stated without any real justification. I'm not saying these choices are impossible to justify, but no real arguments that these are plausible assumptions are given. Again, if this is supposed to be a generative model loosely based on neuroscience, that can be fine, but if it's a theoretical neuroscience model, it's not clear these choices have been fully thought-through, and the experiments should ideally be focused on demonstrating that the model can perform like the brain in some perceptual task.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you very much for your thorough review and valuable suggestions. We hope the following explanations will help alleviate the concerns raised in the weaknesses, questions, and limitations.\n\n**Weaknesses:**\n\n1. As the reviewer correctly pointed out, our primary contribution is on theoretical neuroscience. Previously, Bayesian brain models have mainly focused on the inference process by using a pre-defined generative model and rarely considered how the generative model was learned. Here, we have taken a step further by explicitly modeling the learning process within the generative model framework. This inclusion may potentially give us insight into the learning mechanism in the brain's perceptual system. In Section 5.1, we employed several generative datasets to assess the learning capabilities of different energy models, and these evaluations help us to identify energy models of higher expression power that can learn complex data distributions. We used three datasets. Notably, CIFAR10 bears the closest resemblance to natural scenes the brain encounters. Therefore, we leveraged the representations learned from CIFAR10 to draw comparisons with neural representations in the brain in Section 5.2, and comparisons indicate that a HEE-based generative model may underlie the mechanism of the brain representing the external world information. This forms the rationale behind Section 5.2 and Figure 6.\n\n2. We appreciate the feedbacks on our presentation, and we will improve Section 2, Figure 1, and all other parts of the paper for better clarity. The reason for we choosing Langevin dynamics as an example to introduce our overall framework in Section 2 is that: Langevin dynamics has been used in previous works to successfully model the sampling-based neural inference (see e.g., [1-2]), which helps readers to understand our framework easily (we will add this underlying rationale in the revised manuscript). In Section 4, we also discussed alternative sampling methods the brain might employ, such as second-order Langevin dynamics. \n\n3. Thanks for the comments and suggestions. We will include the analysis details of Figure 6 and enlarge Figure 5 in the revised version. Figure 6 compares the learned representations of the HEE model from CIFAR10 (a natural image dataset) with the real brain representations in several facets, aiming to demonstrate the potential role of a generative framework in neural information representation. \n\n4. We are grateful for the suggestions, and we will elaborate and integrate those short discussions on previous methods isolatedly appeared in several parts of the paper (line 34-38, line 120-23, and line 255-261).\n\n5. As pointed out by the reviewer, our current work lacks a thorough exploration of the model's correspondence with the real neural system. Here, as the first step, our primary goal is to establish a framework, which addresses a particularly challenging problem in theoretical neuroscience, i.e., achieving simultaneous integration of biologically plausible local learning and sampling-based inference within a hierarchical network structure. But, as pointed out by the reviewer, we should at least give some justifications about the potential biological plausibility of model setting, which we will do definitely in the revised manuscript.   \n\n**Questions:**\n\n1. Thanks very much for the insightful suggestions. Indeed, interpreting eq.(17) as extracting the autocorrelation of inputs is biologically more reasonable. We will take this advice in the revised manuscript.  \n\n2. Previous researches have widely discussed about the advantages of hierarchical networks on information representation/learning. Compared to single or double-layer structures, deep networks can better approximate complex probability distributions p(x), enabling deeper layers to extract more abstract features [5], and these features further facilitate downstream tasks in the brain, such as decision making. In our work, Section 4 (lines 172-179) addresses the impact of network depth on convergence speed while keeping the total number of neurons fixed. Notably, an optimal network depth exists, as depicted in Figure 3D. Furthermore, our experiments demonstrate that as the number of layers increases, neural representations exhibit an enhanced linear discriminative capacity for object categories (Figure 6G), while neurons with orientation tuning decrease in number (Figure 6D). We will expand this discussion in the revised version.\n\n**Limitations:**\n\n1. Please note that although our generative model adheres to a Markovian process and can be perceived as a feedforward process by graphical model (Figure 2A & Eq.(7)), the neural implementation of our generative model actually involves recurrent connections between neurons, resembling a recurrent neural network (Figure 2B & Eq.(9-10)). Modeling brain circuits involves neural networks rather than graphical model. This is also true for other generative models [1-4].\n\n2. We appreciate the reviewer for offering two theoretical approaches that could potentially substantiate this hypothesis. Please see **rebuttal to all reviewers (part 2)** for further discussion.\n\n3. In our framework, we need suitable energy models that can express multimodal distributions (as this is essential for capturing the statistics of natural images). Therefore, in Section 5.1, we assessed the expressive powers of different energy models and selected the one which is capable of representing multimodal distributions for the followed study.\n\nWe are deeply grateful for the comments of the reviewer, which has provided significant insights for improving our work. We are committed to addressing these issues comprehensively in the revised version, and hope that the reviewer can raise the score accordingly.\n\n[1] Rodrigo Echeveste,Nature neuroscience,2020\n[2] Agnieszka Grabska-Barwinska, NeruIPS,2013\n[3] Nessler, B,Pfeiffer, PLoS Computational Biology,2013\n[4] Hennequin,NeurIPS,2014\n[5] DiCarlo,Neuron,2013\n"}, {"review_id": "8f7V4nvI8D", "reviewer": "Reviewer_moHy", "summary": "The paper introduces the hierarchical exponential family energy based model as a biologically plausible mechanism for the brain to interpret the external world. The authors first describe the learning/inference and the generation process of the model, the dynamics is local and biologically plausible through introducing a set of fast interneurons for the log-partition functions. They further show that adding adaptation mechanism can accelerate the sampling process. Finally, the authors present extensive numerical results to show that the model achieves good generation quality, and exhibit representations similar to those in the biological visual systems. ", "questions": "1. Eq4: It is a little confusing, before the second equality there is the KL divergence which is defined by Eq.2 with averaging over $z$ and $x$, whereas after the equality it is actually $p_{\\theta}(x,z)$ for a sample of $\\{x,z\\}$. Maybe it would be better to introduce separate notations for the random variables $x$ and $z$ (perhaps capital letters) and the realizations of $x$ and $z$.\n2. Line 99-100: If I understand correctly, both joint generation and marginal generation can generate samples following  $p_{\\theta}(x)$? Why say “However, in order to get the marginal distribution $p_{\\theta}(x)$” here?\n3. What could the fast interneurons $\\epsilon_l$ correspond to in the visual circuitry? \n4. Is there an intuitive explanation why the marginal method always seem to perform better than the joint method for generation? \n5. Fig 6D: As the proportion of both simple and complex cells decrease with the number of layers, do you see neurons with more complicated receptive fields (objects, etc.)? \n", "limitations": "Yes the authors have addressed the limitations.", "rating": 8, "confidence": 3, "soundness": 4, "presentation": 4, "contribution": 3, "strengths": "1. This paper presents a nice combination of theoretical analysis and extensive numerical results to demonstrate that their model exhibits multiple desired properties (including the acceleration effect of the adaptation mechanism, generation quality and representations). \n2. The authors show multiple interesting aspects of their models that are very relevant to the neuroscience community, including the acceleration effect of neural adaptation, and the similarity between the representations in biological visual system and the model, thereby providing new insights on the possible advantages/origins of well-known neural mechanisms/representations in the context of Bayesian inference in the brain. \n", "weaknesses": "The paper is very well written in general, and establishes nice connections between their findings and neural representations in the visual system. However, the comparison with neuroscience is limited to well-known experimental findings, it would be nice if the authors can add a discussion to highlight several possible experimental predictions of their model. ", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We acknowledge the encouraging and valuable comments of the reviewer, and would like to address the concerns of the reviewer in details below.\n\n**Weaknesses:**\n\nWe are grateful for the issues the reviewer highlighted. In the present study, as the first step, we aim to establish a framework, which addresses a particularly challenging problem in theoretical neuroscience, i.e., achieving simultaneous integration of biologically plausible local learning and sampling-based inference within a hierarchical network structure. In our forthcoming work, we will explore in details how model predictions can be experimentally validated. Recently, we have also developed a keen interest in inputs with temporal information and are conducting relevant experiments on stimuli with serial dependence [1]. We will incorporate some discussions about the latest developments in the revised manuscript.\n\n[1] Fischer, Jason, and David Whitney. \"Serial dependence in visual perception.\" Nature neuroscience 17.5 (2014): 738-743.\n\n**Questions:**\n1. Thanks for the suggestion. In the revised version, we will distinguish between random variables and their samples. For example, Eq. (2) will be revised as follows: $\\nabla D_{KL}\\left[p_{true}(\\mathbf{x})\\parallel p_{\\theta}(\\mathbf{x})\\right]=-E_{\\tilde{\\mathbf{x}}\\sim p_{true}(\\mathbf{x})}E_{\\tilde{\\mathbf{z}}\\sim p_\\theta(\\mathbf{z}|\\tilde{\\mathbf{x}})}\\left[\\nabla_\\theta\\ln p_\\theta({\\tilde{\\mathbf{x}},\\tilde{\\mathbf{z}}})\\right]$, where $\\mathbf{x},\\mathbf{z}$ are random vairables and  $\\tilde{\\mathbf{x}}, \\tilde{\\mathbf{z}}$ are samples\n\n2. Indeed, both joint generation and marginal generation can produce samples following $p_{\\theta}(\\mathbf{x})$. We recognize that our expression was unclear. To clarify, joint generation involves simultaneous sampling in the x-space and z-space, leading to x's marginal distribution following $p_{\\theta}(\\mathbf{x})$; whereas, marginal generation first samples in the z-space and then generates x according to $p_\\theta(\\mathbf{x}|\\mathbf{z})$, which also results in the x's marginal distribution conforming to $p_{\\theta}(\\mathbf{x})$. We will clarify this in the revised manuscript.\n\n3. Please see **rebuttal to all reviewers (part 2)**.\n\n4. Joint generation involves simultaneous sampling in both x-space and z-space according to $p_\\theta(x,z)$. On the other hand, marginal generation first samples from z-space according to $p_\\theta(z)$, and after obtaining $z$ samples, proceeds to sample from x-space following $p_\\theta(x|z)$. Considering the time complexities of $O(m)$ for x-space and $O(n)$ for z-space sampling, joint generation has a time complexity of $O(m*n)$, while marginal generation's complexity is $O(m+n)$. This implies that marginal generation is more likely to yield high-quality samples as its sampling space is smaller.\n\n5. We have investigated not only neurons of orientation tuning, but also neurons tuned to high-order features [1] and neurons responsive to image categories. However, the proportion of these neurons was extremely low, constituting less than 1%. Therefore, we chose not to report these results. \n\n    Although only a small number of neurons exhibit selectivity for categories at the single-neuron level, information about categories can still be linearly decoded from the population response. We used neural representations in each layer to train a linear SVM to discriminate the categories of CIFAR-10 images. We found that the classification accuracy increases with the layer hierarchy (Figure 6G), with the final layer achieving an accuracy exceeding 60%, significantly surpassing the random baseline accuracy of 10%. \n\n[1] Julesz, Bela. \"Textons, the elements of texture perception, and their interactions.\" Nature 290.5802 (1981): 91-97.\n"}]

Author rebuttal:
**Reply to all reviewers**

We greatly appreciate the thorough review of our paper by all four reviewers. The positive feedbacks in the reviews have effectively highlighted the contributions of our work. We will provide a summary of this aspect in Part 1. Additionally, we will address the common questions raised about our work in Parts 2 and 3.

**Part 1**

Previously, Bayesian brain models primarily emphasized the inference process, utilizing pre-defined generative models, and rarely delved into how the generative model was acquired. In this work, we have extended beyond this by explicitly incorporating the learning process within the hierarchical generative model framework. This addition could potentially offer insights into the learning mechanism within the brain's perceptual system.

Furthermore, our brain-inspired energy-based model (EBM) also presents a technique for estimating the partition function in EBMs, which is a challenging problem within the machine learning community.

**Part 2**

As reviewers zTHc, tHDe and moHy pointed out, $\tau_u\ll\tau_z$ (line 128) is a strong assumption. In the current paper, we do acknowledge the absence of a discussion on this assumption and interneurons. Therefore, we are providing additional discussion here and will incorporate this aspect into the revised version.

There are various types of interneurons that target on pyramidal cells, comprising approximately 10-20% of the overall neuron population in the cerebral cortex [1]. The interneurons in the HEE model bear the closest resemblance to the Large Basket Cell or Nest Basket Cell [2], which collectively constitute around 50% of interneurons. Their electrophysiological characteristics include fast spiking, non-accommodating, and non-adapting behaviors. These interneurons have also been identified in the visual cortex of ferrets [3], displaying short-duration action potentials (approximately 0.5 ms at half height). This suggests that these neurons have shorter time constants compared to pyramidal cells. 

[1] Therese Riedemann. Diversity and function of somatostatin-expressing interneurons in the cerebral cortex. International Journal ofMolecular Sciences, 20(12):2952, 2019.

[2] Markram, Henry, et al. "Reconstruction and simulation of neocortical microcircuitry." Cell 163.2 (2015): 456-492.

[3] Descalzo, Vanessa F., et al. "Slow adaptation in fast-spiking neurons of visual cortex." Journal of neurophysiology 93.2 (2005): 1111-1118.

**Part 3**

As the reviewer nSfE pointed out, it is better to compare our model with other methods on the partition function estimation. However, we find that for the joint distribution in our model (belonging to the hierarchical exponential-family), implementing amortized generation and implicit generation are challenging. The difficulty arises from the fact that, in our model, besides requiring the partition function during training (Eq. 11), we also need the partition function during sampling (Eq. 8). In contrast, amortized generation and implicit generation do not consider partition function during the sampling process. Thus, if comparing our model with the other two methods, the joint distribution employed by them will be different from the one used in our model, and this makes comparison not in the same condition. We need carefully design experiments to address this issue.

As the reviewer pointed out, our comparison with IEBM lacks controlled variables. In term of connectivity, IEBM employs a CNN architecture, while our model employs a full connection structure. While our model could also be adapted to a CNN-like structure, the operation of weight sharing in CNN is not biologically plausible, as it requires the changes of synapses over a large space are coordinated. Therefore, in our model implementation, we chose a fully connected architecture instead of CNN to maintain better biological realism. Regarding the sampling method, IEBM employs Langevin sampling, while our approach employs Langevin sampling + adaptation. This technique, while not identical to Hamiltonian Monte Carlo (HMC), falls within the same category of sampling methods [1]. Notably, IEBM also discusses the impact of HMC and highlights the challenge in controlling the number of leapfrog simulations during training. In term of parameter count, IEBM utilizes around 5 million parameters, while our model employs 4 million. Additionally, in term of loss function, our approach directly employs log-likelihood, while IEBM incorporates an additional regularization term. Moving forward, in order to enhance the comparison with IEBM, we are considering a strategy that temporarily sets aside biological constraints, rather we will focus solely on experimental outcomes by adopting a CNN architecture, employing Langevin sampling, and ensuring parity in parameter count. By this, we can conduct a direct and informative comparison with IEBM.

[1] Ma, Y. A., Chen, T., & Fox, E. B. (2015). A complete recipe for stochastic gradient MCMC. Advances in Neural Information Processing Systems, 2015-Janua, 2917–2925.


Meta-review: This is a borderline case. Three reviewers voting for acceptance (with scores 5,6,8) , one for rejection (score 3). The main arguments for rejection are that the paper needs to be revised to more fully situate this paper in a theoretical neuroscience context. Reviewers however beliefe that a thoroughly revised draft can address all of the outstanding issues and that the theoretical novelty proposed in this article is an interesting contribution. I encourage the authors to incorporate all the promised changes into the revised version of the manuscript.
Decision: Accept (poster)