Shaping the distribution of neural responses with interneurons in a recurrent circuit model
Authors: David Lipshutz
Eero P Simoncelli
Venue: NeurIPS 2024 poster (2024)
Keywords: Efficient coding
optimal transport
Gaussianization
Hebbian plasticity
gain modulation
neural adaptation
Primary area: neuroscience_and_cognitive_science
TL;DR: Motivated by the efficient coding hypothesis, we propose a normative circuit model with local interneurons that learns to reshape its inputs so that its responses follow a target density.

Abstract: Efficient coding theory posits that sensory circuits transform natural signals into neural representations that maximize information transmission subject to resource constraints. Local interneurons are thought to play an important role in these transformations, shaping patterns of circuit activity to facilitate and direct information flow. However, the relationship between these coordinated, nonlinear, circuit-level transformations and the properties of interneurons (e.g., connectivity, activation functions) remains unknown. Here, we propose a normative computational model that establishes such a relationship. Our model is derived from an optimal transport objective that conceptualizes the circuit's input-response function as transforming the inputs to achieve a target response distribution. The circuit, which is comprised of primary neurons that are recurrently connected to a set of local interneurons, continuously optimizes this objective by dynamically adjusting both the synaptic connections between neurons as well as the interneuron activation functions. In an application motivated by redundancy reduction theory, we demonstrate that when the inputs are natural image statistics and the target distribution is a spherical Gaussian, the circuit learns a nonlinear transformation that significantly reduces statistical dependencies in neural responses. Overall, our results provide a framework in which the distribution of circuit responses is systematically and nonlinearly controlled by adjustment of interneuron connectivity and activation functions.

Shaping the distribution of neural responses with interneurons in a recurrent circuit model
David Lipshutz
Center for Computational Neuroscience, Flatiron Institute diipehutzeflatiioninniiiiii org
Eero Simoncelli
Center for Computational Neuroscience, Flatiron Institute Center Neural Science, New York University eero. simoncelliisseeomom........ edu
Abstract
Efficient coding theory posits that sensory circuits transform natural signals into neural representations that maximize information transmission subject to resource constraints. Local interneurons are thought to play an important role in these transformations, dynamically shaping patterns of local circuit activity to facilitate and direct information flow. However, the relationship between these coordinated, nonlinear, circuit-level transformations and the properties of interneurons (e.g., connectivity, activation functions, response dynamics) remains unknown. Here, we propose normative computational model that establishes such relationship. Our model is derived from an optimal transport objective that conceptualizes the circuit's input-response function as transforming the inputs to achieve an efficient target response distribution. The circuit, which is comprised of primary neurons that recurrently connected to set of local interneurons, continuously optimizes this objective by dynamically adjusting both the synaptic connections between neurons well as the interneuron activation functions. In an example application motivated by redundancy reduction, we construct circuit that learns a dynamical nonlinear transformation that maps natural image data to spherical Gaussian, significantly reducing statistical dependencies in neural responses. Overall, our results provide framework in which the distribution of circuit responses is systematically and nonlinearly controlled by adjustment of interneuron connectivity and activation functions.
Introduction
The problem of transforming a signal into representation with a given target distribution (or within target set of distributions) is classical problem whose origins can be traced back more than two centuries [1]. Many methods statistics, signal processing and machine learning can be interpreted within this context. For example, data whitening is a common preprocessing step that linearly transforms a signal to have identity covariance [2]. Independent component analysis [ICA; 3] is signal processing method that linearly transforms signal
minimize higher-order statistical dependencies in addition to removing second-order dependencies. Nonlinear transformations of skewed or heavy-tailed data approximately Gaussianize their distribution can facilitate statistical analyses [4, 5]. Machine learning methods for density estimation such as Gaussianization [6-10] and normalizing flows [11-13] nonlinearly transform nigh-dimensional signals with complex densities into more tractable representations with approximately Gaussian densities.
38th Conference on Neural Information Processing Systems (NcurIPS 2024).

---

T12
up of aduptive interneuron
92
W1
W1
Signal: (81,
Network
Response: (11,
Wg
Wg
Figure I: Schematic of recurrent circuit with N 2 primary neurons and K 3 interneurons. Left: (81.82) with Ps- Center: Primary (r), receive external feedforward inputs, s, recurrent feedback from an
Scatter plot of a 2D input signals (gray points) circles), with outputs auxiliary population of interneurons (purple circles). -Wn, where n (n1, ni, Tiy) are the internurron outputs. Projection vectors { [W1, Wa, Wa encode feedforward synaptic weights connecting primary neurons to interneurons i 1, 2, with symmetric feedback connections. Inset: The 2th interneuron (here = 1) receives weighted inputs 2; = r. wi, which is fed through the activation function / [0,,, .) and scaled by the gain gi to generate the output == g: f(,,, z.). Right: Scatter plot of the 2D circuit responses r (r1,T2) with Purget-
(black
This problem may also provide framework for constructing and evaluating normative theories of sensory processing. Efficient coding theory posits that sensory systems maximize the information they transmit about sensory signals downstream areas subject resource constraints [14-17]. In one instantiation of this theory, the redundancy reduction hypothesis posíts that sensory circuits transform natural signals into representations minimize eliminate statistical dependencies between coordinates, essentially producing factorized response distributions [14, 15, 18]. In separate, but related, instantiation, sparse coding theory posits that population responses are optimized for sparsity [19-21], which is naturally interpreted as a constraint on the shape of the distribution of responses, In another line of theoretical work, sensory representations are posited maximize the Fisher information of the inputs [22-24], which can be interpreted as statement about the joint distribution of the inputs and responses. Importantly, each of these theories can be formulated as a transformation of the signal into a representation with target distribution that is optimal under information theoretic and metabolic constraints. However, not clear how neural circuits learn or implement these potentially nonlinear transformations.
Neural circuits are typically comprised of populations of primary (excitatory) neurons and local (inhibitory) interneurons. Extensive experimental measurements have led to the idea that local interneuron populations allow neural circuits to flexibly shape patterns of primary neuron responses so as to coordinate or regulate information flow [25-31]. Consequently, local interneurons are natural candidate mechanisms responsible for shaping circuit responses into efficient representations. However, the precise relationship between the physiological properties and connectivity of local interneurons and the coordinated response properties of populations of primary neurons remains unclear.
Several normative mechanistic models have been proposed explain how neural circuits linearly transform their inputs into representation whose distribution lies within a target set (e.g., the set of distributions with identity covariance) [32-36]. These models are derived from optimization objectives for linear redundancy reduction, including (adaptive) decorrelation and ICA, and the circuit parameters (e.g., synaptic weights, response gains) are optimized to match the data distribution. In these cases, the optimization steps correspond to processes such as gain modulation and synaptic plasticity, thus demonstrating how adjustments of circuit parameters according to local signals can optimize global, circuit-level objective for redundancy reduction of the neural responses. While linear transformations can remove second-order statistical dependencies, they cannot remove higher-order statistical dependencies that are prominent in sensory signals [7]. Furthermore, early sensory systems exhibit a host of prominent nonlinear response properties [37-39] that are not well-approximated by linear models.

---

Here, we seek a normative circuit model that nonlinearly transforms inputs to produce responses with a (spherical) target distribution. We develop an optimal transport objective for transforming the input signal into neural representation with a given target distribution, and derive an algorithm (Alg. 1) that can be mapped onto a dynamical model of a neural circuit, Fig. 1. The circuit model is comprised of primary neurons that are recurrently connected to local interneurons and the circuit adapts its responses using a combination of Hebbian synaptic plasticity and interneuron adaptation. Complementary computational roles for Hebbian plasticity and interneuron adaptation emerge from this analysis: (i) synapses are updated according to a Hebbian learning rule to identify projections of the signal that are least aligned with the target distribution (essentially, form of projection pursuit [40]); (ii) gains and nonlinear activation functions of interneurons are adjusted to transform the marginal circuit responses along directions defined by the synaptic weights. Together these operations transform the distribution of responses to approximate the target distribution.
As primary test case motivated by redundancy reduction theory, we apply our algorithm to inputs derived from natural images (responses of local oriented filters, qualitatively similar to those found in the primary visual cortex) and the target distribution is the spherical Gaussian.' We find that the algorithm learns a nonlinear transformation that approximately Gaussianizes the responses and significantly reduces statistical dependencies between coordinates. Overall, our results demonstrate how local interneurons may adjust their connectivity and response properties to nonlinearly reshape the distribution of circuit responses, thus facilitating efficient transmission of information.
2 Circuit objective
Consider neural circuit with N 2 2 primary neurons that transforms input signals E R.V which are distributed according to a density Ps. into circuit responses r E RR (Fig. 1). The inputs may represent direct sensory input (c.g.. the rate at which photons are absorbed by cone) or the weighted sum of multiple inputs (c.g., the postsynaptic current). The responses r represent the firing rate (or the logarithm of the firing rate) of the neuron. For simplicity, we assume that the circuit responses are a deterministic function of the input signals; that is, r = T(s) for a function T: RR + R
2.1 Optimal transport objective
We assume that the objective of the circuit is to transform its inputs s so that the circuit responses T's) follow (spherical) target distribution Ptarget while minimizing the Li=ddsstnce between responses and input signals. Mathematically, this corresponds to an optimal transport problem [42]
min E [I] - s12 + (|] such that T(s) ~ Pharget: T
where the minimization is over suitable class of functions T, A E E is a regularizing term and, unless otherwise noted, expectations are over the input distribution Ps. Note that the choice of X does not affect the optimal solution as E[| T(s) is fixed provided T(s) Ptarget' however, it will affect the optimization algorithm. Assuming Ps is sufficiently regular, the minimum in eq. (1) is the squared W distance between Ps and Prarget and the input-response transformation T is the so-called optimal transport plan that maps the input distribution ys to the target distribution Ptarget-
Our goal is to derive an online algorithm for optimizing the objective in eq. (1) that can be implemented in neural circuit model. We accomplish this by definng a distance between the response distribution Pr and the target distribution Puager that can be estimated in neural circuit, and then solving the optimization problem by incorporating this measure as a constraint, using the method of Lagrange multipliers.
2.2 Measuring the discrepancy between the response distribution and the target distribution
Common candidates for measuring the discrepancy between two distributions include Kullback- Leibler (KL) divergence integral probability metries [43]. The former typically require numerous
Biolooiccall,, we can interpret this as a Gaussian distribution of voltages on cell bodies, which are transformed through an exponential activation function to log- lormall distribution of firing rates that are vositive-valued, heavy-tailed, and similar to those observed in cortex [41].
3

---

samples to estimate the density Pr, whereas neural circuits must operate in the online setting without access to the full history of their responses. Therefore, we use an integral probability metric that quantifies the difference between the random variables when evaluated using constraint functions that can be rapidly estimated online.
We restrict our solution to use constraint functions that compare the marginal response distributions to the marginals of Ptarget denoted Pmarginal. which are all equal under our assumption that Parget is spherical. Such constraint functiions are well matched to signals that are linear mixtures of independent sources, as in generative models for ICA. Even when the signal statistics are not generated according to linear mixture model, the Cramér and Wold theorem [44] suggests that transforming sufficiently many marginals of the response distribution may effectively transform the multivariate response distribution (though the number of marginals required may be quite large). Our motivation for measuring marginal response distributions is due, in part, by our goal of modeling local interneurons, whose inputs are naturally modeled weighted sums of primary neuron responses (i.e., their input distributions are marginals of the primary responses).
To compare the marginal response distributions, we first select a finite set of directions defined by K I unit vectors wi,..., WK RN,, which can be either randomly sampled, chosen based on prior knowledge of the signal statistics, or learned from data using projection pursuit [40]. We then choose a class of scalar functions {h(0,-) -)] parameterized by which defines semi-metric between the marginal of in the direction w and Pmarginal to be max, [E(((f, r w)]l, where
o(0, == h(6,2) Ea-erepta [h(O, 2)].
For example, when {h(0,-)} parameterizes all Lipschitz-1 (indicator) functions, this induces the Wasserstein-1 (total variation) distance between the marginal response distribution and Pmanginal-
Given directions Wi,..., WK and constraint functions {h((f,)) })}, we express the distance between the response density Pr and the standard Gaussian distribution the sum of the marginal distances:
K
dw dw.)
Σ
max Err,, 10(0,, w.)] l,
where W [Wi,..., WK] is the N X K matrix of concatenated unit vectors."
How do we choose the constraint functions (h(0..)) In general, the choice should be well-suited to the input distribution pa and the target distribution Pr. For example, consider the simple case when the inputs follow centered Gaussian distribution with unknown covariance structure, and the target distribution is the spherical Gaussian distribution . I). The Wasserstein-2 distance between marginal distribution of the inputs and the standard normal distribution ((, 1) can be expressed in terms of the difference between the second moments, so a quadratic function h=2) _ suffices. In section 4.2, we consider a parametric class motivated by natural signal statistics.
2.3 Optimization using Lagrange multipliers
We replace the condition r ~ Pragett in eq. (1) with the condition maxw dw.s(p) ) = 0, which we enforce using Lagrange multipliers. This results in the minimax optimization problem
max max max E min L(W, o.g.s.r). W
(2)
0
where C is defined by
L(W (WW.Ogg.sr)) = ||| s|1 + |||r|| Σ grolOi,r w:),
and 0 : (#1...., 0x) is the set of concatenated parameters, dimensional vector of Lagrange multipliers, and the circuit transform is defined by T(s) arg min, L(W, 0. g, r). The maximization over g and 0 effectively minimizes the distances between the marginal distributions of the responses r along the directions WI,..., Ww and Pmarginal whereas the maxinization over the matrix W learns directions along which the marginals of s are least aligned with Pmarginal: essentially performing projection pursuit [40].
(01...., 9k) is K-
"This distance is closely related to (max )sliced Wasserstein metries [45, 46], which quantify the distance between two distributions in terms of the Wasserstein distances between their marginal distributions. Specifically, when K = 1 and h(0,, -) parameterizes the set of all Lipschitz-I functions, then the sliced and max-sliced Wasserstein- 1 distances between pr and Parget are Bu- gN 1) j/dd.. (Pr)] and max dw a (pr).

---

3 Algorithm and circuit implementation
We now derive an online maadient-based algorithm for optimizing the objective in eq. (2), then map the algorithm onto a recurrent neural circuit. Spiking activity operates on a much faster timescale than neural synaptic adaptation mechanisms, so assume that the neural activities equilibrate before the neural activations and synapses are updated.
3.1 Fast recurrent neural dynamics
At each iteration, the circuit receives a stimulus The (discretized) recurrent neural response dynamics (Fig. 1) correspond gradient-descent minimization of L with respect to r:
r r-r+ Yr s pr E new:
n4 = (0,,
represents a leak term, 2; : W; the weighted input to the qth interneuron, f(0,,,) 0o(0;,1 -/0z is the activation function, 9i is a multiplicative gain that scales the output, and n, denotes the output. Notably, the interneuron activation response function f (0,. is parameterized by (g., 0,), which can vary across internuroons, so the interneuron responses are heterogeneous. For each synaptic weights w, connect the primary neurons the tt interneuron and symmetric weights w, connect the 7th intermeuron to the primary neurons. From eq. (4), we see that the neural responses are driven by the signal s, a leak term -r, and recurrent weighted feedback from the interneurons - Wn, where := (MI,...,mkk) nk)).
where 7,
a small constant, f =
Since the neural activities equilibrate before other updates performed, the responses point of C(W, 0, g. 8, .). In general, we do not have a closed-form expression for r; however, if gi 2 0 and ((0,, -) are convex, then L(W, 0, g, S, convex and we can express the transform T(s) = arg min L(W, s,g.s,r).
fixed
(5)
In Appx. A we show that the transform T(-) is invertible whenever 9: 2 0 and (( .) are convex and thus it defines precise relationship between the input distribution Ps and response distribution pr-
3.2 Gain modulation, activation function adaptation and Hebbian plasticity
After the neural activities reach equilibrium, we maximize by taking concurrent gradient-ascent steps with respect to 9t. 0, and w,:
Ag. = 17,6000 zz),
A0, = Veb =4);
Aw, T,,,
where ng 77. nur O are the respective learning rates, which control the relative speeds of gain modulation, neural adaptation and synaptic plasticity, respectively. For example, at the extremes, can fix the gains, activation functions synaptic weights by setting fly 0, r/o 0 The respectively. Notably, while synaptic plasticity [47] and gain modulation [48] are well-studied circuit mechanisms that support learning and adaptation, adjustments nonlinear neural activation functions are not as well established. Nevertheless, there is some emerging evidence that neurons also adapt their activation functions response to changes in their input statistics [49, 50].
The circuit operates online and the updates are local in the sense that the updates to the gain and activation function of the th interneuron depend only variables 8; and 2,. The updates to the synapses w, (and -w;) are proportional (inversely proportional) to the product of the pre- and postsynaptic activities, so they are both local and Hebbian (anti-Hebbian) [51]. Finally, to ensure that the vectors Wi,..., wR have unit norm, we normalize the weights after each update: W; f- w:/ll w:ll- This can be viewed as form of homeostatic plasticity such synaptic scaling [52]. While the feedforward weights w, and feedback weights -w; are constrained to be symmetric, these can be decoupled due the symmetry of the Hebbian learning rule. Both theoretical and empirical evidence of this is shown for a related adaptive whitening circuit in [36, appendix E.2].
3.3 Online algorithm
Combining the neural dynamics, the interneuron adaptation and synaptic plasticity steps yields our online algorithm (Alg. 1), which we write vector-matrix notation by defining the normalization function P(W) == [wi/|walll..... ww///ww/////l
5

---

Algorithm 1: Approximate optimal transport with Hebbian plasticity and interneuron adaptation
1: input: 81,82,--
2: initialize: W, A. g. 14, TJr fg. TJo The
3: for =1,2.... do
4: I; + St
5: while not converged do
6: Ze f W rr
// iitterneuron inputs
7: n. < go f z:):
// intternnuooon outputs // neural responses
8:
It ret+1 1), (st - pr) Wn.);
9:
end while
10:
g g + (0 0 2,);
// gain update
11:
④ 0 + no o z:);
// activation update
12:
W - P(W + T,
// Hebbian homeostatic plasticity
13: end for
3.4 Relation existing algorithms
Algorithm 1 is naturally viewed as nonlinear extension of a number of existing algorithms for linear data whitening that have neural circuit implementations [34, 36, 53, 54]. particular, when the constraint function is quadratic, h(z) $ and the target distribution is the spherical Gaussian, (O I). then the activation function is the identity, f(z 2, and the optimization in eq. (2) enforces that the second moments of the responses match the second moments of (0, I), which corresponds to data whitening. If the gains are fixed (7g = 0) and K 2 N (i.e., synaptic adaptation only), then Alg. 1 corresponds to the adaptive whiteming algorithm presented in [34, 53]. Alternatively, if the synaptic weights are fixed (7) = 0) and 2 N (NN 1)/2 (i.e., interneuron gain adaptation only), then Alg. I corresponds to the adaptive whitening algorithm presented in [54]. Finally, if the gains adapt on a fast timescale and the synapses update on a slow timescale (i.e., 10 >> Thu >> 0), Alg. corresponds to the multi-timescale adaptive whitening algorithm presented in [36].
When the target distribution the spherical Gaussian and W is constrained to be an orthogonal matrix, Alg. is related to existing iterative algorithms for Gaussianization that alternate between (a) orthogonal transformations and (b) marginal Gaussianization of the coordinates [6. 8). In the case that the column vectors of W are orthogonal, Gaussianization along one marginal does not affect the responses along other marginals, allowing these operations to be performed independently of one another. In general, we allow the column vectors of W to be non-orthogonal and potentially overcomplete so that marginal Gaussianization along one basis vector affects the other marginal distributions. Consequently, the algorithm is more complicated to analyze mathematically (e.g.. obtaining convergence guarantees), but more biologically realistic since neural systems are unlikely to have orthonormal synaptic weight vectors.
4
Gaussianization of natural image statistics
We apply our algorithm to the problem of efficient nonlinear encoding of natural signals, specifically oriented filter responses to visual images. Redundancy reduction theories posit that early sensory systems transform natural signal into neural representations with reduced statistical redundancies [14, 15, 17]. In support of this hypothesis, early sensory representations exhibit far less spatial and temporal correlation than natural signals [55, 56) and methods such as linear ICA have been used to derive optimal representations of natural signals that are approximately matched early sensory neural properties [18, 19, 57].
Linear whitening and ICA transforms can eliminate simple forms of statistical dependency, but their responses exhibit higher-order statistical dependencies when applied to natural signals [58]. suggesting that sensory systems more efficiently represent natural signals by implementing nonlinear transforms. Consistent with this, nonlinear phenomenologica models of neural responses (c.g., divisive normalization [37, 39]) effectively reduce these higher-order statistical dependencies [38, 59]. However, the circuit mechanisms that support these nonlinear transformation are unknown.
Example code for experiments be found at https: / h/ttt///gtttttttt/ com/ wwt.///owwtt//////////////totttt /ttpppsstuttt/hapigg

---

A) natural image
B] signal density Fr
c Internurennat. gko.d
D) circstt transfom em Nld
response dessity
0.45, 0.8
0.41 200
10a
20-4
10a
-2
-2-
10
10-3
10-1
109
105
10"
10
10
-2
C.D.2 0.42
0.01,
2
Figure 2: Marginal Gaussianization of local filter responses. A) Three example natural images from the Kodak dataset. B) Histograms of local filter responses (black lines) and fitted generalized Gaussian density (red dashed lines) with scale and shape exponent 8. C) Learned interneuron activations g z). with f(0, 2) defined eq. (6) and learned g and and D) corresponding stimulus-response transforms T(s). The optimal activations transforms are shown as thick gray curves. E) Histograms of circuit responses (black lines) and the Gaussian density (red dashed lines).
We consider the particular target of Gaussian responses, which may be interpreted firing rates, or as their logs (i.e... membrane potentials are Gaussian, and firing rates are exponentiated), From the perspective of coding efficiency and computation, Gaussian representations are appealing for a variety of reasons. First, among distributions with given covariance structure, Gaussian distribu- tions have maximum entropy. Therefore, if metabolic demands are a function of the varianne of the response distribution, then the Gaussian distribution maximizes information transmission under metabolic constraints. Second, compression based on information theoretic objectives often reduces to linear projection when the data distribution is Gaussian [60-62], so Gaussianization facilitates downstream computation. In addition, the efficiency of the representation preserved under orthogonal transformation [7]. Finally, experiments have shown that single neurons in the fly early visual system adaptively Gaussianize their univariate responses [63], and neural populations in early sensory systems decorrelate their multivariate responses [56, 64-67]. Therefore, neural circuit models that nonlinearly transform signals to jointly Gaussianize their responses may offer normative, parsimonious explanations of nonlinear transformations in early sensory systems.
4.1 Description of the input signal
We computed the responses of a local oriented filter [68], roughly matched to typical receptive field selectivity neurons primate visual cortex [69], applied to natural images from the Kodak dataset [70]. These local filter responses can be well approximated by generalized Gaussian distributions of the form (s(s) x exp(-1s (l//// P) where is referred to the scale parameter and is referred to as the shape parameter [71]. Fig. 2AB shows example images and histograms of the local filter responses along with fitted generalized Gaussian distributions whose scale and shape parameters (a, (3) vary across images. While linear methods such as variance normalization are sufficient for rescaling the distribution, adaptive nonlinear transformations required to reshape these heave- distributions.
notorious for their sparse heavy-tailed statistical properties that
Next, we generated 2-dimensional signals from pairs the local filter responses for images at fixed horizontal spatial offsets ranging between d 2 and 64. Contour plots (using kernel density

---

A) signal
B) ZCA
C) responses (K=2)
D) responses (K=3)
E) responses (K
Mi-1.33e-1
Mi-1 00e-1
91e2
MI-5 02e-3
Mi-4.6le-3
Ml-3.936-2
MI=4,0264
Mi
MI-5.6e-3
MI-4.41e33
Miel 2
1.
Mind Min4,90p-3
Min4.ble-3
Mim5.
d=32
Figure 3: Joint Gaussianization of pairs of filter responses, spatially displaced by 2, 8, 32 pixels. Evenly spaced contours for the spherical Gaussian distribution are depicted as dashed red circles. Iso-probability contours for (A) the local filter responses, (B) ZCA whitened local filter responses, and (CDE) learned circuit responses (with K 2, 3, 4 interneurons) depicted (in the respective column) black curves along with the estimated mutual information between coordinates. The learned column vectors of are indicated by the faint gray lines.
estimation) of the local filter response pairs and symmetric (or ZCA) whitened local filter response pairs for natural image (specifically, the top-left image from Fig. 2) are shown for select d in Fig. 3AB. To quantify the statistical dependencies between coefficients, we estimated the mutual information between the pairs of coefficients after discretizing them into bins of width 0.5. In Fig. 4, we plot the estimated mutual information (using bin size 0.5) between the local filter response pairs (blue line) and ZCA whitened local filter response pairs (orange line) for spatial offsets between d = and d = 64. Note that aside from d = 2. the linear ZCA whitening transform does not significantly reduce the mutual information between coordinates. (Similar results have been found when applying linear ICA transforms; see, e.g., [72. Figure 6].) Fig. also suggests that ZCA whitening can even slightly increase the mutual information between coordinates see also, rows d and d 32 of Fig. 3AB-though these effects are quite small and may be a consequence of the discretization step when estimating mutual information.
4.2
Choice of activation functions
How do we choose the family of activation functions {f(,,)) One approach is to choose kernel that can approximate general class of functions. An alternative approach, adopted here, which motivated by the efficient coding hypothesis [14, 15, 17], is to choose family of activation functions that well matched to the marginal statistics natural signals. Fig. 2C, we plot examples of optimal activations for transforming local filter responses from different images (thick gray curves).
Since the marginals of the local filter responses are well-approximated by generalized Gaussian distributions [71], a sensible approach is to identify family of intermeuron activation functions that are optimal for transforming generalized Gaussian distributions with varying (a, B) into the standard Gaussian distribution. When 0, this implies (see Appx. B) that for each choice of scale a and shape 8, there a gain g and parameter such that
9f = ((-),
FT]
where (--) is the edf of (0, 1). However, we define f(0, 2) in terms of the above display, then we do not have closed-form solution for o(0, z) Vod(O, 2). which are both required to implement

---

signal
whitened signal
responses (K=2)
responses (K 3)
responses (K 4)
matual In Indommat (bo)
64
spatial affset (pixels)
Figure 4: Estimated mutual information, with 95% confidence intervals (estimated across 23 images), for original signals (pairs of filter responses), the ZCA whitened signal, and learned circuit responses (with K 2,3, interneurons).
Alg. 1. Instead, we found that F-1 a,3 (-) can be well approximated by the simple algebraic form
f(0, 200 a(0)= + b(0) (=)](=)
(6)
where a()) and b(0) are specified nonnegative functions of 0 > 1. Intuitively, the linear component shapes the marginal density locally around zero, while the higher-order monomial shapes the tails of the marginal distribution. In Appx. B, we show that the monomial activation f(0. 2) sign(>)|2]" is optimal for Gaussianizing scalar signals whose marginal tail densities satisfy (s()) x |x|t-= exp(-ls|29), where y = 1/0. This closely resembles the tail densities of generalized Gaussian densities (when 1), suggesting monomial activations are effective for shaping the tails.
4.3 Marginal density of local filter responses
We first apply Alg. 1 in the scalar setting N K 1 to demonstrate that choice of activation function in eq. (6) is indeed well matched to the shape of heavy-tailed marginals of local filter responses. For each image, we ran Alg. 1 on the local filter responses with μ 0, lcarning rates (Tgg To) (10-5, 5 10-5) 5 and batch size 10 for 10s iterations. Fig. 2CD shows the learned interneuron activation functions of (0,-) and the learned transforms T(-). Fig. 2E shows histograms of the eircuit responses. Compared to the local filter responses, the circuit responses are visually much closer Gaussian. We found that the circuit performs worse when the distribution ps more 'peaked' around zero (i.e., when is sparser and the fitted shape parameter 3 is smaller), evidenced by the mismatch between the response distribution p, and the Gaussian distribution A N((0, 1) near zero the bottom row of Fig. 2E (see Appx. C for more examples). However, even in this case, the interneuron activation and circuit transform are close to optimal (Fig. 2CD, bottom row).
4.4 Joint density of pairs of local filter responses
Next, we apply Alg. in the multivariate setting N = 2, to pairs of spatially offset filter responses. For each image and offset, we ran Alg. with
interneurons, μ 0, learning rates (Ng, no. The) (10-4, 10-6 10-4), and batch size 10 for 10 iterations. Contour plots of the learned circuit responses for one image are shown in Fig. 3CDE; see Appx. for more examples. The mutual information between circuit responses is shown in Fig. 4. We see that K = 3 interneurons significantly reduces the mutual information between circuit responses for spatial offsets less than d = 32. The reduction is much greater than obtained using K = 2 interneurons and about the same as obtained using = 4 interneurons. For spatial offsets greater than 32, the raw filter responses already have low mutual information and the circuit does not not offer any significant reduction.

---

5 Discussion
We derived a novel online algorithm for transforming a signal to approximate a target distribution, using a recurrent neural circuit with Hebbian synaptic plasticity, gain modulation, and adaptation of interneuron activation functions. The model draws inspiration from the extensive neuroscience literature on efficient coding [15], Hebbian synaptic plasticity [47], interneuron function [29] and gain modulation [48], proposing complementary roles for different physiological processes: Heb- bian synaptic plasticity learns stimulus axes that are least matched to the target distribution, and interneurons adapt their gains and activation functions to transform the marginal responses along these axes. The form of the input-output function for the local interneurons -ainsss--innnnnne nineiarnnnllnear with gain modulation closely resembles phenomenological models of neurons [73], and the parameters W, 0, g can potentially be fit to neural recordings and compared with the optimal parameters that can be derived from the signal statistics Ps. Furthermore, our model predicts a relationship between the interneuron activation function and the circuit transform; see Fig. 2CD for an example of an expansive interneuron activation function that corresponds to a compressive circuit transformation.
There are, however, some aspects of our circuit that are not biologically realistic. For example, our model focuses on the role of local interneurons reshaping the response distribution and, for simplicity, assumes that the primary neurons have linear activation functions. A more realistic model would also include nonlinearity and adaptation in the primary neurons. Moreover, our model only includes synaptic connections between primary neurons and interneurons, which may be consistent with some sensory circuits (c.g., olfactory bulb), but cannot account for excitatory- excitatory connections or inhibitory- cuhibitory-nniibitoryy connections in cortical circuits. Finally, the synaptic weights are not sogttconstaiiee, and violate Dale's law. This can be addressed by modifying the objective in eq. 2 so that the optimization is over non-negative weight matrices W 2 O, which will responsessess result in a projected gradient step in Alg. I; however, the circuit responses will generally be less aligned with the target distribution.
limitation of our simulations that we only test our method on two-dimensional inputs, demon- strating that three interneurons are sufficient to dramatically reduce the redundancy in the circuit responses. However, natural signals are generally high- nigh-dimensional and is not clear how the number of interneurons required to effectively reduce redundancy will scale with the dimension of the signal. There is some basis for optimism. For example, visual inputs are highly structured- e.g. statistical dependencies between inputs rapidly decay with the distance between the inputs- so local interneurons only need to connect to neurons with overlapping or adjacent receptive fields, limiting the number of interneurons that are required the dimensionality of the input signal grows [54].
There are a number of existing computational models that also explain how neural circuits can implement nonlinear transformations to efficiently encode their inputs. For example, several neural circuit models implement forms of divisive normalization [74-76], transformation that optimal for efficient encoding of natural signals [38, 59]. In addition, there is body of work on normative spiking models derived from objectives that maximize the information encoded per spike [77-79], which can account for neural adaptation mechanisms such as gain control. Our work differs from these, by proposing a novel framing of sensory circuit computation in terms of transformations of probability distributions, which can be viewed by Laughlin [ [16]. We then demonstrate in a normative circuit model how interneurons can play a eritical role in optimizing this objective by measuring the marginal distribution of circuit responses and adjusting their feedback accordingly.
population level version of the seminal work
Finally, our results may also be relevant beyond the biological setting. Gaussianization and normal- izing flows are active areas of research [10, 13, 80]. We offer a novel continuous- learning solution inspired by neuroscience that learns using combination of weight updates and activation function up- dates (related to trainable activation functions [81]). In low-dimensional settings, when the constraint functions matched to the signal statistics, we show that Gaussianization can be achieved using relatively few parameters. It is of primary interest to understand how the methods introduced here scale to -igh-dimnnsional signals, where the curse of dimensionality presents significant challenges.
10

---

Acknowledgments
We thank Colin Bredenberg and members of the Center for Computational Neuroscience at the Flatiron Institute for helpful feedback on an earlier draft of this work.
References
[1] Gaspard Monge. Mémoire sur la théorie des déblais et des remblais. Mem. Math. Phys. Acad. Royale Sci., pages 666-704, 1781.
[2] Agnan Kessy, Alex Lewin, and Korbinian Strimmer. Optimal whitening and decorrelation. The American Statistician, 72(4):309-314 2018.
[3] Aapo Hyvärinen and Erkki Oja. Independent component analysis: algorithms and applications. Neural Networks, 1344-5):111 1-340,1140010110 2000.
[4] Mike H Hoyle. Transformations: An introduction and bibliography. International Statistical Review/Revue Internationale de Statistique, pages 203-223, 1973.
[5] Remi M Sakia. The Box-Cox transformation technique: review. Journal of the Royal Statistical Society Series The Statistician, 41(2) 1992.
[6] Scott Chen and Ramesh Gopinath. Gaussianization. Advances Neural Information Processing systems, 13, 2000.
[7] Siwei Lyu and Eero P Simoncelli. Reducing statistical dependencies in natural signals using radial gaussianization. Advances in Neural Information Processing Systems, 1008:1009-1016 2008.
[8] Valero Laparra, Gustavo Camps-Valke Valls, and Jesús Malo. Iterative Gaussianization: from ICA to random rotations. IEEE Transactions on Neural Networks, 22(4):537-549, 2011.
[9] Johannes Ballé, Valero Laparra, and Ecro P Simoncelli. Density modeling of images using a gen- eralized normalization transformation. International Conference on Learning Representations, 2016.
[10] Chenlin Meng, Yang Song, Jiaming Song, and Stefano Ermon. Gaussianization flows. In International Conference on Artificial Intelligence and Statistics, pages 4336-4345. PMLR, 2020.
[11] Danilo Rezende and Shakir Mohamed. Variational inference with normalizing flows. In International conference on machine learning, pages 1530-1538. PMLR, 2015.
[12] Laurent Dinh, David Krueger, and Yoshua Bengio. NICE: Non-linear independent components estimation. International Conference on Learning Representations, 2015.
[13] George Papamakarios, Eric Nalisnick, Danilo Jimenez Rezende, Shakir Mohamed, and Balaji Lakshminarayanan. Normalizing flows for probabilistic modeling and inference. The Journal of Machine Learning Research, (111))261722688 2021.
[14] Fred Attneave. Some informational aspects of visual perception. Psychological Review, 61(3): 183-193, 1954.
[15] H B Barlow. Possible Principles Underlying the Transformations of Sensory Messages. In Sensory Communication, pages 216-234. The MIT Press, 1961.
[16] Simon Laughlin. A simple coding procedure enhances neuron's information capacity. Zeitschrift fur Naturforschung C. Journal of Biosciences, pages 910-912, 1981.
[17] Eero P Simoncelli and Bruno A Olshausen. Natural image statistics and neural representation. Annual Review of Neuroscience, 24(1):1193-1216 2001.
[18] Anthony J Bell and Terrence J Sejnowski, The "independent components" of natural scenes are edge filters. Vision research, 37(23):3327-3338 1997.

---

[19] Bruno A Olshausen and David Field. Emergence of simple-cell receptive field properties by learning sparse code for natural images. Nature, 381(6583):607-609 1996.
[20] Peter Foldiak. Sparse coding in the primate cortex. The Handbook of Brain Theory and Neural Networks, 2003.
[21] Bruno A Olshausen and David Field. Sparse coding of sensory inputs. Current Opinion in Neurobiology, 14(4):481 487,, 2004.
[22] Alexandre Pouget, Sophie Deneve, Jean- Cean-Christophe Ducom, and Peter E Latham. Narrow versus wide tuning curves: What's best for a population code? Neural Computation, 11(1): 85-90, 1999.
[23] Kechen Zhang and Terrence J Sejnowski. Neuronal tuning: To sharpen or broaden? Neural Computation, 11 11(1):75-86 84, 1999.
[24] Deep Ganguli and Eero P Simoncelli. Efficient sensory encoding and Bayesian inference with heterogeneous neural populations. Neural computation, 26(10):2133-2134 2014.
[25] TF Freund and G Buzsáki. Interneurons of the hippocampus. Hippocampus, 6(4):347-470) 1996.
[26] Chris McBain and André Fisahn. Interneurons unbound. Nature Reviews Neuroxcience, 2(1): 11-23, 2001.
[27] Henry Markram, Maria Toledo-Rodriguez, Yun Wang, Anirudh Gupta, Gilad Silberberg, and Caizhi Wu. Interneurons of the neocortical inhibitory system. Nature Reviews Neuroscience, (10):793-88)) 807, 2004.
[28] Carl P Wonders and Stewart Anderson. The origin and specification of cortical interneurons. Nature Reviews Neuroscience, 7(9):687-696, 2006.
[29] Adam Kepecs and Gordon Fishell. Interneuron cell types are fit to function. Nature, 505(7483): 318-326, 2014.
[30] Gord Fishell and Adam Kepecs. Internerron types attractors and controllers. Annual Review of Neuroscience, 43:1-30, 2020.
[31] Mattia Chini, Thomas Pfeffer, and lleana Hanganu-Opatz. An increase of inhibition drives the developmental decorrelation of neural activity. Elife. 11:e78811, 2022.
[32] Cristina Savin, Prashant Joshi, and Jochen Triesch. Independent component analysis in spiking neurons. PLoS Compulational Biology, 6(4):e 1000557, 2010.
[33] Paul D King, Joel Zylberberg, and Michael R De Weeses. Inhibitory interneurons decorrelate excitatory cells to drive sparse code formation in a spiking model of V1. Journal of Neuroscience, 33(13):5475-5485 2013.
[34] Cengiz Pehlevan and Dmitri B Chklovskii. normative theory of adaptive dimensionality reduction in neural networks. Advances in Neural Information Processing Systems, 28, 2015.
[35] Nikolai M Chapochnikov, Cengiz Pehlevan, and Dmitri B Chklovskii. Normative and mecha- nistic model of an adaptive circuit for efficient encoding and feature extraction. Proceedings of the National Academy of Sciences, ( 2023.
[36] Lyndon Duong, Eero Simoncelli, Dmitri Chklovski,, and David Lipshutz. Adaptive whitening with fast gain modulation and slow synaptic plasticity. Advances in Neural Informa- tion Processing systems, 2023.
[37] David J. Heeger. Normalization of cell responses in cat striate cortex. Visual neuroscience, 9(2): 181-97, 1992. ISSN 0952-5238. doi: 10.101777009255230000996
[38] Odelia Schwartz and Fero Simoncelli. Natural signal statistics and sensory gain control. Nature Neuroscience, 489)81198255, August 2001.
12

---

[39] Matteo Carandini and David Heeger. Normalization canonical neural computation. Nature Reviews Neuroscience, 13(1):51-62 62, 2012.
[40] Peter Huber. Projection pursuit. The Annals of Statistics, pages 435-475, 1985.
[41] György Buzsáki and Kenji Mizuseki. The log-dynamic brain: how skewed distributions affect network operations. Nature Reviews Neuroscience, 15(4):264-22) 278, 2014.
[42]
Cédric Villani. Optimal Transport: Old and New, volume 338. Springer, 2009.
[43] Bharath K Sriperumbudur, Kenji Fukumizu, Arthur Gretton, Bernhard Schölkopf, and Gert RG Lanckriet. On integral probability metrics, g-divergences and binary classification. arXiv preprint arXiv:0901.200 2269, 2009.
[44] Harald Cramér and Herman Wold. Some theorems on distribution functions. Journal of the London Mathematical Society, 1(4):290-294, 1936.
[45] Julien Rabin, Gabriel Peyré, Julie Delon, and Mare Bernot. Wasserstein Barycenter and its application texture mixing. In Scale Space and Variational Methods Computer Vision: Third International Conference, SSVM 201 I, Ein-Gedi, Israel, May 29-June 2, 2011, Revised Selected Papers 3, pages 435-446 Springer, 2012.
[46] Ishan Deshpande, Yuan-Ting Hu, Ruoyu Sun, Ayis Pyrros, Nasir Siddiqui, Sanmi Koyejo, Zhizhen Zhao, David Forsyth, and Alexander Schwing. Max-sliced Wasserstein distance and its use for GANs. In Proceedings of the IEEE/CVF Conference Computer Vision and Puttern Recognition, pages 10648-10656, 2019.
[47] Natalia Caporale and Yang Dan. Spike liming-depentie plasticity: Hebbian learning rule. Annu. Rev. Neurosci., 31:25 46,254662226 2008.
[48] Katie A Ferguson and Jessica A Cardin. Mechanisms underlying gain modulation in the cortex. Nature Reviews Neuroscience, 21(2):80-92 020000992222000 2020.
[49] Julijana Gjorgjieva, Guillaume Drion, and Eve Marder. Computational implications of biophys- ical diversity and multiple timescales in neurons and synapses for circuit performance. Current Opinion in Neurobiology, 37:44-52, 2016.
[50] Alison I Weber, Kamesh Krishnamurthy, and Adrienne L Fairhall. Coding principles in adaptation. Annual Review of Vision Science, 5:427 4499,44,440 2019.
[51] DO Hebb. The Organization of Behavior: A Neuropsychological Theory, Wiley & Sons, 1949.
[52] Gina Turrigiano. The self-tuning neuron: synaptic scaling excitatory synapses. Cell, 135 (3):422-435, 2008
[53] David Lipshutz, Cengiz Pehlevan, and Dmitri Chklovskii, Interneurons accelerate learning dynamies in recurrent neural networks for statistical adaptation. In The Eleventh International Conference on Learning Representations, 2023.
[54] Lyndon R Duong, David Lipshutz, David J Heeger. Dmitri B Chklovskii, and Eero P Simoncelli. Adaptive whitening in neural populations with gain-modulating interneurons. Proceedings of the 40th International Conference Machine Learning, PMLR, 202:8902-8921, 2023.
[55] Yang Dan, Joseph J Atick, and R Clay Reid. Efficient coding of natural scenes in the lateral geniculate nucleus: experimental test of computational theory. Journal of Neuroscience, 16 (10):3351-3362, 1996.
[56] Xaq Pitkow and Markus Meister. Decorrelation and efficient coding by retinal ganglion cells. Nature Neuroscience, 15(4):628 6655, 2012.
[57] Michael S Lewicki. Efficient coding of natural sounds. Nature Neuroscience, 5(4):356-364 2002.
13

---

[58] Bernhard Wegmann and Christoph Zetzsche. Statistical dependence between orientation filter outputs used a human- siau--iiion--bsee image code. In Visual Communications and Image Processing 90: Fifth in a Series, volume 1360, pages 909-923. SPIE, 1990.
[59] Siwei Lyu. Dependency reduction with divisive normalization: Justification and effectiveness. Neural Computation, 23(11):2942-2973 2011.
[60] Trevor Hastie, Robert Tibshirani, Jerome H Friedman, and Jerome H Friedman. The Elements of Statistical Learning: Data Mining, Inference, and Prediction, volume 2. Springer, 2009.
[61] Ralph Linsker. Self-organization in perceptual network. Computer, 21(3):005-117, 1988.
[62] Gal Chechik, Amir Globerson, Naftali Tishby, and Yair Weiss. Information bottleneck for Gaussian variables. Advances in Neural Information Processing Systems, 16, 2003.
[63] JH Van Hateren. Processing of natural time series of intensities by the visual system of the blowfly. Vision Research, 37(23)33407-3111 1997.
[64] Sonya Giridhar, Brent Doiron, and Nathaniel N Urban. Timescale-dependent shaping of correlation by olfactory bulb lateral inhibition. Proceedings of the National Academy of Sciences, 108(14):5833-584) 2011.
[65] Andrea Benucci, Aman B Saleem, and Matteo Carandini. Adaptation maintains population homeostasis in primary visual cortex. Nature Neuroscience, 16(6):724-729, 2013.
[66] Adrian A Wanner and Rainer W Friedrich. Whitening of odor representations by the wiring diagram of the olfactory bulb. Nature Neuroscience, 23(3):433 442, 2020.
[67] Ariana R Andrei, Alan Akil, Natasha Kharas, Robert Rosenbaum, Krešimir Josié, and Valentin Dragoi. Rapid compensatory plasticity revealed by dynamic correlated activity in monkeys in vivo. Nature Neuroscience, 2023.
(68] Eero Simoncelli and William T Freeman. The steerable pyramid: A flexible architecture for multi-scale derivative computation. In Proceedings., International Conference on Image Processing, volume 3, pages 444-447. IEEE, 1995.
[69] David Field. Wavelets, vision and the statistics of natural scenes. Philosophical Transactions of the Royal Society of London. Series A: Mathematical, Physical and Engineering Sciences, 357(1760);2277-5442 1999.
[70] Eastman Kodak Company, True color Kodak images dataset, 1993. URL https://wmt math. purdue. aadda/llucier/HOTTTT
[71] Eero P Simoncelli. Statistical models for images: Compression, restoration and synthesis. In Conference Record of the Thirty-First Asilomar Conference on Signals, Systems and Computers, volume 1, pages 673-678. IEEE, 1997.
[72] Siwei Lyu and Eero P Simoncelli. Nonlinear extraction of independent components of natural images using radial gaussianization. Neural Computation, 2186))148--159, 2009.
[73] Robbe LT Goris, J Anthony Movshon, and Eero P Simoncelli, Partitioning neuronal variability, Nature Neuroscience, 17(6):858-88) 2014.
[74] Daniel B Rubin, Stephen D Van Hooser, and Kenneth Miller. The stabilized supralinear network: unifying circuit motif underlying multi-input integration in sensory cortex. Neuron, 85(2):402-417.222 2015.
[75] Matthew Chalk, Paul Masset, Sophie Deneve, and Boris Gutkin. Sensory noise predicts divisive reshaping of receptive fields. PLoS Computational Biology. 1366)e1005582, 2017.
[76] Jesús Malo, José Juan from
and Marcelo Bertalmío. Cortical divisive normalization neural dynamics. Journal of Nonlinear Science, 34(2):1-36 2024.
[77] Veronika Koren and Sophie Denève. Computational account spontaneous activity as a signature of predictive coding. PLoS computational biology, 13(1)xe0005355, 2017.
14

---

[78] Alireza Alemi, Christian Machens, Sophie Deneve, and Jean-Jacques Slotine. Learning nonlin- ear dynamics in efficient, balanced spiking networks using local plasticity rules. Proceedings of the AAAI Conference on Artificial Intelligence, volume 32, 2018.
[79] Gabrielle J Gutierrez and Sophie Denève. Population adaptation in efficient balanced networks. eLife, 8:e46926, September 2019.
[80] Vincent Stimper, David Liu, Andrew Campbell, Vincent Berenz, Lukas Ryll, Bernhard Schölkopf, and José Miguel Hernández-Lobato. normflows: A PyTorch package for nor- malizing flows. Journal of Open Source Software, 8(86);3361, 2023.
[81] Andrea Apicella, Francesco Donnarumma, Francesco Isgrò, and Roberto Prevete. A survey on modern trainable activation functions. Neural Networks, 138:14-32, 2021.
[82] Aapo Hyvärinen and Erkki Oja. Simple neuron models for independent component analysis. International Journal of Neural Systems, 7(06))671-687, 1996.
15

---

A Analysis of circuit transform
In this section, analyze the circuit transform T()) in the that μu 2 0 (i.e., A -1), 94 and 6(0,, -) are convex so that L(W, 0, g, 8, is convex and
T(s) = arg min L(W, 0.g.s.r..
When > O, we can solve for the equilibrium responses using the fixed point iteration: pr (x+1) = W(g f(O, W WTr("))),
pa = =,
where f(0, z) == (f(01, 21).... f(Ok, =x)) denotes the vector obtained by applying the function f elementwise the pairs (0,, 2,) and "o" denotes the elementwise (Hadamard) product.
A.1 Invertibility of the transform
the constraint functions are twice continuously differentiable, then (ss) is continuously differen- tiable and we can relate the response density p. the signal density Ps as follows:
Ps (s) (drt(Jr(s)) pr (T(())
where Jr denotes the Jacobian of T with respect s. In general, we don't have a closed form solution for T' however, setting the update in eq. (4) to zero, see that the inverse transform 7 (when it is well-defined) satisfies
K Σ g.f r.:::...::.w..
where we have used the fact that the partial derivative of c00, z) with respect to = is equal to the partial derivative of h(0, 2) with respect to z. It follows that its Jacobian Jy-1 satisfies
(r) μI Wdiag (gia) .xax (r)) WT,
where (r(r) == (0,, w.)/O=². Provided that either (a) μ > or (b) g: > and (00,,) is strictly convex for all i and the column vectors of W span R., then the Jacobian of 7 is positive definite everywhere. This implies that the Jacobian of T is positive definite everywhere, and thus T invertible everywhere.
A.2 Relation to function class
Equation (7) establishes a relationship between the set of parameterizable (inverse) transforms and the function class (h(0,,)) To better understand this relationship, consider the scalar setting N = = in which case the optimal transformation is given by = maaginal F., where Furget and F, are the
cumulative distribution functions (edfs) Pmnarginal and pa respectively. When T F=1 mmrrgiinl F. follows from eq. (7) that h(0,2) satisfies
@h(0, 2) g0 = = F Fimaggina (z) faz.
g O=
(8)
The multidimensional setting is more complicated; however, in the case that the signal is orthogonal mixture of independent sources, we can derive a precise relationship between the signal density and (W, 0, g) Consider the case that the signal is of the form Au, where A is N x N orthogonal mixing matrix and u = (uj,..., ux) has statistically independent coordinates. Suppose
W AT. After left multiplying eq. (7) on both sides by WT, we get
wW r g f(O, f(0.W
If r is Gaussian, then so is z = W and
1...,N. N. Therefore, an optimal solution eq. (2) is when (0.) and {gg} satisfy (0. F. Fmarginal (2)) - p2iz i == 1, I,
u4 μ2; + 9:f (0,, (0,,
(9)
where Fu, is the cumulative distribution function for u4.
16

---

B Activation functions
In this section, we discuss the activation functions used for the experiments carried out in Sec. 4 of the main text.
B.1 Choice of activation functions
How should we choose the family of activation functions (f} As stated in the main text, one approach is the choose family that is well matched to the marginal statistics of the inputs. In this case, the input marginals are well approximated by generalized Gaussian distributions of the form
(s()) (1/(1 exp(-ls/a/()
with varying scale a and shape 3. Here I is the gamma function. Therefore, from eq. (8), we see that an optimal family of activation functions {S +)3 is such that for each choice of scale and shape B, there a gain and parameter 0 such that
gf =) + pcz F ((2),
where Fa,g and & denote the cdfs of the generalized Gaussian distribution and the standard Gaus- sian distribution. In general, defining the family of activation functions directly in terms of the above display leads to challenges implementing Alg. 1 since 6(0, 2) and Vo0(6, 2) are not readily computable.
We instead sought simple algebraic expression that approximates F-1 ((2). First, suppose the 0,8 activation function takes the form /(0, =) <=gn(=)|=|*. When 0, this activation function is optimal for signal whose edf F momomial satisfies
gsign(=)|=11 F-1 monomial +(2).
Rearranging and differentiating with respect to 2, we find that the pdf of the signal Peomomiia satsifies
1
Pmonomial (gsiinn)()sssinn(())) -|
exp (-2²²
V2Ti
Substituting in with gsign(z)]22 we find the following expression for the pdf
Pmoooonial (s)
a 9V2T 2m
1 exp (-21s//0²)
where q = 1/0. This closely resembles aspects of the pdf for the tails the generalized Gaussian distribution. To reshape the distribution local when ss O, we included linear term, resulting in the activation function
[(0,z) = a(0)= (())sign)))
Here a()) and b(0) are nonnegative functions given by
a(0) exp((20 3.85) 1.95
b(θ) 22. 5.9).
The forms of (0) and b(0) were chosen to approximately minimize ming max, [F, (g/ -0)=)))+(())) when F, is the edf of a generalized Gaussian distribution with shape parameter / between 0.2 and
B.2
Activation function updates
Here we derive the updates to the activation functions.
Recall from Alg. that the 0 update is given
by 0 A + no Vo0(0, 2), where
a(0)
2 (z² -1) 1)+ + b(0) + (|=|0+1 C(0 + 1)),
((0,)) = 110 =) EE-A((..) [0)
17

---

and ((p) is the absolute p-moment of a scalar Gaussian random variable:
2
p p+11
C(p) == No.
π
Note that when 0 = 3, the constraint function normalizes the kurtosis of the marginal and the resulting update bears similarity to ones used in neural models of ICA [82]. Differentiating with respect to 0, we obtain the updates
05(0,2 a' a0)
(2² 1) + (0 (0) b(0)
((++)))
(0 1)2
b(0)
(|=|0+] log = (0 -=|===(++++)))
++1
where
(v)
2D p ++1
(log + 23(00
w,(0)
p ++]
(10)
n
(p) is the polygamma function and we have used the fact that the derivative of the gamma function is (p) = (p))
C Experiments
C.1 Experimental set up
We ran our experiments on cluster comprised of 40-core Intel Skylake nodes with 768GB of RAM. Each experiment shown in Fig. 2 (i.e., Gaussianization the coefficients from image) took less than 5 minutes run. Each experiment shown in Fig. took around 5 hours run (we did not optimize the choice of hyper-parameters)
C.2 Additional experimental results
In Fig. 5, we provide additional examples of histograms of local filter responses and optimized responses for various images from the Kodak dataset [70] (see Fig. the main text for detailed description). We note that when the shape parameter /8 of the fitted generalized Gaussian distribution small, the distribution responses has a characteristic dip near zero. This is likely due to the fact that when is small, there is more probability mass concentrated near zero and so small discrepancies between the learned activation function ( 2) and the optimal activation function F F (2(2 when x 0 can lead to large discrepancies between the response distribution and \ 1) near zero. This could potentially be resolved by choosing a parameterization gf 2) that better approximates optimal activation function.
In Fig. 6, we provide additional examples of contour plots local filter responses and optimized responses for three additional images from the Kodak dataset (corresponds to top three rows of Fig. 5, see Fig. 3 of the main text for detailed description). In some of these examples we see that unlike the standard Gaussian distribution (O, I) the learned response distribution is clearly non-monotone with respect to [|]]]- This non-monotonicity likely inherited from non-optimality of the interneuron activation functions that leads to the discrepancies in the scalar case (see Fig. 5).
18

---

mage
218.0
30-1
-2
30"
-2
10-
30
0.05.
10
30-2
2-0000 0
10-1
10-7
10-1
0-0. 0.55
2 0.00, a
10
30"1
10-2
- 0.000 d 2.55
10+1
30-
282
Figure 5: Gaussianization of local filter responses.
19

---

A) signal
B) ZCA
C) responses (K=2)
D) responses (K=3)
E) responses (K=4)
Hi-2 36e
Mi-) 15e1
Ni-2e2
Mi-1 76e-2
MI-9.01e-
Mi-le-2
2002
25
MI-1 1i1e2
1.202
2=1,50022
2
A) signal
B) ZCA
C) responses (K-2)
D) responses (K- 31
E) responses (K-4)
Mi-2.39e-1
Nim2-4te-1
Ni-100t-11
NI-402022
MI-34tb22
Mi-1000-1
Niwg Niw2.196-2
Niw /.Bl
24
Mi-d.50 65c 2
Mind 2
35e
26-2.762 2
M0-25002 2
signal
B) ZCA
responses (K-2)
responses (K-3) 31
(K-4)
M=1.01e-1
Niw1.1le1
Nie4Bie2
NI-194e2
Mi=1:75e22
Nim 2
Nle6.5le2
Ml-14l22
Ml-237e2
Ni-1.000
Mir M
Ni-n.Ate-1
MI-8 00-0-bb--1 40m1
23
Figure 6: Gaussianization of pairs of local filter responses.
20

---

NeurIPS Paper Checklist
1. Claims
Question: Do the main claims made in the abstract and introduction accurately reflect the paper's contributions and scope?
Answer: [Yes]
Justification: Our claims in the abstract and introduction match the theoretical and experi- mental results in the paper.
2. Limitations
Question:
Does the paper discuss the limitations of the work performed by the authors?
Answer: [Yes]
Justification: We include limitations of our work and potential alternative approaches in the Discussion section.
3. Theory Assumptions and Proofs
Question: For each theoretical result, does the paper provide the full set of assumptions and complete (and correct) proof?
Answer: [Yes]
Justification: While our paper does not include proofs, we do include full details of our analytical results.
4. Experimental Result Reproducibility
Does the paper fully disclose all the information needed to reproduce the main ex- perimental results the paper to the extent that affects the main claims and/or conclusions of the paper (regardless of whether the code and data are provided or not)?
Question:
Answer: [Yes]
Justification: We describe the implementation of our algorithm and include link to code for reproducing our results.
5. Open access to data and code
Question: Does the paper provide open access to the data and code, with sufficient instruc- tions faithfully reproduce the main experimental results, described in supplemental material?
Answer: [Yes]
Justification: We link GitHub repository with example code for reproducing the experi- ments described in the paper.
6. Experimental Setting/ettails
Question: Does the paper specify all the training and test details (e.g., data splits, hyper- parameters, how they were chosen, type of optimizer, etc.) necessary to understand the results?
Answer: [Yes]
Justification: We provide full details of our experiments.
7. 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 include error bars.
8. Experiments Compute Resources
Question: For each experiment, does the paper provide sufficient information the puter resources (type of compute workers, memory, time of execution) needed to reproduce the experiments?
Answer: [Yes]
21

---

Justification: We include the amount of compute used and runtime for the experiments (see Appx. C).
9. Code Of Ethies
Question: Does the research conducted in the paper conform, in every respect, with the NeurIPS Code of Ethics https: : http///eurrppp.ppe cccpublic/Ethice.cooocccccccccce http://tehhisiieeeniiniiiieieeee
Answer: [Yes]
Justification: Our research conforms with the NeurlPS Code of Ethics.
10. Broader Impacts
Question: Does the paper discuss both potential positive societal impacts and negative societal impacts of the work performed?
Answer: [NA]
Justification:
We do not foresee an immediate societal impact of our work.
11. Safeguards
Does the paper describe safeguards that have been put in place for responsible release of data models that have high risk for misuse (e.g., pretrained language models, image generators, or scraped datasets)?
Question:
Answer: [NA]
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 cite the Kodak dataset and include a URL to the asset.
13. New Assets
Question: Are new assets introduced in the paper well documented and is the documentation provided alongside the assets?
Answer: [NA]
14. 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 details about compensation (if any)?
Answer: [NA]
15. Institutional Review Board (IRB) Approvals or Equivalent for Research with Human Subjects
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?
Question:
Answer: [NA]
22


Reviews:
[{"review_id": "XWbwj19qcd", "reviewer": "Reviewer_zp8S", "summary": "This paper introduces a method that utilizes multiple inhibitory neurons and the Hebbian learning rule to transform a signal into a representation that conforms to a specified target distribution. Specifically, the model incorporates Hebbian synaptic plasticity to establish connections that optimally match this target distribution. Interneurons within the model are adapted in terms of their gain and activation directions to respond most effectively to the input signals, mirroring the potential optimization mechanisms found in biological neural processing.\n\nUtilizing this approach, the study processes natural images from the Kodak dataset through a wavelet transform. This extracts 2-dimensional signals, specifically pairs of wavelet coefficients for images at fixed horizontal spatial offsets ranging from d = 2 to d = 64. The model then transforms the distribution of this two-dimensional information into a Gaussian distribution. This transformation demonstrates the model's ability to handle complex data structures and align them with statistically predictable patterns, enhancing both the analysis and interpretation of natural images.", "questions": "1. As mentioned in the Weaknesses section, the model is not biologically realistic. Given this, my primary question is: What is the purpose of using neural networks to implement this function? Furthermore, if it does not realistically mimic biological processes, should it be compared with more engineering-oriented approaches?\n\n2. Scaling of Inhibitory Neurons with Input Dimensionality: The model currently demonstrates with two-dimensional inputs, requiring three inhibitory neurons. How does the number of required inhibitory neurons scale as the dimensionality of the input increases? This is crucial for understanding the feasibility and complexity of the model when applied to higher-dimensional data.", "limitations": "Yes, the authors have addressed some of the limitations of their work in the Discussion section of the paper. They specifically mention the model’s lack of biological realism and its primary demonstrations within low-dimensional settings. These acknowledgments align with the checklist guidelines on discussing limitations.", "rating": 5, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 1, "strengths": "1. Clear Presentation: The paper articulates its concepts with exceptional clarity, facilitating a deep understanding of complex models and their applications.\n\n2. Effective Transformation of Distributions: The model is highly effective at converting input signals into specific, targeted distributions, optimizing data for further processing and analysis.\n\n3. Online Operational Capability: One of the model's significant advantages is its ability to operate online. This feature significantly enhances its practicality for real-world applications, allowing for real-time data processing and continuous learning without the need for retraining.", "weaknesses": "1.  A main concern is that the model is not biologically realistic. In addition to the authors' admission that the feedforward synaptic weights $W^T$ and feedback weights $-W$ are symmetrically constrained, there is also a question about the biological plausibility of using pairs of wavelet coefficients as inputs. It's unclear how biological systems would naturally derive such information and whether the distribution of these inputs is something that biological systems need to transform.\n\n2. The setup of online learning with a batch size of 10 also raises questions about its biological feasibility. It's uncertain how biological systems could implement a similar mechanism.\n\n3. The paper only demonstrates the method's effectiveness on two-dimensional inputs, which might not sufficiently prove its efficacy. In real-world scenarios, we often deal with inputs that are high-dimensional ($N>>2$), and the model's performance in such conditions remains untested.\n\n4. The paper incorporates many specific settings, such as the choice of activation functions based on the Gaussianization of scalar signals. These particular choices may limit the generalizability and applicability of the model to different datasets or broader applications.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for your careful reading of our paper and for your thoughtful comments. We are pleased that you find the paper clear and appreciate the models online capabilities! We understand your concerns about the biological realism of our model and have addressed these concerns in our general author rebuttal above. Although there are aspects of our model that are not matched to underlying biological details, we do believe it offers a novel perspective for thinking about circuit computation that can provide a framework for interpreting experimental measurements.\n\nIn response to weakness #4: We actually think this points to an interesting aspect of sensory signals. In some respects, it's quite surprising that the model can effectively Gaussianize two-dimensional signals just by reshaping the response marginals along a few axes. This is likely due to structure of natural signals, and the efficient coding hypothesis posits that sensory circuits take advantage of this as well. Therefore, while we agree there are aspects of our circuit that may not generalize to arbitrary datasets, we do not think this reduces its efficacy as a model of circuit computation.", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "IWB7OY9xKd", "reviewer": "Reviewer_Voaj", "summary": "Authors propose an online algorithm that solves the problem of optimal transport. In particular, assuming a spherical distribution of stimuli, the goal of the algorithm is to generate neural responses such that their distribution best approximates the distribution of stimuli. Authors find that a neural network with excitatory and inhibitory neurons can solve such optimization problem, by learning a non-linear transformation that reduces dependency between neural responses. This is done through Hebbian learning on synapses and by adjusting activation functions of single neurons.", "questions": "Authors motivate their modelling by claiming that it is currently unclear how neural networks could implement nonlinear transformations. Which non-linear transformations are not captured by the listed models?\nThe paper by Alemi et al. AAAI 2018, proposes efficient spiking models that implement non-linear transformation of stimuli at the population level. In Koren & Panzeri, NeurIPS 2022, neural populations perform linear transformation of stimuli on the level of neural populations, but on the level of single neurons these transformations are strongly non-linear. What do authors mean by non-linear transformations? Non-linear on what level? This should be clarified.\n\nAuthors say that approximating target distribution is useful because it facilitates efficient transmission of signals. It is however not entirely clear why making the distribution of responses Gaussian-like facilitates efficient transmission. Is this improvement in efficient transmission attributed to redundancy reduction?\n\nIs inhibition required to solve this optimization problem, e.g., could the same operation be achieved through plasticity of E neurons?", "limitations": "Authors addressed the limitations adequately.", "rating": 5, "confidence": 3, "soundness": 3, "presentation": 2, "contribution": 2, "strengths": "The paper presents is a novel combination of well known techniques and provides an alternative to existing approaches. The work seems technically sound and authors describe some interesting relation to similar approaches. Proposed algorithm is a nonlinear extension of existing algorithms for data whitening using a neural network.", "weaknesses": "The paper has the ambition of finding a neural implementation of generalized whitening such that it could be also implemented by biological neural circuits. However, a major concern is that inhibitory neurons do not obey Dale's law, and the model is therefore not biologically plausible. In general, it is not clear what the paper brings to understanding of biological or artificial neural networks.\nSince the model violates Dale's law, the discussion about the function of interneurons seems out of scope.\n\nMoreover, there are significant parts of the paper that are unclear. A number of details about the methods are given, but the main mechanism that supports learning remains unclear. Authors mention that the circuit \"learns directions that are least matched to the target distribution\", but it remained unclear to me how is this achieved and why this helps to Gaussianize the distribution of output firing rates of principal neurons. \n\nAuthors motivate their model by claiming that divisive normalization has not yet been solved with a neural network. This seems false (see for example Chalk, Masset, Deneve & Gutkin,  PLOS CB 2017). Gain control is also captured by available models. Efficient spiking models and normative models where changes in gain control can be captured by modulating the metabolic cost on spiking in the model's objective (see Guiterrez and Deneve, eLife 2019; Koren and Deneve, PLOS CB 2017). In general, the paper does not sufficiently take into account previous work on efficient coding with spikes that is closely related.  \n\nThere are typos on several places, for example: lines 75, 77, 100, 234", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for your careful reading of our paper. We take your concerns quite seriously and have revised our paper accordingly. Please find our responses to your listed weakness and your questions below.\n\n### Weaknesses\n- We tested a modified version of our algorithm in which we enforce Dale's law. The modified version performs well, although there is a noticeable degradation in performance when compared with the original algorithm. For more details, please see our general author response and the attached PDF.\n- We appreciate your feedback that parts of the paper are unclear. We'll do our best to improve the revised description.\n- Thank you for pointing out these omissions. We take this quite seriously and plan to revise our paper to make it clear that there are a number of existing nonlinear neural circuits for efficiently encoding their inputs (see our general author response).\n- Thank you for pointing out these typos.\n\n### Questions\n- To be precise, we mean that the circuit transform is nonlinear. Specifically, the function $T:{\\bf s}\\mapsto{\\bf r}$, which maps the vector of circuit inputs to the vector of circuit responses, is nonlinear. We will clarify this in our revision.\n- Yes, the primary reason is attributed to redundancy reduction. There are other factors that make Gaussian distributions appealing, but this is the main reason.\n- It is likely that plasticity of excitatory neurons can contribute to reshaping of neural responses; however, redundancy reduction likely requires local inhibitory interneurons. For example, we are unaware how plasticity of EE connections can reduce correlations between two primary neuron responses that receive highly correlated inputs.", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "ef7tcif7Hi", "reviewer": "Reviewer_txmi", "summary": "This paper investigates a crucial question in neuroscience: how do neural circuits convert inputs into a target distribution, specifically focusing on how local interneurons transform natural image statistics into a spherical Gaussian distribution. The authors approach this problem through the lens of optimal transport. By using an integral probability metric, the task of transforming the input distribution into a spherical Gaussian distribution is framed as a minimax optimization problem. This constrained minimax optimization translates to recurrent neural dynamics (representing the minimization part) combined with gain modulation, activation function adaptation, and plasticity (representing the maximization part), with plasticity being local and Hebbian. Notably, the study's principled adjustment of activation functions is a unique contribution. The authors present experimental results on natural images using wavelet transformations.", "questions": "The constraint $g_i > 0$ is not included in Algorithm 1. Would imposing these constraints ($w > 0$, $g > 0$) negatively impact the experimental results?", "limitations": "The authors have adequately addressed the limitation.", "rating": 5, "confidence": 2, "soundness": 3, "presentation": 2, "contribution": 2, "strengths": "Note: As I am not well-versed in wavelet transformations, which are extensively utilized in the experimental section, my review focuses primarily on the theoretical aspects prior to the experiments (Section 4).\n\n**Originality**: This paper tackles a significant question in neuroscience: how neural circuits transform input signals into a target distribution. The authors frame this problem as an optimal transport problem. Utilizing the integral probability metric, they demonstrate that this problem is equivalent to a constrained minimax optimization, which can be interpreted through neural circuit dynamics. The paper's theoretical contributions are both solid and original.\n\n**Quality**: The work is theoretically robust, with well-justified methods and experiments that effectively validate the theoretical model. The use of spherical Gaussian distributions, which might initially appear to be a limitation, is convincingly justified in Section 4.\n\n**Clarity**: The paper is clearly written overall.\n\n**Significance**: This paper introduces a normative model for normalizing inputs, offering potentially valuable insights for the theoretical neuroscience community.", "weaknesses": "- The authors should define $p_{target}$, $p_{marginal}$, and $p_r$ explicitly instead of using descriptive terms for better clarity.\n- The proposed model necessitates different $g$, $\\theta$, and $w$ for different stimuli. While neurons can potentially adjust their gain and activation function (i.e., $g$, $\\theta$) relatively quickly, the need for weight updates to handle varying inputs significantly limits the model's applicability and biological plausibility.\n- As noted in the limitations section, the weights are not sign-constrained. Therefore the distinguishing between excitatory and inhibitory neurons in this model makes no sense.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for your comments. We will revise our paper in accordance with your suggestions.\n\n### Weaknesses\n- We will explicitly define these distributions in our revision.\n- This is a great point. Interestingly, the learned weights $W$ are approximately shared across images, whereas the optimal $g$ and $\\theta$ vary between images, suggesting that there are structural properties that are shared across images. This is consistent with rapid adjustments of $g$ and $\\theta$ and slow adjustments of $W$. We did focus on this point because Duong et al. [NeurIPS 2023] have shown similar results. Specifically, they showed how a multi-timescale (linear) circuit with fast gain adjustments and slow synaptic updates can effectively whiten circuit responses to changing contexts. Here, our focus is on the nonlinear aspects of the circuit so we did not emphasize this point, but we will note it in our revision.\n- We can enforce $W>0$ at the cost of a slight degradation in the performance of the algorithm; see our general author response above and the attached PDF.\n\n### Questions\nIndeed, we did not include the constraint $g_i>0$ in Algorithm 1. However, the optimal solutions all have $g_i>0$ and including the constraint would not have changed the experimental results. As mentioned above, enforcing $W>0$ results in a slight degradation in performance.", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "kd8TsD8WFa", "reviewer": "Reviewer_DTUH", "summary": "This paper proposes a normative recurrent circuit to solve an optimal transport problem, focusing on the problem of Gaussianization of natural image statistics.", "questions": "- Along with the work of Van Hateren, the authors might consider citing Laughlin, \"A simple coding procedure enhances a neuron's information capacity,\" and Juusola and Hardie, \"Light Adaptation in Drosophila Photoreceptors.\" I leave it to their discretion whether or not to do so.\n\n- In that vein, it might be nice to provide further evidence for adaptation of activation functions in biology. If additional references occur to me, I will update my review. \n\n- In Line 271, should \"activations\" be plural? I think this is a typo. \n\n- The legibility of Figure 3A could be improved; as it stands the text is a bit too small and the lines a bit too faint. \n\n- It could be nice to provide further support for the approximate form of the activation functions. Can you numerically evaluate the expression below Line 630?", "limitations": "As noted above, the authors clearly discuss the limitations of their work.", "rating": 8, "confidence": 4, "soundness": 3, "presentation": 4, "contribution": 4, "strengths": "This paper was a pleasure to read. The writing is clear, the formulation of the problem elegant, and the results represent a clear advance relative to past works on whitening circuits.", "weaknesses": "I have only a few small concerns regarding the biological realism of the model, but these are mostly already acknowledged by the authors. One point that I believe warrants further discussion is the plausibility of the activation functions (6). These are somewhat reminiscent of two-sided versions of the rectified-power law activations used in cortical models, no? I also have a handful of miscellaneous questions (see below), but on the whole these do not dampen my enthusiasm.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for your review, we are pleased that you enjoyed reading our work!\n\n### Weaknesses\n\nWe agree the plausibility of the activation functions warrants further discussion. The interneuron activation functions indeed resemble two-sided rectified-power law activations. A potentially more plausible variant would have interneurons responding with a  *rectified* power law, although implementing the same computation would then require twice as many interneurons.\n\nAlternatively, our model provides a precise relation between the input distribution, target distribution and interneuron activation function. Therefore, fixing two of these determines the third. We fixed the input distribution and target distribution and derived the optimal interneuron activation function. One could instead fix the input distribution and interneuron activation function (e.g., based on experimental measurements) and estimate the target distribution. \n\n### Questions\n\n- Thank you for the suggestion. We will include references to the works by Laughlin and Juusola & Hardie in our revision. \n- We will also add a reference to Sharpee et al. \"Adaptive filtering enhances information transmission in visual cortex\", which showed that neurons in cat primary visual cortex adapt to higher-order statistics beyond mean and covariance.\n- Thanks for catching this typo.\n- We will improve the legibility of Figure 3A. \n- Yes, numerical evaluations of the expression below Line 630 are shown as thick gray curves in Figure 2C. To emphasize this, we will add text below Line 630 pointing the reader to the figure.", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}]

Author rebuttal:
Thank you for your careful reading and helpful comments. Here we respond to two important concerns and provide individual responses below.

## Biological realism

Reviewers **zp8S**, **Voaj** and **txmi** listed the biological realism of our model as a primary concern that limits the applicability of our model in understanding neural circuits. We appreciate your concerns, though we still believe that our model represents a useful advance.

First, some of the concerns can be addressed with minor adjustments that do not affect the overall circuit computation. We list these concerns with the most easily addressed concerns first.
1. **Wavelet coefficients.** Apologies for the jargon.  These are responses of a local oriented filters (similar to Gabor filters) applied to natural images. They are qualitatively similar to simple cell responses in primary visual cortex [Field 1987], so they are representative of natural inputs to cortical circuits in the visual cortex. We will edit the text and refer to the inputs as "local filter responses''.
2. **Batch training.** Using a batch size 10 was an optimization choice and is not required. We could have implemented a fully online optimization algorithm with batch size 1. To demonstrate this, we've run Algorithm 1 in the fully online setting and reproduced Figure 3B of the main text; see the attached PDF. 
3. **Symmetry of synaptic weights.** While Algorithm 1 enforces symmetry between the primary neuron-to-interneuron weights $W^\top$ and interneuron-to-primary neuron weights $-W$, this is not required. When the weights are decoupled, they converge asymptotically toward a symmetric solution due to the symmetry of their learning rules. This has been demonstrated both analytically and numerically in previous works; see, e.g, appendix E.2 of [Duong et al., NeurIPS 2023].
4. **Scalability.** We have strong reason to believe that the number of inhibitory interneurons required will scale reasonably with the dimension of the input. Visual inputs are highly structured; e.g., statistical dependencies between inputs rapidly decay with the distance between the inputs. Therefore, local interneurons only need to connect to neurons with overlapping or adjacent receptive fields, which greatly reduces the number of interneurons that are required as the size of the model increases.
5. **Violation of Dale's law.** This is the most serious concern (although note that Dale's law may not be as hard of a constraint as previously believed; see [Saunders et al. 2016, Granger et al. 2020]). There are two possible solutions. One approach is to seek a circuit implementation of the algorithm that obeys Dale's law. This could potentially be achieved by introducing additional (excitatory or inhibitory) neural populations, though the details of this construction need to be worked out. An alternative approach is to replace the current optimization problem with a *constrained* optimization problem that respects Dale's law and can be optimized via a projected gradient descent algorithm. As a preliminary test of this, we optimized the constrained optimization problem (in the context of Gaussianization of visual inputs) and find that the solution still significantly reduces the mutual information between neural responses (more than linear whitening, but less than the unconstrained model that violates Dale's law); see the attached PDF.

Second, we believe that our normative model strikes an effective balance between simplicity and biological realism that provides a useful framework for interpreting experimental measurements of neural circuits. In particular:
1. Our model conceptualizes the circuit computation as transforming the input distribution to achieve a target distribution (or set of distributions). This can potentially be tested experimentally by measuring distributions of circuit responses (and how these change when the input distribution shifts). We think this is especially useful when recording from large populations of neurons and there is some experimental work along these lines [Benucci et al. 2013].
2. We establish a link between the circuit level computation and physiological processes such as neural adaptation in local interneurons. This can be experimentally tested by measuring how interneuron FI curves adapt in response to changes in the distribution of circuit inputs.

## Relation to existing work

Reviewer **Voaj** pointed out that our original submission does not cite salient works related to circuit models that implement nonlinear transformations to efficiently encode information. We take this point quite seriously and will work to improve our description of the relationship of our work to related literature.  The following will be added to the paper:
> There are a number of existing computational models that explain how neural circuits can implement nonlinear transformations to efficiently encode their inputs. For example, there are several neural circuit models that implement forms of divisive normalization [Rubin et al. 2015, Chalk et al. 2017, Malo et al. 2024], a transformation that is optimal for efficient encoding of natural signals [Schwartz & Simoncelli 2001, Lyu 2011]. In addition, there is a body of work on normative spiking models derived from objectives which maximize the information encoded per spike [Koren & Deneve 2017, Alemi et al. 2018, Gutierrez & Deneve 2019], which can account for neural adaptation mechanisms such as gain control. Our work differs from these, by proposing a novel framing of sensory circuit computation in terms of transformations of probability distributions, which can be viewed as a population level version of the seminal work by Laughlin. We then demonstrate in a normative circuit model how interneurons can play a critical role in optimizing this objective by measuring the marginal distribution of circuit responses and adjusting their feedback accordingly.

Meta-review: The present study starts from the optimal transport objective function, and derives recurrent circuit dynamics with nonlinear interneurons to optimize the objective function. It links the normative computation and canonical cortical microcircuit, and technically it is a nonlinear extension of a recent circuit model for linear data whitening (Duong, et al., ICML 23; NeurIPS 23). Although the manuscript distinguishes the present study with Duong's papers, it would be better to discuss how much we gain by considering nonlinear interneurons from the perspective of computational performance and neuroscience. Lastly, please revise the manuscript accordingly based on the reviewers' feedback, including the discussions of Dale's law, symmetry of connection weight, scalability, etc.
Decision: Accept (poster)