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

Authors: eryn  sale
wen-hao  zhang
Venue: NeurIPS_cc_2024_Conference
Published: 2024-05-14 20:47:54

Abstract: Accumulating evidence suggests stochastic cortical circuits can perform sampling-based Bayesian inference to compute the latent stimulus posterior. Canonical cortical circuits consist of excitatory (E) neurons and types of inhibitory (I) interneurons. Nevertheless, nearly no sampling neural circuit models consider the diversity of interneurons, and thus how interneurons contribute to sampling remains poorly understood. To provide theoretical insight, we build a nonlinear canonical circuit model consisting of recurrently connected E neurons and two types of I neurons including Parvalbumin (PV) and Somatostatin (SOM) neurons. The E neurons are modeled as a canonical ring (attractor) model, receiving global inhibition from PV neurons, and locally tuning-dependent inhibition from SOM neurons.
We theoretically analyze the nonlinear circuit dynamics and analytically identify the Bayesian sampling algorithm performed by the circuit dynamics. We found a reduced circuit with only E and PV neurons performs Langevin sampling, and the inclusion of SOM neurons with tuning-dependent inhibition speeds up the sampling via upgrading the Langevin into Hamiltonian sampling. Moreover, the Hamiltonian framework requires SOM neurons to receive no direct feedforward connections, consistent with neuroanatomy. Our work provides overarching connections between nonlinear circuits with various types of interneurons and sampling algorithms, deepening our understanding of circuit implementation of Bayesian inference.

Main paper: 
 introduction: 
 The brain lives in a world of uncertainty and ambiguity and thus has to infer unobserved world states. The Bayesian inference is a normative framework to implement inference, and extensive studies have suggested that the brain's perception is consistent with the Bayesian inference, forming the concept of Bayesian brain [1] , including, e.g., visual processing [2] , multisensory integration [3] , decision-making [4] , sensorimotor learning [5] , etc. Studying how neural circuits in the brain realize Bayesian inference has been an active topic in neuroscience [6] [7] [8] . Many neural circuit models of Bayesian inference have been developed with distinct representational and algorithmic mechanisms, e.g., parametric-based representation [4, [8] [9] [10] [11] and sampling-based representation [12] [13] [14] [15] [16] [17] [18] [19] .Despite a large body of neural circuit models of Bayesian inference, there are still gaps between our current Bayesian circuit models and canonical recurrent circuits in the cortex. One obvious distinction is previous Bayesian circuit models haven't considered the rich diversity of neuronal types in the cortex, especially inhibitory interneurons. The canonical cortical microcircuit contains three major types of inhibitory (I) interneurons [20] [21] [22] [23] [24] , including Parvalbumin (PV), Somatostatin (SOM), and Vasoactive Intestinal Peptide (VIP) neurons (Fig. 1A ). These interneurons have different electrical properties, stimulus-tuning profiles, distinct connectivity, and synaptic modulations with other neurons. For example, PV neurons are weakly tuned to stimulus [20] , and multiplicatively modulate E neurons in a way called divisive normalization by sending axons to E neurons' cell body [21, 22, 25, 26] . In contrast, SOM neurons are tuned to the stimulus, and module E neurons in an additive way via sending axons to the distal dendrites of E neurons [22, 25, 26] . Another gap comes from the nonlinearity of cortical circuit dynamics, which impedes the analytical understanding of circuits' Bayesian inference algorithms. Most earlier studies relied on numerical methods to analyze the algorithm in nonlinear circuits (e.g., [15] [16] [17] 19] ), or considered linear neural dynamics to obtain analytical solutions (e.g., [12, 18, 27, 28] ). We still lack a comprehensive understanding about how the nonlinear recurrent circuit dynamics with diversity of interneurons implement Bayesian inference.To gain insight into Bayesian computation in a nonlinear recurrent circuit with types of interneurons, we build a canonical recurrent circuit model consisting of excitatory (E) and two types of interneurons, including PV and SOM neurons (Fig. 1B ), and investigate how the model implements sampling-based Bayesian inference. The E neurons are modeled as a rate-based ring (attractor) model that emerges the tunings over a 1D stimulus such as orientation. The E neurons receive internal Poisson-like variability mimicking stochastic spike generation, which provides a variability source to drive sampling [18] . For simplicity, the PV neurons in the model are not tuned to the stimulus as a limiting case of their weak tuning found in experiments [20] , and provides global inhibition to E neurons via divisive normalization to ensure stability [21, 25] . In contrast, the SOM neurons have stimulus tunings and provide locally tuned inhibitory feedback to E neurons in an additive way [20, 29] . The circuit model with the above connectivity successfully reproduces multiplicative and additive modulations on E neurons' tunings from PV and SOM neurons respectively (Fig. 1G-H ) [26] .We perform theoretical analysis on the nonlinear recurrent circuit dynamics, and analytically identify the sampling algorithm adopted in the circuit. We find the reduced circuit with only E and PV neurons can implement Langevin sampling in the stimulus feature manifold. The tuning-dependent inhibitory feedback from SOM speeds up the sampling by upgrading the Langevin sampling into Hamiltonian sampling. And the two types of interneurons have different effects on sampling speed. Moreover, we find that Hamiltonian sampling requires SOM neurons not to receive feedforward sensory inputs, consistent with neuroanatomy with few feedforward synapses targeting SOM neurons [22, 24] . The nonlinear circuit model with fixed weights can flexibly sample posteriors with different uncertainties, if located in the linear input-output regime. At last, the circuits can be extended to sample multivariate stimulus posteriors and bimodal posteriors. 
the recurrent neural circuit with various types of interneurons: 
 The cerebral cortex is a repetition of the canonical neural circuit composed of multiple types of neurons (Fig. 1A ), including excitatory (E) neurons and three major types of inhibitory (I) interneurons (PV, SOM, and VIP; classified via biomarkers [30] ). To study how sampling-based Bayesian inference is implemented by the canonical neural circuit composed of various types of neurons, we build a recurrent neural circuit model consisting of E neurons and two types of I interneurons of PV and SOM neurons (Fig. 1B ). The model doesn't include VIP neurons, which will form our future research (see Discussion). The basic wiring diagram of the proposed circuit model is consistent with the structure of the canonical cortical circuit (Fig. 1A ). In the model, the E neurons are selective for a 1D periodic stimulus feature z ∈ (-π, π], e.g., the orientation moving direction. Denote θ j as the preferred stimulus feature of the j-th E neuron, and the preferred stimulus features of all N E E neurons, {θ j } N E j=1 uniformly cover the whole range of feature space z (Fig. 1C and F ). This setting is the same as the canonical ring network model that has been widely used in modeling cortical circuits (e.g., [17, [31] [32] [33] [34] [35] [36] ). Mathematically, in the continuum limit (θ j → θ) corresponding to an infinite number of neurons, the dynamics of the E neurons is [34, 37] ,EQUATIONwhere u E (θ, t) and r E (θ, t) are the synaptic input and firing rate respectively of the E neuron preferring z = θ. X denotes neuronal types with E, F and S representing E neurons, sensory feedforward inputs, and the SOM neurons respectively. τ is the time constant of synaptic input, The model result is qualitatively similar to the experimental data from [26] . and ρ = N E /2π is the neuronal density covering the space of stimulus feature z. E neurons receive an internal Poisson variability mimicking stochastic spike generation (Eq. 1, last term), and F E the Fano factor of the injected variability whose value will be adjusted to make the Fano factor of E neurons' activities at the order of 1. ξ(θ, t) is standard Gaussian white noise, i.e., ⟨ξ(θ, t)ξ(θ ′ , t ′ )⟩ = δ(θ -θ ′ )δ(t -t ′ ) with δ(t -t ′ ) being the Dirac delta function. The internal Poisson variability provides the variability for the circuit to draw random samples from the posterior [18] , which may arise from the E and I balance in the cortex as a complex network effect [38] [39] [40] .Recurrent connection kernels. W Y X (θ) denotes the recurrent connection kernel from neurons with type X to those with type Y , which are modeled as Gaussian functions in the model (Fig. 1D ),EQUATIONwhere w Y X controls the peak strength of the recurrent weight, and a the connection width across the stimulus feature space. The kernels W Y X (θ) connecting different neuronal types can have different peak weights w Y X (w Y X > 0 or < 0 regards to E or I synapses respectively). Moreover, different W Y X (θ) may have different connection widths a XY , although most of them have the same width a XY = a unless noted otherwise (see Supplementary Information (SI.) Sec. 6.1). In Eq. ( 1), the symbol * denotes the spatial convolution, i.e., (W * r)(θ) = W(θ -θ ′ )r(θ ′ )dθ ′ , which implies the translation-invariance of the connection weight between neurons in the stimulus feature space.Sensory feedforward inputs. The recurrent circuit model receives sensory feedforward input r F (θ, t) (Eq. 1) randomly evoked from a stimulus feature z in the world. Given a stimulus feature z, we assume the feedforward inputs r F are conditionally independent Poisson spikes with Gaussian tuning (Fig. 2A-B ), which has been widely used before (e.g., probabilistic population code [4, 8, 9] ).EQUATIONwhere λ F (θ|z) is the mean firing rate. In simulating our rate-based model, the r F is approximated as a continuous Gaussian random variable with multiplicative noise to mimic the Poisson statistics.PV neurons. The stimulus orientation weakly modulates the PV neurons [20, 26, 41] , hence, for simplicity, we consider PV neurons in the model are not tuned for stimulus features and only provide global unstructured inhibition to E neurons to keep stability. Moreover, it was suggested PV neurons provide divisive normalization (DN) to modulate E neurons' responses via shunting inhibition [21, 25, 42] . Hence the proposed circuit model absorbs PV neurons' effects in the divisive normalization of E neurons which has been widely used in circuit models [21, 25, 43, 44] ,EQUATIONwhere the DN acts as an activation function transferring the instantaneous synaptic input u E (θ, t) of E neurons into their firing rate r E (θ, t).[x] + = max(x, 0) denotes negative rectification. The integral [u E (θ, t)] 2 + dθ ′ ≡ r P V , reflects PV neurons globally summing all E neurons' activities. w EP is the global inhibition strength characterizing the inhibitory weight from PV to E neurons. SOM neurons. It is suggested that SOM neurons linearly modulate E neurons' responses, in contrast to multiplicative modulation from PV to E neurons [26] . Therefore the model considers the E neurons receive additive synaptic inputs from SOM neurons ((W ES * r S )(θ, t), Eq. 1). Furthermore, SOM neurons are tuned to a stimulus feature with a strength comparable to E neurons [20] , unlike the weak tunings of PV neurons. Thus, the dynamics of SOM neurons are governed by,EQUATIONwhere the g S (scalar) controls the "gain" of SOM neurons and is set as a fixed value across the study. Two features about SOM neurons' connectivity are worth noting (Eq. 5). First, the model doesn't include recurrent inhibitory connections between SOM neurons, due to few mutual inhibitions between SOM neurons [22, 24] . This simplification will only change the effective gain g s of SOM neurons without affecting our conclusions of the circuit algorithm. Second, the proposed circuit model allows the existence of the feedforward connections to SOM neurons (W SF in Eq. 5), even if they are rare in reality [22, 24] . The reason for allowing this rare connection is we want to test whether the Bayesian sampling theory has the power to constrain it.Overall, the proposed circuit model is consistent with most canonical recurrent circuit models in the field (e.g., [31] [32] [33] [34] [35] [36] ). The simplifications considered above reserve the main characteristics of neuronal connectivity and response properties observed in experiments, especially interactions between E neurons and interneurons. For example, it was found enhancing PV neurons will multiplicatively modulate E neurons' tuning curves (Fig. 1G , right), whereas SOM neurons modulate E neurons' tuning additively (Fig. 1H , right) [26] . Both effects are successfully reproduced in the proposed circuit model (Fig. 1G-H , left; see details in SI. Sec. 6.5). 
from recurrent circuit dynamics to bayesian inference: 
 The proposed recurrent circuit dynamics (Eqs. 1 and 5) is supposed to implement Bayesian inference by computing the stimulus posterior based on a received feedforward input, p(z|r F ) ∝ p(r F |z)p(z).(6) It regards the stage from external stimulus z to the feedforward input r F as the probabilistic generative process (Fig. 2A ). Implementing Bayesian inference requires the recurrent circuit to store the generative model. Rather than defining a generative model and deriving its neural circuit implementation as in many previous studies, here we ask the question the other way around: if the proposed recurrent circuit model based on neurophysiology could do Bayesian inference (Eq. 1A), what generative model is stored in the circuit and what stimulus posteriors are computed by the circuit? Furthermore, what is the Bayesian inference algorithm adopted by the circuit dynamics?Stimulus likelihood. The stochastic feedforward input from the stimulus feature z (Eq. 3) naturally specifies the stimulus likelihood that can be calculated as a Gaussian likelihood (see SI. Sec. 2.1),EQUATIONwhere the Gaussian stimulus likelihood function comes from the the Gaussian profile of feedforward input tuning λ F (θ|z) (Eq. 3) [9] . The mean µ z and the precision Λ of the stimulus likelihood can be read out from r F via a linear decoder called population vector [45, 46] ,EQUATIONGeometrically, µ z is regarded as the location of r F in the stimulus feature space, and Λ is proportional to the input spike count (Fig. 2B-C ). In this way, a single snapshot of r F parametrically conveys the whole stimulus likelihood function p(r F |z) [9] .Subjective prior. We suppose the recurrent circuit utilizes its stored subjective prior to compute the subjective stimulus posteriors, with a mild assumption that the subjective prior in the circuit matches the objective prior in the outside world. Nevertheless, the subjective circuit prior remains unknown at this point. Next, we theoretically analyze the circuit dynamics to identify the stored subjective prior in the circuit and find out the circuit algorithm of Bayesian sampling. 4 The Bayesian sampling in the stochastic circuit dynamics 
theoretical analysis of the neural dynamics: 
 We theoretically analyze the nonlinear circuit dynamics to investigate how it can implement Bayesian sampling. We perform perturbative analysis of the nonlinear circuit dynamics, identify the lowdimensional stimulus feature manifold (subspace) in the high-dimensional population response space, and eventually study the circuit dynamics on the stimulus feature manifold. Given a feedforward input r F (Eq. 3), it can be checked that the synaptic input u X (θ) and the firing rate r X (θ) of neurons X in the equilibrium attractor states are both Gaussian profiles (Fig. 1E-F ; see SI. Sec. 3.1) , ⟨u X (θ)⟩ = U X exp[-(θ -zX ) 2 /4a 2 ], ⟨r X (θ)⟩ = R X exp[-(θ -zX ) 2 /2a 2 ], (X = E, S) (9)where ⟨•⟩ denotes the average over different realizations. U X and R X denote the height of the population synaptic input and firing rate respectively, and can be analytically computed (Eq. S37). The position of population activity on the stimulus feature space is zX = µ z is the same as the location µ z of the feedforward inputs (SI. Sec. 3.1). Intuitively, this is because the recurrent circuit is homogeneous along the stimulus feature space, i.e., all neurons are uniformly distributed on the stimulus feature space and the recurrent connections are translational invariance.With the corruption of sensory noises and the internal Poisson variability, the instantaneous neural responses will deviate from the equilibrium attractor state (Eq. 9). We treat each instantaneous response as a perturbation from its equilibrium state, i.e., u X (θ, t) = ⟨u X (θ)⟩ + δu X (θ, t), (X = E, S), and the relaxation dynamics of the perturbation δu X (θ, t) can be derived [47] . Then performing eigen-analysis of the perturbation dynamics we analytically find out the stimulus feature manifold (subspace), which is specified by its (unnormalized) eigenvector [35, 47, 48] ,EQUATION)Previous studies have shown the stimulus feature eigenvector has the largest eigenvalue in the perturbation dynamics [47] . Therefore, we project the dynamics of E and SOM neurons (Eqs. 1 and 5) onto their respective stimulus feature eigenvectors, where the projection is computing the inner product between the neuronal responses and the eigenvector, i.e., ⟨ϕ(θ), f (θ)⟩ = ϕ(θ)f (θ)dθ, with f (θ) representing the left or right handed side of Eqs. (1 and 5). The projection yields the dynamics of the E and SOM neurons on the stimulus feature manifold (see details in SI. Sec. 3.3)EQUATIONEQUATIONThe approximation comes from ignorance of some negligible nonlinear terms. z E and z S denote the instantaneous positions of neural activities at the stimulus feature manifold at time t (Fig. 2E ). µ z is the observed stimulus feature conveyed by the feedforward input (Eq. 8). τX = τ U X (X = E, S)is the time constant of the circuit dynamics on the stimulus feature manifold, where U X is the peak value of population synaptic input (Eq. 9). The coefficients g XY in the above equation denote the coupling strength between neural response positions, where g XY ∝ w XY R Y with R Y the peak firing rate of the pre-synaptic neural population (Eq. S48). Theσ 2 E = 8aF E /(3 √ 3π)is the variance of the internal variability on the stimulus feature manifold coming from the internal Poisson variability (Eq. 1, last term). σ E is a constant value irrelevant with feedforward inputs and network responses.Since the circuit contains multiple types of neurons, we first analyze a reduced circuit dynamics on the stimulus feature manifold without SOM neuron (Fig. 2A , setting g ES to zero in Eq. 11) and then study how the SOM neurons affect the sampling dynamics. Without SOM neurons (only E and PV neurons), the E neurons' dynamics on the stimulus feature manifold is,EQUATIONwhich is a first-order Langevin dynamics. This motivates the possibility that the E and PV neurons can implement Langevin sampling in the stimulus feature manifold. If true, the instantaneous position of E activity, z E (Fig. 2D ), can be regarded as a stimulus feature sample generated by the circuit.In theory, the Langevin sampling of a posterior p(z|r F ) ∝ p(r F |z)p(z) corresponds to performing stochastic ascent on the log posterior surface [49, 50] ,EQUATIONwhere τ z is the sampling time constant controlling the sampling speed. The 2nd row is obtained by substituting the Gaussian likelihood (Eq. 7).Uniform subjective prior. Comparing Eqs. ( 14) and ( 13), we can identify a stored uniform stimulus prior p(z) in the reduced circuit model stores . This is because the gradient of a uniform prior is ∇ ln p(z) = 0, and then the Langevin sampling dynamics (Eq. 14) reduces to a form similar to the circuit dynamics on the stimulus manifold (Eq. 13). The uniform stimulus prior comes from homogeneous neurons in the circuit, i.e., neurons are uniformly distributed along the stimulus feature space, and the translation-invariant connection profile (Eq. 2). It implies that the circuit needs to break the symmetry of homogeneous neurons to store a non-uniform prior (see Discussion).Condition for realizing Langevin sampling. Utilizing Langevin dynamics to sample the posterior requires the drift and diffusion coefficients to share the same time constant τ z (Eq. 14). To satisfy this requirement, the E dynamics on the stimulus feature manifold (Eq. 13) should have g EF /σ 2 E = Λ/2. With the expressions of g EF and σ 2 E (Eq. S48), the feedforward connection weight w EF should be,EQUATIONIntuitively, larger internal variability F E increases the sampling variance, and it requires a larger feedforward weigh w EF to compensate for sampling variance increase to match with the posterior variance. To verify our theoretical derivation (Eq. 15), we search whether there is an optimal value of the feedforward weight w EF allowing the circuit without SOM neurons to sample the posterior (likelihood). Indeed, we find once the recurrent circuit model is set with that optimal w EF , the reduced circuit model (consist of E and PV neurons) with all parameters fixed can flexibly sample the likelihood with various uncertainties (Fig. 2G ). A characteristic of Langevin sampling is the crosscorrelation function of samples (Eq. 14) exponentially decays with time which can be calculated as ρ(∆t) = exp(-g EF ∆t/τ E ) (SI. Eqs. S7). To verify whether the reduced circuit performs sampling with the Langevin dynamics as suggested by Eq. ( 13), we estimate the cross-correlation function of stimulus sample z E generated by the network, which indeed exhibits an exponential form, and our theoretical calculation ρ(∆t) predicts the actual cross-correlation function well (Fig. 2F ).The SOM neurons augment the dimensionality of the circuit dynamics on the stimulus feature manifold from the first order to the second order (Eqs. 11 and 12). In principle, the second-order Hamiltonian sampling dynamics motivates us to explore whether and how the circuit can implement Hamiltonian-like sampling [50, 51] . The computational benefit from the increased complexity of Hamiltonian sampling compared with Langevin sampling is accelerating the sampling speed by generating less correlated samples over time [50] . To investigate such a possibility, we write one commonly used Hamiltonian sampling dynamics [50, 52, 53] ,EQUATIONy is the auxiliary variable representing momentum in the Hamiltonian sampling. It can speed up sampling by increasing sampling step size when close to the posterior center, helping the sampling trajectories move toward the other side of the posterior (Fig. 3B ; details in SI. Sec. 2.1). β is the friction strength determining how fast the momentum will decay to zero. We next bridge the circuit dynamics (Eqs. 11 and 12) with the Hamiltonian sampling dynamics (Eq. 16).Intuitively, the sample from SOM neurons z S (Eq. 12) resemble the auxiliary variable y in the Hamiltonian dynamics (Eq. 16). Nevertheless, there are several gaps between the two: First, in the Hamiltonian sampling, the stimulus sample z is purely driven by the auxiliary variable y (Eq. 16), whereas in the circuit model z E receives both z S and µ z (Eq. 11). Second, Hamiltonian sampling injects variability into the auxiliary dynamics y (Eq. 16), while the variability (from stochastic spike generation) is injected into the dynamics of z E . To bridge the gap between the circuit model and sampling dynamics, we assume that the circuit is conducting a mixture of Langevin (Eq. 14) and Hamiltonian sampling (Eq. 16), and thus we split the E neurons' sampling dynamics into two parts,EQUATIONwhere α L ∈ [0, 1] denotes the proportion of feedforward input contributed by Langevin sampling. We see the y S resembles the auxiliary variable y in Hamiltonian sampling (Eq. 16), which implies the auxiliary variable in the circuit is a mixture of samples from E neurons z E and SOM neurons z S . Given the definition of y S , the stochastic dynamics of y S can be derived as a form similar to the one in the Hamiltonian sampling (Eq. 16; details in SI. Sec. 4.1).EQUATIONβ y , β E , β S and σ y are functions of the coefficients in Eq. ( 17) (detailed expressions at Eq. S53).It can be checked that equilibrium distribution of the mixed dynamics is the posterior (details of Fokker-Planck approach in SI. Sec. 2.5).Conditions for realizing mixed sampling. To utilize the full circuit model with SOM neurons (Eqs. 17-18) to implement Bayesian sampling, the coefficients in the circuit model should satisfy the relations of coefficients required in the Langevin sampling (Eq. 14) and the Hamiltonian sampling (Eq. 16). First, setting up the Langevin sampling part in the circuit model (Eq. 17, blue) requires α L g EF /σ 2 E = Λ/2, which finally leads to,EQUATIONwhose value is 1/α L times the feedforward weight in the Bayesian sampling circuit without SOM neurons (Eq. 15). Second, realizing the Hamiltonian sampling part in the circuit model (Eq. 17, orange; and Eq. 18) yields three conditions shown below,EQUATION)The conditions in Eqs. ( 19) and ( 20) combined enable the full circuit model with SOM neurons to implement Bayesian sampling. An important insight from Eq. ( 20b) is that the SOM neurons should not receive feedforward inputs directly (w SF = 0 in Eq. 5; Fig. 1A , dashed line removed). This theoretical result is consistent with the anatomy that SOM neurons receive much fewer feedforward synapses than other types of neurons [23, 24] . In addition, substituting the detailed expression of coefficients into Eq. ( 20), we find there is a low-dimensional manifold in the circuit model's connection weight space for the circuit to sample the posterior p(z|r F ) (SI. Sec. 4.2),EQUATIONG(α L ) is a nonlinear function of α L which specifies the proportion of Langevin sampling in the circuit dynamics, which remains invariant with feedforward input and network activities (SI. Eq. S65). U X and R X are the height of the population synaptic input and firing rate of neuronal populations X (Eq. 9). Eq. ( 21) implies that the sampling in the full circuit model is robust, without the need to fine-tune recurrent weights. To verify the theoretical result (Eq. 21), we fix the weight from E to SOM neurons, w SE , and search whether there is a line manifold in the two-dimensional parameter space of w ES and w EF for the circuit correctly sample the stimulus posterior. Indeed, Fig. 3E numerically confirms the line manifold of weights under which the sampling distribution matches the posterior. Moreover, the introduction of SOM interneruons makes the cross-correlation of sample z E decay faster (Fig. 3C ), suggesting speeding up sampling (Fig. 3D ).Flexible sampling posteriors in the linear regime. Moreover, the circuit model with fixed weights should be able to sample the posteriors with different uncertainties. In the circuit model, the posterior uncertainty is determined by the feedforward input rate R F (Eq. 3), where a larger input rate leading to smaller posterior uncertainty (Eq. 8). We find that when the nonlinear circuit model is located at the linear regime, it can flexibly sample posteriors with different uncertainties. To see this effect, we can change the feedforward input rate by multiplying a gain factor g, i.e., R F → gR F that change the likelihood precision Λ → gΛ (Eq. 8 and S5). If the circuit is at the linear regime, the peak value of synaptic input, U X → gU X , and the population firing rate, R X → gR X , are both multiplied with the gain factor g being applied to the feedforward input. And then it can be checked Eq. ( 21) is still satisfied. This theoretical result is confirmed by numerical simulation (Fig. 3F ) which shows the precision of circuit's stimulus samples increases with feedforward rate R F and aligns well with the likelihood. Here we adjust the recurrent E-to-E weight w EE to set the circuit has an approximately linear response at the range of the feedforward input rate. Fig. S1 shows if the network deviates from the linear regime, the circuit's sampling distribution will deviate from the likelihood.We further investigate how quantitative measures of sampling, e.g., sampling speed and temporal correlation of samples, will be affected by interneurons such as the inhibitory feedback weight. In principle, both sampling speed and temporal correlation can be revealed by the eigenvalues of the circuit sampling dynamics (Eqs. [11] [12] . When organizing the circuit sampling dynamics into a matrix form (Eqs. 11-12), ż = -Mz + µ, with z = (z E , z S ) ⊤ and µ lumping terms exclusive z E or z S in Eqs. (11 -12) , the eigenvalues of the circuit sampling dynamics are (SI. Eq. S69),λ ± = tr(M) ± tr(M) 2 -4 det(M) ≜ tr(M) ± √ ∆,where tr(M) = τ -1 The sampling speed is limited by the smallest real part of eigenvalues, i.e., Re(λ -) (Eq. 22) , in that the KL divergence from the distribution of samples up to time t to the stationary distribution (Fig. 3D ) decays exponentially as exp[-Re(λ -)t] [54] .EQUATIONPV neurons. PV's inhibition weight w EP modulates the eigenvalues by decreasing the common factor τ E = τ U E (w EP ) (Eq. 22) , where U E (w EP ), the peak value of E population synaptic input (Eq. 9), decreases with w EP . Hence, stronger PV inhibition increases the slowest eigenvalue λ -(Fig. 4A ) and leads to faster sampling, exhibited by the faster decay of the temporal correlation between samples (Fig. 4A , blue and red). Moreover, the multiplicative modulation from PV neurons will not induce the imaginary part of eigenvalues (Fig. 4B , bottom), i.e., no oscillation between samples. SOM neurons. The SOM inhibitory weight w ES will have non-monotonic effects on sampling speed measured by the real part of the slowest eigenvalue. There is a value of w ES to maximize the sampling speed (Re(λ -) (Fig. 4C , top). Moreover, SOM's inhibition will induce temporal oscillations between samples, i.e., emerging imaginary part of eigenvalue λ -(Fig. 4C, bottom ). This oscillation is confirmed by the cross-correlation of samples (Fig. 4A , purple). Moreover, to mimic neural experimental data analysis, the oscillation induced by SOM neurons can be revealed by the power spectrum analysis of the local field potential (LFP) (SI. Eq. S90). 
sampling complex posteriors in canonical recurrent circuits: 
 High-dimensional stimulus posteriors. As a proof of concept example, we consider sampling bivariate stimulus posteriors by coupled circuits (Fig. 5A ) with each circuit the same as Fig. 1B . And only E neurons across circuits are coupled. Each circuit m (m = 1, 2) receives a feedforward input generated by a latent stimulus feature z m , and will sample z m . Hence, the number of coupled circuits equals to the stimulus feature dimension. We found the coupled circuits store an associative (subjective) prior, i.e., p(z 1 , z 2 ) ∝ exp[-Λ s (z 1 -z 2 ) 2 /2] (Fig. 5C ), with the prior precision Λ s increasing with inter-circuit coupling weights (Fig. 5D ; Eq. S85). Math analysis is presented in SI. Sec. 5. Concatenating the samples generated by E neurons in two circuits (with the same readout as Fig. 3A-B ), we can obtain the 2D posterior sampled from the coupled circuit (Fig. 5B ). Similarly, the SOM interneurons speed up sampling, reflected by the sampling trajectories with SOM transverse over a wider posterior region than the one without SOM in the same period (Fig. 5B ).Bimodal stimulus distributions. We also use the circuit (Fig. 1C ) to sample uni-variate bimodal posteriors (Fig. S2 ), in response to superpositions of two feedforward inputs with each generated by a latent stimulus. The circuit can sample a bimodal distribution, where samples jump between two modes alternatively, due to the bi-stability in the circuit dynamics (Fig. S2C and E ). In contrast, without SOM neurons (Fig. 2A ), the circuit can only sample a uni-modal distribution, as a uni-modal approximation of the bimodal one (Fig. S2F ). Due to the space limit, we haven't comprehensively linked the circuit' bimodal sampling distribution with posteriors, which will form our future work.Question: Does the paper provide open access to the data and code, with sufficient instructions to faithfully reproduce the main experimental results, as described in supplemental material?Answer: [Yes] Justification: Codes are uploaded and shared for reproduction and further experiments.Guidelines:• The answer NA means that paper does not include experiments requiring code. • Providing as much information as possible in supplemental material (appended to the paper) is recommended, but including URLs to data and code is permitted. 
conclusion and discussion: 
 The present study investigates how canonical recurrent circuits with diverse inhibitory interneurons implement Bayesian sampling. The nonlinear circuit model consists of E neurons, and two types of interneurons neurons (PV and SOM neurons). PV and SOM neurons have distinct tuning properties and modulations on E neurons: PV neurons have weak stimulus tuning and multiplicatively modulate E neurons, while SOM neurons have stimulus tuning and send additive inputs to E neurons. Through theoretical analysis and numerical simulations, we find the reduced circuit with only E and PV neurons implements Langevin sampling in the stimulus manifold (subspace) to compute the stimulus posterior.The SOM neurons accelerate the circuit sampling speed by upgrading the Langevin sampling into Hamiltonian sampling. We further find Hamiltonian sampling requires the SOM neurons to receive no feedforward connections, consistent with neuroanatomy. We also investigate how inhibitory strength from two types of interneurons affects the sampling speed. Our work is one of the earliest studies investigating the Bayesian sampling algorithm in canonical nonlinear recurrent circuits with diverse interneurons, and provides new insight into interneurons in Bayesian computation.Comparison with other work. The computational mechanism of accelerated sampling by SOM neurons is similar to previous studies considering structured inhibitory feedback to E neurons (e.g., [14, 17, 53, 55] ), while some notable differences exist. Earlier studies considered the sampling in neural response space where the posterior dimension is the same as the number of neurons [14, 17] , whereas the current circuit samples in the stimulus feature manifold (subspace) in the neural response space. On the other hand, although other previous studies considered sampling in the stimulus feature manifold similar to the current model [53, 55] , their structured inhibitory feedback comes from the biophysical mechanism within single neurons, e.g., spike frequency adaptation [53] , and potassium channels [55] . From a neurobiology perspective, SOM neurons can be modulated by VIP neurons (Fig. 1A ), whereas intracellular mechanisms are hardly modulated [53, 55] , suggesting accelerated sampling from SOM neurons might be more flexible than single neuron mechanism in reality. Moreover, the proposed circuit model is similar to a recent one [18] in terms of utilizing internal Possion variability to draw samples. However, the current circuit model is nonlinear while the other study considered a linear circuit model [18] . Lastly, the coupled circuit sampling of bivariate posteriors was also considered in [34] , but which didn't figure out the circuit's sampling algorithms.Limitations and extensions of the model. The proposed circuit model doesn't include VIP neurons that exclusively target SOM neurons (Fig. 1A ). Our future work will incorporate them and study their effects on circuit sampling. Based on our current conclusion, it is likely that VIP neurons act as a "knob" to modulate the sampling speed, depending on task needs, by changing the activation level ("gain") of SOM neurons. Moreover, we find the proposed canonical circuit stores a uniform stimulus prior, due to the homogeneity of neurons distributing on the stimulus manifold. The homogeneous neuron simplification has been widely used in neural coding and continuous attractor networks (e.g., [31, 34, 35, 56] ), which simplifies the math analysis without altering results substantially. Nevertheless, the circuit has to break the neuronal homogeneity to store a non-uniform prior [57] , and certainly, cortical neurons are heterogeneous. The neuronal heterogeneity can be realized by manipulating the translation-invariant recurrent connection matrix (Eq. 2), e.g., introducing randomness (zero mean with certain variance) on recurrent weight which has also been widely used in (chaotic) Excitation and Inhibition (E/I) balanced networks ( [38] [39] [40] 58] ). 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, which is statistically equivalent to the injected multiplicative variability in our rate-based network (Eq. 1, last term). Moreover, the proposed circuit model only infers a static stimulus, and we will extend to infer a dynamic stimulus described by a hidden Markov model in the future. All of these form our future research.Question: Does the paper specify all the training and test details (e.g., data splits, hyperparameters, how they were chosen, type of optimizer, etc.) necessary to understand the results?Answer: [NA] Justification: The present paper doesn't study the learning problem.Guidelines:• The answer NA means that the paper does not include experiments.• The experimental setting should be presented in the core of the paper to a level of detail that is necessary to appreciate the results and make sense of them. • The full details can be provided either with the code, in appendix, or as supplemental material. 
experiment statistical significance: 
 Question: Does the paper report error bars suitably and correctly defined or other appropriate information about the statistical significance of the experiments?Answer: [Yes] Justification: We define the error bars in Fig. 1G-H , and the define statistics of simulating the proposed circuit model in text corresponding to Fig. 2 3 4 , including the mean, correlation, empirical distribution, and the KL divergence. The Supplementary Information has details of how we obtain these statistics.Guidelines:• The answer NA means that the paper does not include experiments.• The authors should answer "Yes" if the results are accompanied by error bars, confidence intervals, or statistical significance tests, at least for the experiments that support the main claims of the paper.• The factors of variability that the error bars are capturing should be clearly stated (for example, train/test split, initialization, random drawing of some parameter, or overall run with given experimental conditions). • The method for calculating the error bars should be explained (closed form formula, call to a library function, bootstrap, etc.) • The assumptions made should be given (e.g., Normally distributed errors). • It should be clear whether the error bar is the standard deviation or the standard error of the mean. Justification: It can be found at the SI. Sec. ??.Guidelines:• The answer NA means that the paper does not include experiments.• The paper should indicate the type of compute workers CPU or GPU, internal cluster, or cloud provider, including relevant memory and storage. • The paper should provide the amount of compute required for each of the individual experimental runs as well as estimate the total compute. • The paper should disclose whether the full research project required more compute than the experiments reported in the paper (e.g., preliminary or failed experiments that didn't make it into the paper). 
code of ethics: 
 Question: Does the research conducted in the paper conform, in every respect, with the NeurIPS Code of Ethics https://neurips.cc/public/EthicsGuidelines? 
broader impacts: 
 Question: Does the paper discuss both potential positive societal impacts and negative societal impacts of the work performed?Answer: [No] Justification: This is a theoretical study of basic neuroscience research and will not have direct societal impacts.Guidelines:• The answer NA means that there is no societal impact of the work performed.• If the authors answer NA or No, they should explain why their work has no societal impact or why the paper does not address societal impact. • Examples of negative societal impacts include potential malicious or unintended uses (e.g., disinformation, generating fake profiles, surveillance), fairness considerations (e.g., deployment of technologies that could make decisions that unfairly impact specific groups), privacy considerations, and security considerations. • The conference expects that many papers will be foundational research and not tied to particular applications, let alone deployments. However, if there is a direct path to any negative applications, the authors should point it out. For example, it is legitimate to point out that an improvement in the quality of generative models could be used to generate deepfakes for disinformation. On the other hand, it is not needed to point out that a generic algorithm for optimizing neural networks could enable people to train models that generate Deepfakes faster. • The authors should consider possible harms that could arise when the technology is being used as intended and functioning correctly, harms that could arise when the technology is being used as intended but gives incorrect results, and harms following from (intentional or unintentional) misuse of the technology. • If there are negative societal impacts, the authors could also discuss possible mitigation strategies (e.g., gated release of models, providing defenses in addition to attacks, mechanisms for monitoring misuse, mechanisms to monitor how a system learns from feedback over time, improving the efficiency and accessibility of ML). 
safeguards: 
 Question: Does the paper describe safeguards that have been put in place for responsible release of data or models that have a high risk for misuse (e.g., pretrained language models, image generators, or scraped datasets)?Answer: [NA] Justification: We propose a theoretical model of the neural circuit in the brain and it will not impose any risk.Guidelines:• The answer NA means that the paper poses no such risks.• Released models that have a high risk for misuse or dual-use should be released with necessary safeguards to allow for controlled use of the model, for example by requiring that users adhere to usage guidelines or restrictions to access the model or implementing safety filters. • Datasets that have been scraped from the Internet could pose safety risks. The authors should describe how they avoided releasing unsafe images. • We recognize that providing effective safeguards is challenging, and many papers do not require this, but we encourage authors to take this into account and make a best faith effort.12. Licenses for existing assets Question: Are the creators or original owners of assets (e.g., code, data, models), used in the paper, properly credited and are the license and terms of use explicitly mentioned and properly respected?Answer: [Yes]Justification: We adapt two figures from a published experimental paper (Fig. 1G-H ) and we clearly state where the figures are from.Guidelines:• The answer NA means that the paper does not use existing assets.• The authors should cite the original paper that produced the code package or dataset.• The authors should state which version of the asset is used and, if possible, include a URL. • The name of the license (e.g., CC-BY 4.0) should be included for each asset.• For scraped data from a particular source (e.g., website), the copyright and terms of service of that source should be provided. • If assets are released, the license, copyright information, and terms of use in the package should be provided. For popular datasets, paperswithcode.com/datasets has curated licenses for some datasets. Their licensing guide can help determine the license of a dataset. • For existing datasets that are re-packaged, both the original license and the license of the derived asset (if it has changed) should be provided. • If this information is not available online, the authors are encouraged to reach out to the asset's creators. 
new assets: 
 Question: Are new assets introduced in the paper well documented and is the documentation provided alongside the assets?Answer: [NA] Justification: This is a theoretical neuroscience study where the math formulation of the neural circuit model and its theoretical analysis presented in the paper are our main results.Guidelines:• The answer NA means that the paper does not release new assets.• Researchers should communicate the details of the dataset/code/model as part of their submissions via structured templates. This includes details about training, license, limitations, etc. • The paper should discuss whether and how consent was obtained from people whose asset is used. • At submission time, remember to anonymize your assets (if applicable). You can either create an anonymized URL or include an anonymized zip file. 
crowdsourcing and research with human subjects: 
 Question: For crowdsourcing experiments and research with human subjects, does the paper include the full text of instructions given to participants and screenshots, if applicable, as well as details about compensation (if any)?Answer: [NA] Justification: This is a theoretical neuroscience study.Guidelines:• The answer NA means that the paper does not involve crowdsourcing nor research with human subjects. • Including this information in the supplemental material is fine, but if the main contribution of the paper involves human subjects, then as much detail as possible should be included in the main paper. • According to the NeurIPS Code of Ethics, workers involved in data collection, curation, or other labor should be paid at least the minimum wage in the country of the data collector. 
institutional review board (irb) approvals or equivalent for research with human subjects: 
 Question: Does the paper describe potential risks incurred by study participants, whether such risks were disclosed to the subjects, and whether Institutional Review Board (IRB) approvals (or an equivalent approval/review based on the requirements of your country or institution) were obtained?Answer: [NA] Justification: This is a theoretical neuroscience study.Guidelines:• The answer NA means that the paper does not involve crowdsourcing nor research with human subjects.• Depending on the country in which research is conducted, IRB approval (or equivalent) may be required for any human subjects research. If you obtained IRB approval, you should clearly state this in the paper. • We recognize that the procedures for this may vary significantly between institutions and locations, and we expect authors to adhere to the NeurIPS Code of Ethics and the guidelines for their institution. • For initial submissions, do not include any information that would break anonymity (if applicable), such as the institution conducting the review. 
 
 Appendix: Justification: We confirm our research is abide by the Code of Ethics.Guidelines:• The answer NA means that the authors have not reviewed the NeurIPS Code of Ethics.• If the authors answer No, they should explain the special circumstances that require a deviation from the Code of Ethics. • The authors should make sure to preserve anonymity (e.g., if there is a special consideration due to laws or regulations in their jurisdiction). 


Decision: This paper considers how canonical cortical circuits with multiple interneuron types can implement sampling algorithms. The reviewers broadly agreed a strength of this paper is in its tight connection between the physiology and the sampling algorithm, which will be of interest to the Bayesian inference community. In light of this, I recommend acceptance.

There are some concerns that the paper only presents simple posteriors. I recommend the authors utilize the extra page allotted for accepted papers to include results with more complex priors.

Reviews: review: paper_summary: paper_summary: This paper studied how the introduction of two inhibitory neuron population affects Bayesian inference in firing rate models with additive noise terms. Analytical derivation and simulations are done with circular 1D input variable (orientation). They found the inclusion of SOM leads to faster inference.
strength_weakness: strength_weakness: 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.
questions: questions: 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)
limitations: limitations: Limitations are adequately addressed.
paper_summary: paper_summary: The paper proposes a dynamical model for how neurophysiologically realistic circuits of pyramidal excitatory neurons, with two types of inhibitory interneurons, could perform Bayesian inference in a generative model of a uniform prior (over a bounded domain such as orientation angles) and a Gaussian likelihood. Analyzing the proposed neuronal circuit, the paper finds that its dynamics implement Hamiltonian Monte Carlo, and that removing one type of inhibitory interneuron results in Langevin Monte Carlo.  The paper then claims to replicate certain neurophysiological findings, specifically regarding L2/3 of the laminar circuit but also, according to some experiments, of L5.
strength_weakness: strength_weakness: 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.  However, the authors have addressed this weakness in their responses and demonstrated that their model can sample from multimodal distributions.
questions: questions: 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?
limitations: limitations: 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.
paper_summary: paper_summary: The paper is build on previous work on Bayesian sampling to provide incremental but novel analytical insights about a neuronal circuit implementing bayesian inference. The circuit comprises excitatory neurons and two types of inhibitory neurons, PV and SOM, and uses the ring architecture. Authors find that the dynamics of such a circuit performs a Bayesian sampling algorithm and analyse analytically the properties of such sampling. They find that the circuit with E neurons and PV interneurons performs Langevin sampling while including SOM interneurons, the circuit performs Hamiltonian sampling.
strength_weakness: strength_weakness: The biggest weakness of the paper is its convoluted style of presentation - the paper is hard to read. Authors could dedicate more space to give intuitive explanations about the meaning of the main results. For example, authors could give some intuition about Langevin and Hamiltonian sampling and why having different types of interneurons gives different types of sampling. The meaning of the main results of the paper remains largely unclear. 

The authors emphasise that their results are relevant to a canonical circuit and describe the canonical circuit comprising three types of interneurons, however, their model only incorporates two types of interneurons. I understand this might be a work in progress. Nevertheless, I find it would be better to avoid mentioning the "canonical circuit" if the theory does not actually describe the canonical circuit.

Further, authors use the assumption that PV neurons are not tuned to the stimulus and justify this with findings of one empirical study. It seems rather that tuning of PV neurons is an open question, with some studies reporting relatively strong tuning of PV interneurons in the cortex (see e.g. Runyan et al., Neuron 2010,  Moore and Wehr, J. Neurosci. 2013). While it is justifiable to use the assumption of untuned PV interneurons, I find it better to present it as an assumption supported by some of the empirical literature.
questions: questions: Unless I misunderstood something, the model is by design limited to encoding of a single stimulus feature. An increasing amount of literature in neuroscience shows that biological neural networks (in particular in the cortex), typically represent several stimulus- and behaviour-related variables (Or at least, such variables can be reliably decoded from neural activity in vivo; see for example Stringer et al. Nature 2019, Stronger et al. Science 2019, Mussall et al. Nature 2019). These observations include the same neuron being responsive to multiple features of the stimulus or of the behaviour. Can authors comment on that? 

Does Bayesian sampling implemented by the model imply a static stimulus and a convergence of network activity to such stimulus? How does the network deal with a stimulus that has rapid shifts in dynamics or that changes the mean?

Can the model deal with non-Gaussian stimuli?

Tuning curves of neurons of the same cell type are homogeneous  across neurons (have the same shape). Is this assumption required for analytical computations? Could there be a computational reason that would instead allow heterogeneity of tuning across neurons of the same cell type?
limitations: limitations: Limitations are adequately addressed.
paper_summary: paper_summary: This paper investigates how canonical recurrent circuits in the cortex can implement sampling algorithms. The authors show that circuits with only E and PV neurons implement Langevin dynamics, while including SOM neurons enables a mixture of Langevin and Hamiltonian sampling.
strength_weakness: strength_weakness: 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).

Minor:
1. Figure 3 caption does not include panel F.
2. 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?).
questions: questions: I am curious about how the authors would construct a distributed version of the proposed circuit that can sample from high-dimensional stimulus posterior.
limitations: limitations: See weaknesses

Comments: title: Cortical regardance of our circuit model
Thanks!
Calculating the equilibrium distribution via Fokker-Planck approach
Dynamical system theoretical analysis, and our contribution
Invariant equilibirum distribution with mixing Langevin and Hamiltonian
Thank you
Follow-up question
comment: Thanks very much for your positive reply! 
We will include the new results in the main text of the revised manuscript.
I thank the Authors for their thorough reply.
The added results are important to be included in the main text. Assuming they will be included, I increase my score.
Thanks for your reply. We will improve the manuscript further by adding more background about the canonical circuit, interneurons, and laminar layers, in order to make our manuscript more friendly to machine learning readers.

To explain the cortical regardance of our circuit model, let's talk about some background of the cerebral cortex. If we flatten the cerebral cortex, it is a 3D neural sheet with $x-y$ axes denoting the positions on the cortical surface, and the $z$ axis denoting the depth perpendicular to the cortical surface. If we move along the $x-y$ directions, we have different cortical regions, whereas moving along the $z$ direction, we go through 6 laminar layers. 

For the __cortical column__, it is regarded as the cylinder structure perpendicular to the cortical surface, i.e., the set of neurons spanning the whole $z$ cortical axis but only a tiny part along the $x-y$ axes. Typically, the neurons within a cortical column encode a similar stimulus feature, e.g., the neuronal column preferring a $0$ degree of moving direction. Typically there are hundreds of excitatory neurons (E) within a column. \
Furthermore, a __cortical hypercolumn__ is regarded as a set of cortical columns (spanning more on $x-y$ axes compared with one column) containing a full set of preferred stimulus feature values, e.g., all columns preferring all moving directions. 
The size of the hypercolumn in the primary visual cortex is about $1$ mm $\times$ $1$ mm along the $x-y$ axes. \
At last, a __cortical region__ is composed of many cortical hypercolumns. For example, each hypercolumn processes the visual inputs at a particular visual location by having all neurons prefer different moving directions.

Therefore, our model is regarded as the microcircuit within one laminar layer of a cortical hypercolumn, rather than a model including interactions across laminar layers.
We mostly regard our model as the layer 2/3 within a hypercolumn in the primary visual cortex, since most of the experimental evidence comes from there. 
Specifically, $u_E(\theta)$ in Eq. (1) in the paper regards as the excitatory (E) neurons in layer 2/3 of a cortical column preferring direction $\theta$, which mathematically captures the mean response of all E neurons there. Then $\{ u_E(\theta)\}_\theta$ captures the responses of E neurons of the whole hypercolumn.
My thanks to the authors for their additional analyses and results.  I agree with them that the bottom-up, neuroscience-first approach is important to this paper, and so am satisfied to have the additional results included in supplementary material for the final manuscript.

However, a last request: could the authors please clarify throughout the manuscript when they are referring to a canonical intralaminar microcircuit vs the more typical usage of a canonical laminar (eg: across the column) microcircuit?  This clarification was important for my understanding of the whole paper and needs to come across in the text.

From there I can see through to raising my score.
Thanks very much for your reply! 

> About the last point on heterogeneous tuning curves leading to chaotic dynamics, do I understand correctly that having heterogeneous tuning curves in the Eq.(1) instead of homogeneous ones, the stochastic term in the Eq.(1) could potentially be removed or replaced by a deterministic term?

Yes, this is exactly what we suggested.

In computer simulations of sampling algorithms, we call a built-in __pseudo-random__ function whose underlying dynamics is also __chaotic__. Conceptually, this is similar to the chaotic E/I balanced network dynamics. 
In simulating a neural sampling network model, although we also call a random function and inject the variability into the model, we are concerned that real neural circuits probably do not have a statistically stable source of variability like the random function in programming. And it is less likely a stable variability can magically or trivially emerge in neural circuits.
Therefore we think how the noise (or variability) used in sampling can be generated in neural circuits is an important issue at the neural implementation level, but which has not been paid sufficient attention in the field. 

Due to the page limit, we don't have sufficient space to emphasize the importance of the multiplicative variability (Eq. 1 last term) in our study. Please find our below reply.

> Unless I overlooked something, I find that the variable \xi (\theta,t) in the Eq.(1) is not clearly defined. Also, how does the type of stochastic process in the last term of Eq.(1) influence the type of sampling resulting from the network?

Sorry for the confusion. The $\xi(\theta, t)$ is a standard Gaussian white noise with zero mean and unit variance, and satisfying 
$\langle \xi(\theta, t)\xi(\theta', t') \rangle = \delta(\theta - \theta') \delta(t - t')$. 
In the rate-based network model, we use the (continuous) Gaussian white noise to approximate the multiplicative Poisson variability that results from stochastic spike generation (lines 83-87). This approximation is reflected by that the variance of $\xi(\theta,t)$ is proportional to instantaneous synaptic input $u_E(\theta,t)$ (negatively rectified with $[\cdot]_+$ in Eq. 1).

This __multiplicative__ variability is essential for the network model with __fixed weights__ to flexibly sample posteriors with different uncertainties (Fig. 2G, 3F and Eq. 21; Please check the significance of this function in Global Rebuttal and our reply to Reviewer __RifV__). That is, if we replace the multiplicative noise with an additive noise, e.g., $\sigma \xi(\theta, t)$, our network model with fixed weights cannot sample posteriors with different uncertainties.

To intuitively explain the underlying mechanism, let's review the Langevin sampling dynamics (a copy of Eq. S33)
$$\frac{dz}{dt} = \tau_L^{-1} \nabla \log p(z) + (\tau_L/2)^{-1/2} $$
which requires the __drift and diffusion terms both share the same factor $\tau_L$__. How this stringent requirement could be achieved in a recurrent neural circuit model has not been investigated before, to our best knowledge. 
With the multiplicative noise in Eq. (1), the samples generated by $E$ neurons are governed by Eq. (13) (copied below),

$$ \frac{dz_E}{dt} = \tau_E^{-1} g_{EF}(\mu_z - z_E) + \sigma_E \tau_E ^{-1/2} \xi_t $$

The common factor in drift and diffusion terms is $\tau_E = \tau U_E$ (defined in line 185), which is proportional to $U_E$, the magnitude of the synaptic inputs of E neurons (defined in Eq. 9). 
Intuitively, this is because (the full math derivations are in Eqs. S23-S30)

- Higher neuronal response (larger $U_E$) will have larger inertia, analogous to a heavier object having larger inertia in kinematics, which is reflected by a large value of $\tau_E$ in the drift term.
- Due to the multiplicative noise (Eq. 1, last ter,), larger $U_E$ will incur larger nosies in single neurons. Then the projected noise on the $z$ dynamics (the low-dim attractor manifold) will be higher. 

Combined, the multiplicative noise effectively links the drift and diffusion coefficients in the Langevin dynamics to share a common factor.
> 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.

This is a very important question we are interested in exploring. Regardless of the biological plausibility, one possibility is introducing a __gating mechanism__ on the recurrent weights between subnetworks, similar to gated recurrent networks. Given a value of the gating variable, the network stores a correlation prior. Combining all values of the gating variable, the network can store a __mixture__ of correlation priors, which may emerge in high-order prior statistics. This is a possibility we would like to explore in the near future.

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

Thanks very much for your positive reply!
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.
I thank the Authors for their reply. 

About the last point on heterogeneous tuning curves leading to chaotic dynamics, do I understand correctly that having heterogeneous tuning curves in the Eq.(1) instead of homogeneous ones, the stochastic term in the Eq.(1) could potentially be removed or replaced by a deterministic term?

Unless I overlooked something, I find that the variable \xi (\theta,t) in the Eq.(1) is not clearly defined. Also, how does the type of stochastic process in the last term of Eq.(1) influence the type of sampling resulting from the network?
I figured it out, this can be proved using the Fokker-Planck equation. It can be shown that the joint stationary distribution is $p(y,z)\propto \exp(-\frac{y^2}{2} + U(z))$. Scale separation method that the authors gives critically depends on the assumption that y is faster than z.
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.
>$\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.
Dear reviewer:

I just realized the display of Markdown at OpenReview is different from my laptop and the math equation numbers disappeared in OpenReview, which may confuse you about which equation is Eq. (3) in our last reply. Hence we revised our math derivation in the last reply and we hope the current derivation addresses your question about why the $p(z)$ is unchanged. If you feel our derivation by separating time scales of $z$ and $y$ is limited, we can provide a more rigorous version.

In summary, by using time scale separation of $z$ and $y$, the $z$ dynamics in the Hamiltonian sampling is effectively a Langevin sampling dynamics of $p(z)$. Then by mixing two Langevin dynamics sampling the same equilibrium distribution, the equilibrium distribution remains unchanged. 

Yes, in our derivation $\nabla \log p(z)$ is general and can be nonlinear of $z$.
>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).
#### 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,
$$
\begin{align}
\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)
\end{align}
$$
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,
$$
\begin{align}
\dot{z} = - \tau_z^{-1}\beta^{-1} \nabla U(z) + (\tau_z \tau_y)^{-1/2} \xi_t \quad\quad (3)
\end{align}
$$
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)$.
Thanks for the authors' further reply, the last sentence is what I was asking, I said in the original comments that the linear case is straightforward. If it only works for the linear case, that brings another problem, as $\nabla \log p(z)$ is linear iff p(z) is Gaussian, does that mean any distribution that the circuit can sample from is at best mixture of Gaussian? Natural stimulus is known for its high-order statistics, while Gaussian is completely described by its first two moments. Although I don't think one should a priori assume the capability of the neural circuits (that's the job of experimentalists), I am very skeptical that neural circuits can only sample from distributions of Gaussian nature. Why consider nonlinear neural circuits when you just need to sample from Gaussian? 

I also read the authors' rebuttal to all reviewers
> Studied a nonlinear, biologically plausible dynamics but without analytical results of the nonlinear dynamics (e.g., Benjamin, ICLR 2024; Chen, NeurIPS 2023; Echeveste, Nat. Neurosci., 2020). 

This is not true, Chen, NeurIPS 2023 and Echeveste, Nat. Neurosci., 2020 both have substantial analytical results. The diffusion model used in Benjamin, ICLR 2024 is also heavily backed by analytical results. Maybe the authors can be more specific on what their contribution is? For example, the Recurrent connection kernels proposed is also studied in e.g. Dong et al., 2022, where it is called CANN, what is the difference between CANN and the recurrent connection kernels?
Thanks for the reviewer's reply! 

In Ref. 14, Eqs. 14-17 stated mixing the Langevin and Hamiltonian sampling dynamics will not perturb the stationary distribution, although they didn't provide rigorous math proof. 
Unlike Ref. 14 which considered the classical Hamiltonian dynamics, the present study considered a modified Hamiltonian dynamics with friction whose mathematical form is the same as the Eq. 9 in Ref. 50. Then we mix this modified Hamiltonian dynamics with the Langevin dynamics, which yields the Eqs. S34-S35.. 

We are sure mixing two __linear__ SDEs with the same equilibrium distribution will not change the equilibrium distribution, which can be easily proved by writing down the Fokker-Planck equation of the SDE. Nevertheless, we are not sure whether this is still held for nonlinear SDEs.
Could the author be more specific about where in Ref. 14 and 50 this mixture of Langevin and Hamiltonian is used?