Shaping the distribution of neural responses with interneurons in a recurrent circuit model

Authors: david  lipshutz
eero p simoncelli
Venue: NeurIPS_cc_2024_Conference
Published: 2024-05-13 12:52:36

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.

Main paper: 
 introduction: 
 The problem of transforming a signal into a representation with a given target distribution (or within a target set of distributions) is a classical problem whose origins can be traced back more than two centuries [1] . Many methods in statistics, signal processing and machine learning can be interpreted within the context of this problem. 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 a signal processing method that linearly transforms a signal so as to minimize higher-order statistical dependencies in addition to removing second-order dependencies. Nonlinear transformations of skewed or heavy-tailed data to approximately Gaussianize their distribution can facilitate statistical analyses [4, 5] . Machine learning methods for density estimation such as Gaussianization [6] [7] [8] [9] [10] and normalizing flows [11] [12] [13] nonlinearly transform high-dimensional signals with complex densities into more tractable representations with approximately Gaussian densities. This problem may also lie at the core of sensory processing. Efficient coding theory posits that sensory systems maximize the information they transmit about sensory signals to downstream areas subject to resource constraints [14] [15] [16] [17] . In one instantiation of this theory, the redundancy reduction hypothesis posits that sensory circuits transform natural signals into representations to minimize or eliminate statistical dependencies between coordinates, essentially producing factorized response distributions [14, 15, 18] . In a separate, but related, instantiation, sparse coding theory posits that population responses are optimized for sparsity [19] [20] [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 to maximize the Fisher information about the inputs [22] [23] [24] , which can be interpreted as a 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, it is 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 information flow [25] [26] [27] [28] [29] [30] [31] . Consequently, local interneurons are natural candidate substrates responsible for shaping circuit responses into efficient representations. However, the precise relationship between the physiological and anatomical properties of local interneurons and the coordinated response properties of populations of primary neurons remains unclear.Several normative mechanistic models have been proposed to explain how neural circuits can linearly transform their inputs into a representation whose distribution lies within a target set (e.g., the set of distributions with identity covariance) [32] [33] [34] [35] . These models are derived from optimization objectives for linear redundancy reduction, including (adaptive) decorrelation and ICA, and the circuit parameters (e.g., gains, synaptic weights) are optimized to match the data distribution. 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 a 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 [36] [37] [38] that are not captured by these models.Here, we seek a normative model of how circuits nonlinearly transform their inputs to produce responses with a (spherical) target distribution. Starting from an optimal transport objective for transforming the input signal into a neural representation with a given target distribution, we 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, projection pursuit [39] ); (ii) interneuron gains and nonlinear activation functions 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 a primary test case motivated by redundancy reduction theory, we apply our algorithm to the case that the inputs are derived from natural images (using local oriented filters qualitatively similar to those found in the primary visual cortex) and the target distribution is the spherical Gaussian. 1 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. 
circuit objective: 
 Consider a neural circuit with N ≥ 2 primary neurons that transforms input signals s ∈ R N , which are distributed according to a density p s , into circuit responses r ∈ R N (Fig. 1 ). The inputs may represent a direct sensory input (e.g., the rate at which photons are absorbed by a cone) or the weighted sum of multiple inputs (e.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 : R N → R N .We explore the possibility that the circuit objective is to transform its inputs s so that the circuit responses r = T (s) follow a (spherical) target distribution p target while minimizing the L 2 -distance between responses and input signals. Mathematically, this corresponds to an optimal transport problem [41] minEQUATIONwhere the minimization is over a suitable class of functions T , λ ∈ R is a regularizing term and, unless otherwise noted, expectations are over the input distribution p s . Note that the choice of λ ∈ R does not affect the optimal solution as E[∥T (s)∥ 2 ] is fixed provided T (s) ∼ p target ; however, it will affect the optimization algorithm. Assuming p s is sufficiently regular, the minimum in eq. ( 1) is the squared Wasserstein-2 distance between p s and p target and the input-response transformation T is the so-called optimal transport plan that maps the input distribution p s to the target distribution p target .Our goal is to derive an online algorithm for optimizing the objective in eq. ( 1) that can be implemented in a neural circuit model. We accomplish this by defining a distance between the response distribution p r and the target distribution p target that can be estimated in a neural circuit, and then solving the optimization problem by incorporating this measure as a constraint, using the method of Lagrange multipliers.Common candidates for measuring the discrepancy between two distributions include Kullback-Leibler (KL) divergence or integral probability metrics [42] . The former typically require numerous samples to estimate the density p r , 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 p target , denoted p marginal , which are all equal under our assumption that p target is spherical. Such constraint functions 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 a linear mixture model, the Cramér and Wold theorem [43] 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 as 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 ≥ 1 unit vectors w 1 , . . . , w K ∈ R N , which can be randomly sampled, chosen based on prior knowledge of the signal statistics, or learned from data using projection pursuit [39] . We then choose a class of scalar functions {h(θ, •)} parameterized by θ, which defines a semi-metric between the marginal of r in the direction w and p marginal to bemax θ |E[ϕ(θ, r • w)]|, where ϕ(θ, z) := h(θ, z) -E z∼pmarginal [h(θ, z)].For example, when {h(θ, •)} parameterizes all Lipschitz-1 (indicator) functions, this induces the Wasserstein-1 (total variation) distance between the marginal response distribution and p marginal . Given directions w 1 , . . . , w K and constraint functions {h(θ, •)}, we express the distance between the response density p r and the standard Gaussian distribution as the sum of the marginal distances:d W,ϕ (p r ) := K i=1 max θi |E r∼pr [ϕ(θ i , r • w i )] |,whereW := [w 1 , . . . , w K ] is the N × K matrix of concatenated unit vectors. 2How do we choose the constraint functions {h(θ, •)}? In general, the choice should be well-suited to the input distribution p s and the target distribution p r . For example, consider the simple case when the inputs follow a centered Gaussian distribution with unknown covariance structure, and the target distribution is the spherical Gaussian distribution N (0, I). The Wasserstein-2 distance between a marginal distribution of the inputs and the standard normal distribution N (0, 1) can be expressed in terms of the difference between the second moments, so a quadratic function h(z) = 1 2 z2 suffices. In section 4.2, we consider a parametric class motivated by natural signal statistics.We replace the condition r ∼ p target in eq. ( 1) with the condition max W d W,ϕ (p r ) = 0, which we enforce using Lagrange multipliers. This results in the minimax optimization problemEQUATIONwhere L is defined byEQUATIONand θ := (θ 1 , . . . , θ K ) is the set of concatenated parameters, g := (g 1 , . . . , g K ) is a Kdimensional vector of Lagrange multipliers, and the circuit transform is defined by T (s) = arg min r L(W, θ, g, s, r). The maximization over g and θ effectively minimizes the distances between the marginal distributions of the responses r along the directions w 1 , . . . , w K and p marginal , whereas the maximization over the matrix W learns directions along which the marginals of s are least aligned with p marginal , essentially performing projection pursuit [39] . 
algorithm and circuit implementation: 
 We now derive an online gradient-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 or synaptic adaptation mechanisms, so we assume that the neural activities equilibrate before the neural activations and synapses are updated.For each iteration, the circuit receives a stimulus s. The (discretized) recurrent neural response dynamics (Fig. 1 ) correspond to gradient-descent minimization of L with respect to r:EQUATIONwhere η r > 0 is a small constant, µ := 1 + λ represents a leak term,z i := r • w i is the weighted input to the i th interneuron, f (θ i , •) := ∂ϕ(θ i ,•)/∂z is the activation function, g i is a multiplicative gain that scales the output, and n i denotes the output. Notably, the interneuron activation response function g i f (θ i , •) is parameterized by (g i , θ i ), which can vary across interneurons, so the interneuron responses are heterogeneous. For each i, synaptic weights w i connect the primary neurons to the i th interneuron and symmetric weights -w i connect the i th interneuron 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 n := (n 1 , . . . , n K ).Since the neural activities equilibrate before other updates are performed, the responses r are a fixed point of L(W, θ, g, s, •). In general, we do not have a closed-form expression for r; however, if g i ≥ 0 and ϕ(θ i , •) are convex, then L(W, θ, g, s, •) is convex and we can express the transform asEQUATIONIn Appx. A, we show that the transform T (•) is invertible whenever g i ≥ 0 and ϕ(θ i , •) are convex and it defines a precise relationship between the input distribution p s and response distribution p r .After the neural activities reach equilibrium, we maximize L by taking concurrent gradient-ascent steps with respect to g i , θ i and w i :∆g i = η g ϕ(θ i , z i ), ∆θ i = η θ ∇ θ ϕ(θ i , z i ), ∆w i = η w n i r, where η g , η θ , η w ≥ 0 are the respective learning rates, which control the relative speeds of gain modulation, neural adaptation and synaptic plasticity, respectively. For example, at the extremes, we can fix the gains, activation functions or synaptic weights by setting η g = 0, η θ = 0 or η w = 0, respectively. Notably, while synaptic plasticity [46] and gain modulation [47] are well-studied circuit mechanisms that support learning and adaptation, adjustments of nonlinear neural activation functions are not as well established though there is some emerging evidence that neurons also adapt their activation functions in response to changes in their input statistics [48, 49] .The circuit operates online and the updates are local in the sense that the updates to the gain and activation function of the i th interneuron only depend on variables θ i and z i . The updates to the synapses w i (-w i ) are proportional (inversely proportional) to the product of the pre-and postsynaptic activities, so they are both local and Hebbian (anti-Hebbian) [50] . Finally, to ensure that the vectors w 1 , . . . , w K have unit norm, we normalize the weights after each update: w i ← w i /∥w i ∥. This can be viewed as form of homeostatic plasticity such as synaptic scaling [51] . While the feedforward weights w i and feedback weights -w i are constrained to be symmetric, these can be decoupled due to the symmetry of the Hebbian learning rule. Both theoretical and empirical evidence of this is shown for a related adaptive whitening circuit in [52, appendix E.2].Combining the neural dynamics, the interneuron adaptation and synaptic plasticity steps yields our online algorithm (Alg. 1), which we write in vector-matrix notation by defining the normalization function W ← P (W + η w r t n ⊤ t ) ; // Hebbian + homeostatic plasticity 13: end forP (W) := [w 1 /∥w 1 ∥, . . . , w K /∥w K ∥].Algorithm 1 is naturally viewed as a nonlinear extension of existing algorithms for linear data whitening that have neural circuit implementations [34, [52] [53] [54] . In particular, when the constraint function is quadratic, h(z) = 1 2 z 2 , and the target distribution is the spherical Gaussian, N (0, I), then the activation function is the identity, f (z) = z, and the optimization in eq. ( 2) enforces that the second moments of the responses match the second moments of N (0, I), which corresponds to data whitening. If the gains are fixed (η g = 0) and K ≥ N (i.e., synaptic adaptation only), then Alg. 1 corresponds to the adaptive whitening algorithm presented in [34, 53] . Alternatively, if the synaptic weights are fixed (η w = 0) and K ≥ N (N + 1)/2 (i.e., interneuron gain adaptation only), then Alg. 1 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., η g ≫ η w > 0), Alg. 1 corresponds to the multi-timescale adaptive whitening algorithm presented in [52] .When the target distribution is the spherical Gaussian and W is constrained to be an orthogonal matrix, Alg. 1 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 it is more biologically realistic since neural systems are unlikely to be constrained to have orthonormal synaptic weight vectors. 
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. 3 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 correlations 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 to early sensory neuron responses [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 can more efficiently represent natural signals by implementing nonlinear transforms. Consistent with this, nonlinear phenomenological models of neural responses (e.g., divisive normalization [36, 38] ) effectively reduce these higher-order statistical dependencies [37, 59] . However, the circuit mechanisms that support these nonlinear transformations are unknown. We consider the case that circuits nonlinearly transform their inputs to produce Gaussian responses. Gaussian responses may be interpreted as firing rate responses, or as their logs (e.g., membrane potentials are Gaussian, and firing rates are exponentiated). From an efficient coding and computation perspective, Gaussian representations are appealing for a variety of reasons. First, among distributions with given covariance structure, Gaussian distributions have maximum entropy. Therefore, if metabolic demands are a function of the (co)variance 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] [61] [62] , so Gaussianization can be viewed as form of linearization that facilitates downstream computation. In addition, the efficiency of the representation is 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] [65] [66] [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.We computed the responses of a local oriented filter [68] , which captures receptive field selectivity of neurons in primate visual cortex [69] , applied to natural images from the Kodak dataset [70] . These local filter responses are notorious for their sparse heavy-tailed statistical properties that can be well approximated by generalized Gaussian distributions of the form p s (s) ∝ exp(-|s/α| β ), where α is referred to as 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 (α, β) vary across images. While linear methods such as variance normalization are sufficient for rescaling the distribution, adaptive nonlinear transformations are required to reshape these heavy-tailed distributions. 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 between the local filter response pairs (blue line) and ZCA whitened local filter response pairs (orange line) for spatial offsets between d = 2 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. 4 also suggests that ZCA whitening can even slightly increase the mutual information between coordinates-see also, rows d = 8 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.How do we choose the family of activation functions {f (θ, •)}? One approach is to choose a kernel that can approximate a general class of functions. An alternative approach, adopted here, which is motivated by the efficient coding hypothesis [14, 15, 17] , is to choose a family of activation functions that is well matched to the marginal statistics of natural signals. In 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 a family of interneuron activation functions that are optimal for transforming generalized Gaussian distributions with varying (α, β) into the standard Gaussian distribution. When µ = 0, this implies (see Appx. B) that for each choice of scale α and shape β, there is a gain g and parameter θ such thatgf (θ, •) = F -1 α,β • Φ(•), where Φ(•) is the cdf of N (0, 1). However, if we define f (θ, z) in terms of the above display, then we do not have a closed-form solution for ϕ(θ, z) or ∇ θ ϕ(θ, z), which are both required to implement Alg. 1. Instead, we found that F -1 α,β • Φ(•) can be well approximated by the simple algebraic formEQUATIONwhere a(θ) and b(θ) are specified nonnegative functions of θ > 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 (θ, z) = sign(z)|z| θ is optimal for Gaussianizing scalar signals whose marginal tail densities satisfy p s (s) ∝ |s| q-1 exp(-|s| 2q ), where q = 1/θ. This closely resembles the tail densities of generalized Gaussian densities (when α = 1), suggesting monomial activations are effective for shaping the tails.We first apply Alg. 1 in the scalar setting N = K = 1 to demonstrate that our 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, learning rates (η g , η θ ) = (10 -5 , 10 -5 ) and batch size 10 for 10 5 iterations. Fig. 2CD shows the learned interneuron activation functions gf (θ, •) and the learned transforms T (•). Fig. 2E shows histograms of the circuit responses. Compared to the local filter responses, the circuit responses are visually much closer to Gaussian. We found that the circuit performs worse when the distribution p s is more 'peaked' around zero (i.e., when s is sparser and the fitted shape parameter β is smaller), as evidenced by the mismatch between the response distribution p r and the Gaussian distribution N (0, 1) near zero in 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).Next, we apply Alg. 1 in the multivariate setting N = 2. For each image and spatial offset, we ran Alg. 1 with K = 2, 3, 4 interneurons, µ = 0, learning rates (η g , η θ , η w ) = (10 -4 , 10 -6 , 10 -4 ), and batch size 10 for 10 6 iterations. Contour plots of the learned circuit responses for one image are shown in Fig. 3CDE ; see Appx. C for more examples. The mutual information between circuit responses are 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 K = 4 interneurons. For spatial offsets greater than d = 32, the local filter responses already have low mutual information and the circuit does not significantly reduce the mutual information between responses. 
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 neural activation functions. Our model draws inspiration from the extensive neuroscience literature on efficient coding [15] , Hebbian synaptic plasticity [46] , interneuron function [29] and gain modulation [47] . The model proposes complementary roles for different physiological processes: Hebbian synaptic plasticity learns directions that are least matched to the target distribution and interneurons adapt their gains and activation functions to transform the marginal responses along these directions.Our circuit model captures many features of biological circuits, including interneurons, neural adaptation and Hebbian synaptic plasticity, providing a bridge between neural anatomy and physiology and a computational objective inspired by efficient coding. In particular, our model suggests precise relationships between marginal statistics of the signal and interneuron input-output functions. The form of the input-output function for the local interneurons-linear-nonlinear with gain modulationclosely resembles phenomenological models of neurons [73] , and the parameters W, θ, g can potentially be fit to neural recordings and compared with the optimal parameters that can be derived from the signal statistics p s . Furthermore, our model predicts a relation 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 also aspects of our circuit that are not biologically realistic. For example, our model focuses on the role of local interneurons in 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 is consistent with some sensory circuits (e.g., olfactory bulb), but cannot account for excitatory-excitatory connections or inhibitoryinhibitory connections in cortical circuits. Finally, the synaptic weights are not sign-constrained, which violates 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 ≥ 0, which will result in a projected gradient step in Alg. 1; however, the circuit responses will generally be less aligned with the target distribution.A limitation of our experiments is that we only test our method on two-dimensional inputs, demonstrating that three interneurons are sufficient to dramatically reduce the redundancy in the circuit responses. However, natural signals are often very high-dimensional and it is not clear how the number of interneurons required to effectively reduce redundancy will scale with the dimension of the signal. We are optimistic that the number of interneurons required will scale reasonably. 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, which greatly reduces the number of interneurons that are required as the dimension of the input signal grows.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, there are several neural circuit models that implement forms of divisive normalization [74] [75] [76] , a transformation that is optimal for efficient encoding of natural signals [37, 59] . In addition, there is a body of work on normative spiking models derived from objectives which maximize the information encoded per spike [77] [78] [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 as a population level version of the seminal work by Laughlin [16] . 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.Finally, our results may also be relevant beyond the biological setting. Gaussianization and normalizing flows are active areas of research [10, 13, 80] . We offer a novel continuous-learning solution inspired by neuroscience that learns using a combination of weight updates and activation function updates (related to trainable activation functions [81] ). In low-dimensional settings, when the constraint functions are 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 high-dimensional inputs, where the curse of dimensionality presents significant challenges.In this section, we discuss the activation functions used for the experiments carried out in Sec. 4 of the main text.How should we choose the family of activation functions {f (θ, •)}? As stated in the main text, one approach is the choose a 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 formp s (s) = β 2αΓ(1/β) exp(-|s/α| β )with varying scale α and shape β. Here Γ is the gamma function. Therefore, from eq. ( 8), we see that an optimal family of activation functions {f (θ, •)} is such that for each choice of scale α and shape β, there is a gain g and parameter θ such thatgf (θ, z) + µz = F -1 α,β • Φ(z), where F α,β and Φ denote the cdfs of the generalized Gaussian distribution and the standard Gaussian distribution. In general, defining the family of activation functions directly in terms of the above display leads to challenges implementing Alg. 1 since ϕ(θ, z) and ∇ θ ϕ(θ, z) are not readily computable.We instead sought a simple algebraic expression that approximates F -1 α,β • Φ(z). First, suppose the activation function takes the form f (θ, z) = sign(z)|z| θ . When µ = 0, this activation function is optimal for a signal whose cdf F monomial satisfies gsign(z)|z| θ = F -1 monomial • Φ(z). Rearranging and differentiating with respect to z, we find that the pdf of the signal p monomial satsifiesp monomial (gsign(z)|z| θ )gθ|z| θ-1 = 1 √ 2π exp - 1 2 z 2 .Substituting in with s = gsign(z)|z| θ , we find the following expression for the pdfp monomial (s) = q g √ 2π |s/g| q-1 exp - 1 2 |s/g| 2q ,where q = 1/θ. This closely resembles aspects of the pdf for the tails of the generalized Gaussian distribution. To reshape the distribution local when s ≈ 0, we included a linear term, resulting in the activation functionf (θ, z) = a(θ)z + b(θ)sign(z)|z| θ .Here a(θ) and b(θ) are nonnegative functions given by a(θ) = exp((2θ -3.85) 1.95 ) b(θ) = exp(θ 2.32 -5.9).The forms of a(θ) and b(θ) were chosen to approximately minimize min θ max z |F s (gf (θ, z))-Φ(z)| when F s is the cdf of a generalized Gaussian distribution with shape parameter β between 0.2 and 1.Here we derive the updates to the activation functions. Recall from Alg. 1 that the θ update is given by θ ← θ + η θ ∇ θ ϕ(θ, z), whereϕ(θ, z) := h(θ, z) -E z∼N (0,1) [h(θ, z)] = a(θ) 2 (z 2 -1) + b(θ) θ + 1 (|z| θ+1 -C(θ + 1)),and C(p) is the absolute p-moment of a scalar Gaussian random variable:C(p) := E z∼N (0,1) [|z| p ] = 2 p π Γ p + 1 2 .Note that when θ = 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 θ, we obtain the updates∂ϕ(θ, z) ∂θ = a ′ (θ) 2 (z 2 -1) + (θ + 1)b ′ (θ) -b(θ) (θ + 1) 2 (|z| θ+1 -C(θ + 1)) + b(θ) θ + 1 |z| θ+1 log |z| -C ′ (θ + 1) ,whereEQUATIONψ (0) (p) is the polygamma function and we have used the fact that the derivative of the gamma function is Γ ′ (p) = Γ(p)ψ (0) (p).We ran our experiments on a cluster comprised of 40-core Intel Skylake nodes with 768GB of RAM. Each experiment shown in Fig. 2 (i.e., Gaussianization of the coefficients from 1 image) took less than 5 minutes to run. Each experiment shown in Fig. 3 took around 5 hours to run (we did not optimize the choice of hyper-parameters).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. 2 of the main text for a detailed description). We note that when the shape parameter β of the fitted generalized Gaussian distribution is small, the distribution of 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 gf (θ, z) and the optimal activation function F -1 α,β •Φ(z) when z ≈ 0 can lead to large discrepancies between the response distribution and N (0, 1) near zero. This could potentially be resolved by choosing a parameterization of gf (θ, z) that better approximates the optimal activation function.In Fig. 6 , we provide additional examples of contour plots of 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 a detailed description). In some of these examples we see that unlike the standard Gaussian distribution N (0, I) the learned response distribution is clearly non-monotone with respect to ∥r∥. This non-monotonicity is likely inherited from non-optimality of the interneuron activation functions that leads to the discrepancies in the scalar case (see Fig. 5 ). Biologically, we can interpret this as a Gaussian distribution of voltages on cell bodies, which are transformed through an exponential activation function to a positive-valued, heavy-tailed log-normal distribution of firing rates that have been observed in cortex[40].This distance is closely related to (max-)sliced Wasserstein metrics[44,45], which quantify the distance between two distributions in terms of the Wasserstein distances between their marginal distributions. Specifically, when K = 1 and h(θ, •) parameterizes the set of all Lipschitz-1 functions, then the sliced and max-sliced Wasserstein-1 distances between pr and ptarget are E w∼Unif(S N -1 ) [d w,ϕ (pr)] and maxw d w,ϕ (pr).Example code for our experiments can be found at https://github.com/dlipshutz/shaping. 
 
 Appendix:  


Decision: 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.

Reviews: review: paper_summary: paper_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.

Utilizing 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.
strength_weakness: strength_weakness: 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.

2. 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.

3. 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.

4. 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.
questions: 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?

2. 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: 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.
paper_summary: paper_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.
strength_weakness: strength_weakness: 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.
Since the model violates Dale's law, the discussion about the function of interneurons seems out of scope.

Moreover, 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. 

Authors 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.  

There are typos on several places, for example: lines 75, 77, 100, 234
questions: 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?
The 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.

Authors 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?

Is inhibition required to solve this optimization problem, e.g., could the same operation be achieved through plasticity of E neurons?
limitations: limitations: Authors addressed the limitations adequately.
paper_summary: paper_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.
strength_weakness: strength_weakness: - The authors should define $p_{target}$, $p_{marginal}$, and $p_r$ explicitly instead of using descriptive terms for better clarity.
- 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.
- 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.
questions: 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: limitations: The authors have adequately addressed the limitation.
paper_summary: paper_summary: This paper proposes a normative recurrent circuit to solve an optimal transport problem, focusing on the problem of Gaussianization of natural image statistics.
strength_weakness: strength_weakness: 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.
questions: 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.

- 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. 

- In Line 271, should "activations" be plural? I think this is a typo. 

- 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. 

- 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: limitations: As noted above, the authors clearly discuss the limitations of their work.

Comments: title: Thank you for the discussion
comment: Thank you!
I updated the official score.
Likewise, we've appreciated your questions and engagement during the discussion period. 

Request: can you please edit your official score so that your score increase is reflected in the reported average? Thank you!
Thank you for your thoughtful questions and comments, and for your engagement during the discussion period. Based on the discussion, we plan to revise our manuscript to expand our discussion on our model's **biological realism** and its relation to **existing models**. If any additional questions or comments come to mind, please let us know.
Thank you to the authors for your thoughtful reply, which has deepened my understanding of this work. I will raise my score to a 5.
Authors have replied thoroughly to my questions. I have appreciated the discussion around this paper.
The most compelling example is the olfactory bulb, an early stage of olfactory processing in vertebrates, where excitatory mitral cells form dendrodendritic connections with local inhibitory granule cells [Shepherd 2003], leading to symmetric *connectivity* matrices (though not necessarily symmetric *weight* matrices). Within cortical circuits, it's unclear, though somatostatin interneurons that form local, reciprocal connections with pyramidal cells are promising candidates. We view the symmetry constraint as a testable prediction, which can potentially be measured, for example, in recently reported connectomics datasets from mouse visual cortex [Schneider-Mizell bioRxiv 2024].
I have carefully read the comments from other reviewers and appreciate the authors' detailed responses. I have a minor question regarding your discussion on the "Symmetry of synaptic weights" where you mention, "When the weights are decoupled, they converge asymptotically toward a symmetric solution due to the symmetry of their learning rules." While this conclusion is derived from both analytical and numerical analyses, I wonder if such symmetric connections are also observed in real biological systems, especially considering that this section addresses the question of biological plausibility.
I've read the other reviewers' comments and appreciate the authors' thoughtful responses. I'll be keeping my score unchanged.
Thank you for the score bump. We're glad that our proposed revisions have improved the clarity of our work and we will incorporate them in our revised paper.

In response to your questions:

1. This would not change the learned response distribution as this is set by $p_\text{target}$; however it would change the learning algorithm. Specifically, replacing the $L^2$ (quadratic) regularizer $\lambda||T({\bf s})||^2$ in equation (1) with an $L^1$ (linear) regularizer of the form $\lambda|T({\bf s})|$ would encourage sparser neural responses. However, to compensate, the interneurons would adapt to offset the $L^1$ regulariser. Now if the goal is to encourage sparsity, this could built into $p_\text{target}$. For example, rather than choosing $p_\text{target}$ to be Gaussian, it could be chosen to have heavier tails. In this case adding an $L^1$ regularizer to the objective would help the circuit achieve this goal.
2. The main reason for normalization is to prevent the weights from diverging (or collapsing), which is a notorious problem when the synapses update according to a Hebbian update without a homeostatic compensation mechanism; see, e.g., [Abbott & Nelson "Synaptic plasticity: Taming the beast", Nature Neuro. 2000]. As you suggest, this can be achieved by normalizing every $N>1$ steps. Alternatively, this step can be replaced by a dynamic normalization process similar to the one shown in appendix E.1 of [Duong et al. NeurIPS 2023].
Thank you for your reply to my questions. In light of improvements made in the rebuttal I am increasing my score to 5.

I appreciate that you implemented and tested the model that obeys Dale's law - this is an important proof of concept in spite of  reduced performance of the Dalian network. Also, I find that your work can be better appreciated in light of clarifications and improvements described in the general rebuttal. I suggest that these improvements and clarifications are carefully incorporated into the revised paper, in particular, the points about the batch training and symmetry of synaptic weights seem important. 

I have two more questions. 
1) The optimisation problem is formulated with a quadratic regulariser. Is it instead possible to use a linear regulariser? If so, do solutions obtained with a linear regulariser differ from those obtained with a quadratic regulariser?

2) The model uses normalisation of weights after each update. I suppose this is necessary for convergence and I do not have a problem with it. Nevertheless, my question is what are the consequences of not normalising the weights? Also, can good performance be achieved if weights are normalised only once every n>1 updates instead of every update?
Thanks for your thoughtful response to my comments and those of the other reviewers. I maintain my positive assessment. 

A comment: I agree with Reviewer Voaj that it is important to acknowledge past work on normalization and predictive coding with spiking networks. However, I am familiar with these works, and in my opinion the mathematical framing of the present manuscript is certainly sufficiently novel and interesting to warrant publication. I'm not too concerned with violations of Dale's law; this is already an interesting first step.