The Bayesian sampling in a canonical recurrent circuit with a diversity of inhibitory interneurons

We appreciate the reviewer's feedback and positive comments about the math rigor and biological plausibility of our work.

Two major limitations of the sampling circuit proposed in the paper are its low-dimensionality and uniform prior. Overall, the authors only demonstrate sampling of 1-D Gaussian distributions, due to the linear drift. While I understand that it is the norm of papers in neural sampling theory, recent works have shown that recurrent circuits are capable of sampling from much more complicated distributions [1,2]. The authors should discuss how the proposed circuits can scale in more detail.

[1] Lyo, Benjamin and Cristina Savin. “Complex priors and flexible inference in recurrent circuits with dendritic nonlinearities.” bioRxiv (2023): n. pag. [2] Chen, Shirui, et al. "Expressive probabilistic sampling in recurrent neural networks." Advances in Neural Information Processing Systems 36 (2024).

Thank you, we very much appreciate you bringing this issue. We will include the discussion of recent papers that sampled more complicated posteriors in the revised manuscript, and do some comparison.

In addition, the proposed circuit with SOM neurons can sample multi-modal and high-dim posteriors. Please see the results in the rebuttal PDF and our response in Global Rebuttal. We are happy to include the results to strengthen our paper in the revised version with one more page, if you encourage us to do that. The reason for not including them is because of the page limit and the completeness of the manuscript, where we sacrificed the inference task complexity to present the reasoning of building the model and math analysis.

Compared with the two papers mentioned by the reviewers, our study presented a simpler inference task, compensated by the analytical results of identifying sampling algorithms in the circuit, and biological plausibility (diverse interneurons, reproducing experimentally observed tuning curve). Our rebuttal PDF shows the proposed circuit model can also sample multi-modal and high-dim posteriors. We wish the reviewer could re-evaluate the significance of our model framework and analytical methodology.

Minor: Figure 3 caption does not include panel F.

Apologies for the typo. We will add it to the revised manuscript.

Minor: Although it is straightforward for the linear case (S36), for the mixture of general Langevin and Hamiltonian dynamics (S34 & S35), it is worthwhile to include references or show that the marginal distribution of z is indeed p(z) (is that true in the general case?).

The mixture of Langevin and Hamiltonian (Eqs. S34 & S35) was used by Ref. 14 and 50, and it is guaranteed to sample the correct $p(z)$. We will emphasize this property in the revised manuscript.

I am curious about how the authors would construct a distributed version of the proposed circuit that can sample from high-dimensional stimulus posterior.

Please see our results of extending the proposed circuit model to sample high-dim posteriors in the Global Rebuttal and attached Rebuttal PDF.

To sample high-dim posteriors, we can extend the current model into many coupled networks with each the same as the one in the current manuscript. Each network receives a neural input generated from a latent 1D stimulus, and samples the corresponding 1D latent stimulus. As a whole, the coupled networks distributivity sample high-dim posteriors, where the dimension of the posterior is determined by the number of networks in the coupled circuits. The recurrent weights between subnetworks store a high-dim correlation prior of stimuli. In this way, the coupled networks are similar to the network presented in Zhang et al., Nat. Comms. 2023 (see its Fig. 6), but each subnetwork in Zhang's study didn't consider SOM neurons.

Therefore, when we mix the above Eq. (3) with another Langevin sampling dynamics (e.g., Eq. S33 in the Supplementary), the mixed Langevin and Hamiltonian dynamics will still sample the same distribution... And then the dynamics can also sample distributions apart from Gaussian distributions.

This is exactly what I think is sketchy here, I don't think it holds true in general. In particular, your SDE is in z and y, but $\nabla \log p(z)$ can be nonlinear in $z$. Can you elaborate?

Thanks for giving us a chance to elaborate on our contribution.

The theoretical/analytical results we meant is dynamical system analysis of the nonlinear neural dynamics, e.g., perturbative analysis and bifurcation theory, whose details are included in Supp. Sec. 3. Without these analyses, we cannot analytically obtain the sampling dynamics in the low-dim stimulus manifold (Eqs. 11-12) embedded into the high-dim neural dynamics. We agree that Chen 23, Benjamin 23, and Echeveste 20 have elegant theoretical results, but they are not the dynamical system analyses performed in the present study. For example, our theoretical analysis corresponds to theoretically analyzing the recurrent dynamics in Eq. 8 in Chen 23 to find its attractor states, eigen-spectrum analysis, dimensionality-reduction etc.

What do we gain from neural dynamical system analysis?

The dynamical system analysis (Supp. Sec. 3) especially the analytical results (Eqs. 11-12) significantly gain our understanding of fine structure about sampling neural dynamics. For example, our study's nonlinear circuit dynamics is similar to Eq. 8 in Chen 23 in principle, both of which are Langevin dynamics. After our dynamical system analysis, we find the sampling dynamics embedded in the stimulus manifold $z$ can vary from Langevin to Hamiltonian sampling, although the neural dynamics are still Langevin. In addition, the proposed circuit with Langevin neural dynamics has the potential to implement Levy sampling in the stimulus manifold, if we increase the inhibitory feedback weight from SOM to E neurons further and also inject noises into SOM's dynamics, with a similar math mechanism as Dong, NeurIPS 2021 and Sakaguchi, 2021, J. Physical Society of Japan. Basically, we feel the nonlinear neural dynamics has a rich repertoire to implement various sampling algorithms by using dynamical bifurcations, i.e., the proposed neural dynamics are located at three different bifurcation regimes in realizing Langevin, Hamiltonian, and Levy sampling. These fine dynamical mechanisms to implement neural sampling cannot be revealed without analytical results of dynamical system analysis (the results Eqs. 11-12 derived from Eqs. 1-4) as presented in the current study.

The contributions and differences compared with earlier studies on the computational side.

1. Previous studies (Chen 23, Benjamin 23, and Echeveste 20) all considered neural dynamics directly sample in the neural response $r$ space, while the neural dynamics in the present study considers the sampling in a low-dim stimulus manifold $z$ embedded in the neural response $r$ space, similar to Dong 22.

2. The proposed circuit with fixed weights can sample all posteriors $p(z|r_F)$ under different observation $r_F$ which can have different uncertainties (Fig. 2G, 3F; Eq. 21). This is realized by embedding the low-dim sampling in the high-dim neural dynamics and the multiplicative variability in neural dynamics (Eq. 1, last term), where the overall neural firing rate in a population can carry the precision of distributions, i.e., the drift coefficient $g_{XY} \propto w_{XY} R_Y$ in Eqs .(11-12) is proportional to population firing rate $R_Y$ in line 188 . The flexibility of sampling a family of posteriors by neural dynamics with fixed parameters is a more stringent requirement (see reasons in Global Rebuttal), while it is unknown whether Chen 23 and Benjamin 23 could achieve this. Echeveste 20 had this numerical result (its Fig. 2) but it did not perform dynamical system analysis as in the present study so its underlying dynamical system mechanism is not clear. Dong 22 considered a similar network model to ours, while we are concerned that Dong's network cannot achieve this either, because the variance of injected variability (\sigma_V in its Eq. 13) needs to be re-adjusted based on instantaneous feedforward input (the 5th line after Eq. 20 in Dong 22).

$\dot z_t = \lambda\nabla U(z) + 2\lambda \epsilon_t$ Its equilibrium distribution is $p(z)\propto \exp[U(z)]$.

This is incorrect. $p(z)\propto \exp[U(z)/(2\lambda)]$. So even in your simplified version of the proof, the mixture does not sample from the correct distribution. EDIT: I realize that this must be a typo of the authors. But the following paragraph still describes my doubt.

Furthermore, I am also not sure how valid the approximation is, your S34 and S35 sure do not mix the Langevin and Hamiltonian in this way. S34 and S35 is a 2D SDE, and when you write down the Fokker-Planck equation, it is not obvious at all what is the equilibrium distribution.

I am actually okay with it if the authors cannot provide rigorous proof. But I would be very interested if the authors know any.

I also recommend the authors think more about how the recurrent weights between subnetworks store not only a high-dim correlation prior of stimuli but also other high-order statistics in order to represent high-dim complex posterior.

Given the other additions that the author contributed, I am happy to recommend acceptance.

Invariant equilibrium distribution by mixing Langevin and Hamiltonian dynamics

We performed some theoretical analysis to show the mixture of Langevin and Hamiltonian dynamics will not change the equilibrium distribution. We start by copying the equations of Hamiltonian dynamics with friction (Eq. S13-14) here,

$ \tau_z \dot{z} = y; \\ \tau_y \dot{y} = -\nabla U(z) - \beta y + (2\beta \tau_z)^{1/2} \eta_t, \quad\quad (1) $

where $U(z) = -\log p(z)$ is the logarithm of the distribution being sampled.

To intuitively and quickly show the equilibrium distribution of $p(z)$, we consider a case where the time scale of $y$ (determined by the time constant $\tau_y/ \beta$) is much faster than the time scale of $z$, i.e., $y$ changes much faster than $z$. And then we could approximately treat the $z$ is fixed during the time scale of $y$, and then we utilize time separation to analyze the two dynamics separately (this strategy was also used in the paragraph around Eqs. 11-12 in Ref. 50). The equilibrium distribution of $y$ can be solved as $p(y) = \mathcal{N}[y|-\beta^{-1}\nabla U(z), \tau_z/\tau_y]\quad\quad (2)$ Next, we consider the change of $z$. At the time scale of $z$, we can approximately treat the $y$ as a sample from its equilibrium distribution, i.e., $y\sim p(y)$, and write $y$ as $y = - \beta^{-1}\nabla U(z) + (\tau_z/\tau_y)^{1/2} \xi_t$ Substituting the above equation back to the Eq. (1) shown above,

$ \dot{z} = - \tau_z^{-1}\beta^{-1} \nabla U(z) + (\tau_z \tau_y)^{-1/2} \xi_t \quad\quad (3) $

As long as $\beta = 2\tau_y$, the above $z$ dynamics embedded in the Hamiltonian dynamics can be treated by performing Langevin sampling dynamics to sample the distribution $p(z) \propto \exp[-U(z)]$.

Then consider another Langevin dynamics to sample $p(z) \propto \exp[-U(z)]$(e.g., Eq. S33 in the Supplementary), $\dot{z} = [-\tau_L^{-1} \nabla U(z) + (\tau_L/2)^{-1/2} \eta_t]$ and mix it with the above Eq. (3), $\dot{z} = [-\tau_z^{-1}\beta^{-1} \nabla U(z) + (\tau_z \tau_y)^{-1/2} \xi_t] + [-\tau_L^{-1} \nabla U(z) + (\tau_L/2)^{-1/2} \eta_t] = -[\tau_z^{-1}\beta^{-1} + \tau_L^{-1}] \nabla U(z) + [(\tau_z \tau_y)^{-1} + (\tau_L/2)^{-1}]^{1/2} \epsilon_t$

Denoting the drift coefficient $\tau_z^{-1}\beta^{-1} + \tau_L^{-1} \equiv \lambda$ and using the above condition $\beta = 2\tau_y$, we can find the diffusion coefficient above is $\sqrt{(\tau_z \tau_y)^{-1} + (\tau_L/2)^{-1}} = \sqrt{2\lambda}$. Therefore the $z$ dynamics in the mixed Langevin and Hamiltonian dynamics is $\dot{z}_t = \lambda \nabla U(z) + \sqrt{2\lambda} \epsilon_t$ Its equilibrium distribution is $p(z)\propto \exp[-U(z)]$. Hence the mixture of Langevin and Hamiltonian dynamics can also sample the same distributions $p(z)$. Since $U(z)$ is general in the above derivation, the mixed dynamics can sample other types of distributions than Gaussian.

To sample non-Gaussian distributions, we need to modify the recurrent connection profile (Eq. 2). Please find our reply to Reviewer K6Di (in the section "non-Gaussian stimuli").

PS: sorry that our wording confused you. You are right that an SDE can only sample a Gaussian $p(z)$ if its drift term is linear with $z$, whereas the linear in our last reply meant the drift term is linear over $\nabla \log p(z)$.

Sorry for our typos. We have carefully checked our last reply and have corrected math typos, including the diffusion term in the last equation (should be $\sqrt{2\lambda} \epsilon_t$), and redefine the $U(z) = - \log p(z)$ in the Eq. (1) in the last reply by adding a minus sign to make it consistent with conventional definitions.

Calculating the equilibrium distribution $p(z)$ in the mixed dynamics via Fokker-Planck approach

Here we will provide a rigorous analysis showing the Eqs. (S34-S35) can correctly sample the desired distribution $p(z)$. Due to limited space, we rely on derivations used in Chen, ICML 2014 to explain the math mechanism. We start by defining the (equilibrium) distribution to be sampled and the corresponding Hamiltonian (Eq. S8), $\pi(z,y) = \exp[-H(z,y)]; \quad H(z,y) = U(z) + K(y) = -\ln p(z) + y^2/2$ And the Hamiltonian sampling with friction (Eqs. S13-S14) is (removing the last two terms of the Langevin step in Eq. S34)

$$\tau_z \dot{z} = y; ~~ \tau_y \dot{y} = - \nabla U(z) - \beta y + (2\beta \tau_z)^{1/2} \eta_t \quad\quad (R1)$$

which is analogous to the Eq. (9) in Chen, 2014. To emphasize the main mechanism, we consider a simple case where $\tau_y = \tau_z = 1$. For general cases of $\tau_y$ and $\tau_z$, we can set $\beta$ appropriately to sample the $p(z)$ (Eq. 20 in manuscript). Denoting by the vector $Z = (z,y)^\top$, we could convert the Eq. (R1) into the matrix notation,

$$\frac{dZ}{dt} = - \begin{pmatrix} 0 & -1 \\\\ 1 & \beta \end{pmatrix} \begin{pmatrix} \nabla U(z) \\\\ y \end{pmatrix} + \sqrt{2} \begin{pmatrix} 0 & 0 \\\\ 0 & \beta^{1/2} \end{pmatrix} {\boldsymbol \eta}_t = - (D + G) \nabla H(Z) dt + (2D)^{1/2} {\boldsymbol \eta}_t \quad \quad (R2)$$

where $G = \begin{pmatrix} 0 & -I \\\\ I & 0 \end{pmatrix}$ is an anti-symmetric matrix, $D = \begin{pmatrix} 0 & 0 \\\\ 0 & \beta \end{pmatrix}$, and $\nabla H(Z) = [\partial_z H(z,y), \partial_y H(z,y)]^\top = (\nabla U(z), y)^\top$. Now, the Eq. (R2) is the same as the unnumbered equation in Theorem 3.2 in Chen 2014. Similar to the Eqs. (20) and (21) in the Supplementary in Chen 2014, the Fokker-Planck equation of the above dynamics is

$$\partial_t p_t(Z) = \nabla^\top \{ (D+G)[\nabla H(Z)p_t(Z)]\} + \nabla^\top[D\nabla \partial_t p(Z)] \quad \quad (R3)$$

Then utilizing the math property of the anti-symmetric $G$, we have $\nabla^\top [G \nabla p_t(Z)] = - \partial_z \partial_y p_t(z,y) + \partial_z \partial_y p_t(z,y) = 0$, then inserting this zero term into the last term of Eq. (R3),

$$\nabla^\top[D\nabla \partial_t p(Z)] = \nabla^\top[(D+G)\nabla \partial_t p(Z)]$$

Substituting the above equation back to Eq. (R3), and taking the common term $(D+G)$ out of the bracket in the RHS, the Eq. (R3) can be converted into

$$\partial_t p_t(Z) = \nabla^\top \{ (D+G)[\nabla H(Z)p_t(Z) +\nabla \partial_t p(Z)] \} \quad \quad (R4)$$

It can be easily checked that $\pi(z,y) \propto \exp[-H(z,y)]$ is indeed the equilibrium distribution of the Eq. (R2) by calculating $\nabla H(Z) \exp[-H(z,y)] +\nabla \exp[-H(z,y)] = 0$ Upon this point, we complete the reasoning showing the Hamiltonian dynamics with friction can sample $\pi(z) \propto \exp[-U(z)]$. And above is our restatement in explaining the math derivations underlying Theorem 3.2 in Chen 2014.

Next, we mix the Hamiltonian dynamics (Eq. R2) with the Langevin dynamics of $z$ (as we did in Eqs. S34-S35), then the Eq. (R2) becomes

$$\frac{dZ}{dt} = - \begin{pmatrix} \tau_L^{-1} & -1 \\\\ 1 & \beta \end{pmatrix} \begin{pmatrix} \nabla U(z) \\\\ y \end{pmatrix} + \sqrt{2} \begin{pmatrix} \tau_L^{-1/2} & 0 \\\\ 0 & \beta^{1/2} \end{pmatrix} {\eta}_t = - (D' + G) \nabla H(Z) dt + (2D')^{1/2} \eta_t \quad \quad (R5)$$

where the $D' = \begin{pmatrix} \tau_L^{-1} & 0 \\\\ 0 & \beta \end{pmatrix}$. Comparing Eq. (R5) with Eq. (R2), we find mixing the Langevin dynamics only changes the notation of matrix $D$ (appearing in both drift and diffusion terms) into a new matrix notation $D'$. Therefore the corresponding Fokker-Planck equation of the Eq. (R5) will be similar with Eq. (R4) but just replacing the $D$ by $D'$

$$\partial_t p_t(Z) = \nabla^\top \{ (D'+G)[\nabla H(Z)p_t(Z) +\nabla \partial_t p(Z)] \}$$

And then the equilibrium distribution of the mixed dynamics (Eq. R5) is also $\pi(z,y)$, i.e., mixing the Langevin dynamics with Hamiltonian dynamics (Eqs. S34-35) doesn't change the equilibrium distribution $p(z)$.

We plan to incorporate a more detailed version of the above math derivation in a new Section 2.4 in our Supplementary to strengthen our theoretical analysis.

We appreciate the reviewer's feedback and positive comments about math analysis, and the insight of SOM neuron's role in inference presented in our manuscript

The biggest weakness of the paper...

Thanks for your suggestions; we will revise our manuscript accordingly by adding more content about intuitive descriptions. For example, we can add text related to Fig. 3C to explain why the structured inhibition from SOM neurons could speed up sampling by decreasing the temporal correlations of samples. Intuitively, Hamiltonian sampling comes with oscillation, and meanwhile, the structured inhibition from SOM will bring oscillations of stimulus samples computed by E neurons.

The authors emphasis that their results...

Thank you for your suggestion. We will weaken our wording in the revision. We will just use the wording "canonical circuit" in describing the Fig. 1A and call our proposed model a circuit with diverse types of interneurons. In addition, we can also remove the "canonical" in the title.

Further, authors use the assumption that PV neurons...

Thanks for the suggestion and we agree that the tuning strength of PV is under debate. We will revise our text by emphasizing the untuned PV interneurons in the model is an assumption. This simplifies our theoretical analysis, while reserving a certain amount of biological plausibility, via multiplicative modulation on E tunings (Fig. 1G). Specifically, we will weaken our text in lines 104-105 into
"Some experiments found the stimulus orientation weakly modulates the PV neurons, hence, for simplicity, we assume PV neurons in the model are not tuned to stimulus features, providing only global unstructured inhibition to E neurons for stability (see Discussion for the complexity of PV neurons' tunings)"
In Discussion we will cite Runyan 2010 and Moore 2013 to discuss the experimental discovery of some PV neurons with sharp tuning, remind readers that the tuning of PV neurons is debatable, and clearly state that global inhibition from PV is an assumption considered in the study.

Unless I misunderstood something....

Due to the page limit, our submitted manuscript only presents the 1D posterior sampling. The proposed circuit with SOM neurons can sample multi-modal and high-dim posteriors. Please find the results in the rebuttal PDF and our response in Global Rebuttal. We are happy to include the results to strengthen our paper in the revised version with one more page.

Does Bayesian sampling...

Yes, the current manuscript only considers a static stimulus and thus the stimulus samples embedded in the spatiotemporal neuronal activity will converge to the 1D stimulus posterior density, rather the just the point value of the stimulus.

The proposed circuit model, with or without SOM neurons, can be used to infer the changing stimulus $z_t$, i.e., inferring the instantaneous posterior $p(z_t| r_t ... r_1)$ based on all history of neural inputs $(r_t, ..., r_1)$ in a hidden Markov model. Then the variance of the transition prior probability $p(z_t|z_{t-1})$ will be stored in the peak recurrent weight between E neurons, $w_{EE}$ (Eq. 2), decreasing the variance of $p(z_t|z_{t-1})$ (we have analytical results about this). This is because a larger $w_{EE}$ will incur a stronger bump response of E neurons (dark blue in Fig. 1E) which moves slower over the stimulus space, corresponding to a transition probability in MCMC with smaller variance.

Can the model deal with non-Gaussian stimuli?

Yes, please see the attached Rebuttal PDF where the circuit can sample bimodal posteriors. In principle, the proposed framework can sample exponential family distributions (with a similar mechanism to Ma et al., Nat. Neurosci. 2006 and Zhang et al., Nat. Comms. 2023) or a mixture of exponential family distributions, where the spatial profile of recurrent weights (Eq. 2) acts as the natural parameter of distributions and determines the type of sampling distribution in the circuit. The Gaussian distribution comes from the Gaussian profile of recurrent weights (Eq. 2). If we change the recurrent weight profile into other shapes, e.g., von Mises function, the circuit model can sample von Mises distributions.

To sample distributions out of the exponential family, we may consider passing the responses of the current circuit model into a feedforward network which can rescale the Gaussian distribution into arbitrary distributions. This mechanism was used in a recent study, Chen et al., NeurIPS 2023 mentioned by Reviewer RifV.

Tuning curves of neurons of the same cell type are homogeneous...

Yes, the homogeneity simplifies the math analysis significantly, without altering our conclusion. Real cortical circuits have heterogeneity among neurons, and the homogeneity assumption has been widely used in many computational neuroscience models as a justifiable simplification, especially in the work of continuous attractor networks, e.g., Ben-Yi Shai 1995; Wu, Neural Computation 2008; Khona, Nat. Rev. Neurosci., 2022, etc.

Heterogeneous neuronal tunings could be realized by introducing randomness (zero mean with certain variance) into the recurrent connection matrix (Eq. 2), which has also been widely used in (chaotic) Excitation and Inhibition (E/I) balanced networks (Vreeswijk, Science 1996; Neural Computation 1998; Litwin-Kumar, Nat. Neurosci., 2012; Rosenbaum, Nat. Neurosci., 2017, etc).

A potential function of heterogeneity from random recurrent weights is that this puts the spiking networks into the chaotic regime where the network internally generates Poisson variability (last term in Eq. 1). We think the injected multiplicative variability in our rate-based network (Eq. 1) captures the chaotic Poisson variability from spiking network dynamics. Then, stronger random recurrent weights will induce large heterogeneity, corresponding to a larger Fano factor $F_E$ of the injected variability in Eq. 1.

Thanks for the reviewer's appreciation of our math analysis and the biological plausibility of the model.

The paper is quite dense, often unclear for a machine-learning audience, and ultimately only addresses a uniform-Gaussian joint density that cannot capture natural stimuli or behavior.

We will improve our manuscript to make it more friendly to machine learning audiences. We are sorry that the detailed justification of the biological plausibility of the circuit model and the theoretical analysis might appear less interesting to machine learning audiences, whereas they are probably necessary for a complete theoretical neuroscience study.

The proposed circuit with SOM neurons can sample multi-modal and high-dim posteriors. Please find the results in the rebuttal PDF and our response in Global Rebuttal. We are happy to include the results to strengthen our paper in the revised version with one more page, if you encourage us to do that. The reason for not including them is because of the page limit and the completeness of the manuscript, where we sacrificed the inference task complexity to present the reasoning of building the model and math analysis.

Where do the authors propose that their circuit ought to be found, neuroanatomically? How does it fit into the laminar structure proposed by existing cortical microcircuit theories?

The micro-circuit structure (Fig. 2A) is canonical and ubiquitous in all cortical brain regions in mammals. In recent years, most experiments utilized mice's visual cortex to investigate this canonical micro-circuit, because we have high-throughput recording and perturbing tools on mice.

In terms of laminar layers, the canonical micro-circuit (Fig. 2A) typically exists in both layers 2/3 (superficial layer) and layer 5 (deep layer), although the connection strength among neurons varies across laminar layers and across brain regions. The Fig. 7 in Ref. [24] in the submitted manuscript is a good summary of this micro-circuit structure in both superficial and deep layers.

See weaknesses. The paper's chief limitation is considering an overly simplified generative model based on previous work on probabilistic population codes, neglecting the observation that priors can be decoded (under supervised learning settings and at above-chance rates, etc.) almost ubiquitously from cortical recordings and are modulated by top-down feedback activity.

We agree that the derived generative model considered in our theoretical neuroscience study is simple compared with machine learning studies. It's worth noting the difference in research philosophy of the presented study is that we didn't use a top-down approach where we "assign" a simple 1D generative model and "design" a circuit model to implement the inference. Instead, we ask the question the other way around from a bottom-up view. That is, given the recurrent circuit models with types of interneurons that have been extensively studied in neuroscience in recent years (acknowledged by the Reviewer RifV by saying this is "the norm of papers in neural sampling theory"), what kind of generative model and probabilistic inference can be implemented by the biological circuit. This philosophy was detailed in Lines 139-143 and Lines 152-156 of the submitted manuscript.

In this way, even if the discovered 1D sampling in the proposed circuit model is simple, we feel the result still has its merit to be reported, because

1. no previous study has proposed what kind of inference can be achieved in the circuit with types of interneurons;
2. the simple 1D inference will encourage people to further investigate how the canonical circuit could do more complex inference.

For the results uploaded in the Global Rebuttal PDF, we found the circuit with SOM neurons could sample multi-modal stimulus. However, due to the page limit, we didn't include them in the submitted manuscript. And we are happy to include them in the revised manuscript if you would prefer.

We appreciate the reviewer's positive feedback on our writing and math analysis!

My main reservation with the paper is on its significance - I think the paper (at least in the way it is written) interests people within the specific subfield of Bayesian sampling in simplified circuit models for its technical merits; the result of SOM leads to faster inference for simple 1D input may interest SOM researchers in neuroscience. I'm just debating how much of these merit a NIPS publication.

Understanding how neural circuits in the brain do computation is one of NeurIPS's interests, which is also acknowledged by the Reviewer K6Di (strengths). We believe this adventure would bring mutual benefit to both neuroscience and machine learning (ML) communities, in that it not only helps us understand the brain, but also has the potential to be a new building block for ML tasks (Of course, our current model is still a "baby" and is not able to deal with real tasks). There are lots of famous examples of brain-inspired computations, e.g., Neocognition by Fukushima, Conv Net by Yan LeCun, etc.

Even if the manuscript only presents a 1D posterior sampling, its analytical methodology and research philosophy (analytically identified circuit representation and algorithm) further our understanding of the uncertainty representation in neural networks, and how different types of interneurons contribute to sampling, which, we believe, will have long-term impacts on ML. In contrast, conventional perceptron-based neural network models ignored neural fluctuations when they were developed decades ago, which causes uncertainty representation in artificial neural networks, a fundamental question in ML. Moreover, the diversity of interneuron types may provide new possibilities to implement/link with the complicated components used in conventional RNNs used in ML such as LSTM. At last, the 1D sampling is the preface of our series of work, and we have also considered extending the circuit model to deal with more complicated posteriors (please see our Global Rebuttal and Rebuttal Figure).

Have the authors considered having a mixture of two orientations as input and study how the introduction of I-diversity resolves the potential interference? (mentioning it not to ask for more experiments for this submission but just as a general future direction)

Yes, the proposed circuit with SOM neurons can sample multi-modal and high-dim posteriors. Please find the results in the rebuttal PDF and our response in Global Rebuttal. We are happy to include the results to strengthen our paper in the revised version with one more page, if you encourage us to do so. The reason for not including them is because of the page limit and the completeness of the manuscript, where we sacrificed the inference task complexity to present the reasoning of building the model and math analysis.