Unconditional stability of a recurrent neural circuit implementing divisive normalization
Authors: Shivang Rawat
David Heeger
Stefano Martiniani
Venue: NeurIPS 2024 poster (2024)
Keywords: Recurrent Networks
Theoretical Neuroscience
Dynamical Systems
Normalization
Primary area: neuroscience_and_cognitive_science
TL;DR: We prove unconditional stability of a neurodynamical model and demonstrate its performance on sequential datasets.

Abstract: Stability in recurrent neural models poses a significant challenge, particularly in developing biologically plausible neurodynamical models that can be seamlessly trained. Traditional cortical circuit models are notoriously difficult to train due to expansive nonlinearities in the dynamical system, leading to an optimization problem with nonlinear stability constraints that are difficult to impose. Conversely, recurrent neural networks (RNNs) excel in tasks involving sequential data but lack biological plausibility and interpretability. In this work, we address these challenges by linking dynamic divisive normalization (DN) to the stability of "oscillatory recurrent gated neural integrator circuits'' (ORGaNICs), a biologically plausible recurrent cortical circuit model that dynamically achieves DN and that has been shown to simulate a wide range of neurophysiological phenomena. By using the indirect method of Lyapunov, we prove the remarkable property of unconditional local stability for an arbitrary-dimensional ORGaNICs circuit when the recurrent weight matrix is the identity. We thus connect ORGaNICs to a system of coupled damped harmonic oscillators, which enables us to derive the circuit's energy function, providing a normative principle of what the circuit, and individual neurons, aim to accomplish. Further, for a generic recurrent weight matrix, we prove the stability of the 2D model and demonstrate empirically that stability holds in higher dimensions. Finally, we show that ORGaNICs can be trained by backpropagation through time without gradient clipping/scaling, thanks to its intrinsic stability property and adaptive time constants, which address the problems of exploding, vanishing, and oscillating gradients. By evaluating the model's performance on RNN benchmarks, we find that ORGaNICs outperform alternative neurodynamical models on static image classification tasks and perform comparably to LSTMs on sequential tasks.

Unconditional stability of a recurrent neural circuit implementing divisive normalization
Shivang Rawat 1.2 David Heeger
Stefano Martiniani 1,2.5 4.2.6
Courant Institute of Mathematical Sciences, NYU
Center for Soft Matter Research, Department of Physics. NYU Department of Psychology. NYU Center for Neural Science, NYU
Simons Center for Computational Physical Chemistry, Department of Chemistry, NYU
(sr6364, david. heeger, sm7683)@nyu. edu
Abstract
Stability in recurrent neural models poses significant challenge, particularly in developing biologically plausible neurodynamical models that can be seamlessly trained. Traditional cortical circuit models notoriously difficult to train due to expansive nonlinearities in the dynamical system, leading to an optimization problem with nonlinear stability constraints that are difficult impose. Conversely, recurrent neural networks (RNNs) excel in tasks involving sequential data but lack biological plausibility and interpretabillty.. In this work, address these challenges by linking dynamic divisive normalization (DN) to the stability of "oscillatory recurrent gated neural integrator circuits" (ORGaNICs), a biologically plausible recurrent cortical circuit model that dynamically achieves DN and that has been shown to simulate wide range of neurophysiological phenomena. By using the indirect method of Lyapunov, we prove the remarkable property of unconditional local stability for an arbitrary-dimensional ORGaNICs circuit when the recurrent weight matrix is the identity. We thus connect ORGaNIC's to a system of coupled damped harmonic oscillators, which enables us derive the circuit's s energy function, providing normative principle of what the circuit, and individual neurons, aim to accomplish. Further, for generic recurrent weight matrix, we prove the stability of the 2D model and demonstrate empirically that stability holds in higher dimensions. Finally, we show that ORGaNICs can be trained by backpropagation through time without gradient clipping/scaling, thanks to its intrinsic stability property and adaptive time constants, which address the problems of exploding, vanishing, and oscillating gradients. By evaluating the model's performance on RNN benchmarks, we find that ORGaNICs outperform alternative neurodynamical models static image classification tasks and perform comparably to LSTMS on sequential tasks.
Introduction
Deep neural networks (DNNs) have found widespread use in modeling tasks from experimental systems neuroscience. The allure of DNN-based models lies in their ease of training and the flexibility they offer in architecting systems with desired properties [1-3]. In contrast, neurodynamical models like the Filson-Cowan 4 the Stabilized Supralinear Network (SSN) (3) are more biologically plausible than DNNs, but these models confront considerable training challenges due to the lack of stability guarantees for nigh-dimensional problems. Training recurrent neural networks (RNNs), by comparison, is more straightforward thanks to ad hoc regularization techniques like layer nor-
38th Conference on Neural Information Processing Systems (NcurIPS 2024).

---

malization, batch normalization, and gradient clipping/scaling, which help stabilize training without imposing strict stability constraints. Conversely, neurodynamical models require enforcing hard stability constraints while maintaining biological plausibility. In lower dimensions, it is relatively straightforward to derive constraints on model parameters that ensure a dynamically stable sys- tem ( However, for nigh-dimensional systems, this becomes significantly more challenging, as integrating these hard constraints into the optimization problem is more complex { [). Stability is generally advantageous in DNNs, as it is linked to improved generalization, mitigation of exploding gradient problems, increased robustness to input noise, and simplified training techniques [10].
The divisive normalization (DN) model was developed to explain the responses of neurons in the primary visual cortex (V1) (1), and has since been applied to diverse cognitive processes and
neural systems [15 24)) Therefore, DN has been proposed as canonical neural computation (25) that
linked to many well-documented physiological [26] 27 and psychophysical (28) 29 phenomena.
DN models various neural processes: adaptation [30]3] attention (32), automatic gain control [3]]
decorrelation, and statistical whitening [34]. The denning characteristic of DN is that each neuron's
response is divided by weighted sum of the activity of pool of neurons (Eq, 2 below) like when normalizing the length of a vector. Due to its wide applicability and ability to explain a variety of neurophysiologica phenomena, we argue that this characteristic should be central to any neurodynamical model. Both the Wilson-Cowar and SSN models have been shown to approximate DN responses (5 35). but only approximately in certain parameter regimes.
Normalization techniques have been extensively adopted for training DNNs, demonstrating their ability to stabilize, accelerate training, and enhance generalization (36-38). Divisive normalization can be viewed as a comprehensive normalization strategy, with batch and layer normalization being specific instances (39]. Models implementing DN have shown superior performance compared to common normalization methods (Batch (36), Layer [37], Group (40) in tasks such as image recogni- tion with convolutional neural networks (CNNs) [41] and language modeling with RNNs (39 42]
Despite the foundational role of these techniques in deep learning algorithms, their implementation is ad hoc, limiting their conceptual relevance. They serve as practical solutions addressing the limitations of current machine learning frameworks rather than offering principled insights derived from understanding cortical circuits.
has been proposed that DN achieved via a recurrent circuit gated neural integrator circuits (ORGaNICs) are vate-based recurrent neural circuit models that implement DN dynamically via recurrent amplification (47 (8). Since ORGaNICs' response follows the DN equation at steady-state, its steady-state response captures the full range of aforementioned neural phenomena explained by DN (11.44). ORGaNICs have further been shown to simulate key time-dependent neurophysiological and cognitive/perceptual phenomena under realistic biophysical constraints [47 48). Additional phenomena not explained by DN (49) can in principle be integrated into the model. In this paper. however, we focus the effects of DN on the dynamical stability of ORGaNICs. Despite some empirical evidence that ORGaNICs are highly robust, the question of whether the model is stable for arbitrary parameter choices, and thus whether it can be robustly trained on ML tasks by
47
Oscillatory recurrent
(BPTT), remains open.
Here, we establish the unconditional stability - applicable across all parameters and inputs - of a multidimensional
nennon-typess ORGaNICs model when the recurrent weight matrix is the identity. We prove this result, detailed in Section 4 by the indirect method of Lyapunov: we perform linear stability analysis around the model's analytically-known normalization fixed point and reduce the stability problem to that of a high- -ighesinnni mechanical system, whose stability is defined in terms of a tractable quadratic eigenvalue problem. We then address the stability of the model with an arbitrary recurrent weight matrix in Section 5 While the indirect method of Lyapunov becomes intractable for such a system, we provide proof of unconditional stability for a two-dimensional circuit with an arbitrary recurrent weight and offer empirical evidence supporting the claim of stability for high-dimensional systems.
ORGaNICs can be viewed as bicaly plausible extensions of Long Short Term Memory units (LSTMs) (3) and Gated Recurrent Units (GRUs) (50), RNN architectures that have been widely used in ML applications (3 (155)44). The main differences are that ORGaNICs operate in continuous time and have built-in dynamic normalization (via recurrent gain modulation) and built-in attention (via input gain modulation). Thus, expect that ORGaNICs should be able to solve relatively sophisticated tasks (47.. Here, we demonstrate (Section 6 that by virtue of their intrinsic

---

stability, ORGaNIC' can be trained on sequence modeling tasks by BPTT, in the same manner as traditional RNNs (unlike SSN that instead requires costly specialized training strategies (55]), despite implementing power-law activations [5]. Moreover, we show that ORGaNICs trained by naive BPTT (i.e., without gradient clipping/scaling or other ad hoc strategies) achieve performance comparable to LSTMs on the tasks that we consider, despite no systematic hyperparameter tuning.
2 Related Work
Trainable biologically plausible neurodynamical models: There have been several attempts to develop neurodynamical models that mimic the function of biological circuits and that can be trained on cognitive tasks. Song et al. (56) incorporated Dale's law into the vanilla RNN architecture, which was successfully trained across variety of cognitive tasks. Building on this, Soo et al. [57) developed a technique for such RNNs to learn long-term dependencies by using skip connections through time. ORGaNICs is model that is already built on biological principles and learn long-term dependencies intrinsically by tuning the (intrinsic or effective) time constants, therefore it does require the method used in [57]. Soo et al. (55) introduced novel training methodology (dynamics- growth) for SSNs and demonstrated its utility for tasks involving static (time- independent) stimuli. However, this training approach is costly and difficult scale (because SSNs, unlike ORGaNICS are not unconditionally stable), and its applicability on tasks with dynamically changing inputs remains unclear.
Dynamical systems view of RNNs: The stability of continuous- time RNNs has been extensively studied and discussed in comprehensive review by Zhang al. [58). Recent advancements have fo- cused on designing architectures that address the issues of vanishing and exploding gradients, thereby enhancing trainability and performance. A central idea in these designs is achieve better trainability and generalization by ensuring the dynamical stability of the network. Moreover to avoid the problem of vanishing gradients the key idea is to constrain the real part of the eigenvalues of the linearized dynamical system to be close to zero, which facilitates the propagation and retention of information over long durations of time. Chang al. (59) and Erichson et al. (60) achieve this by imposing an antisymmetric constraint on the recurrent weight matrix. Mcanwhile, Rusch et al. [6] 62] propose an architecture based on coupled damped harmonic oscillators, resulting in a second-order system of ordinary differential equations that behaves similarly to how ORGaNICs behave in the vicinity of the normalization fixed point, as we show in Section 4 Despite their impressive performance on various sequential data benchmarks, these models lack biological plausibility due to their use of saturating nonlinearities (instead of normalization) and unrealistic weight parameterizations.
Model description
In its simplest form, the two-neuron-types ORGaNICs model 47 48 can be written as,
with neurons of each type
O y -oy=-y+bozze (1-at) at) (w, (Vy+-Vy>
(I)
Ta 2 à - bj 00 W ((y +y @a+?)
where y E IR" and E R. are the membrane potentials (relative an arbitrary threshold potential that we take to be 0) of the excitatory (y) and inhibitory (a) neurons, evolving according to the dynamical equations defined above with ý and á denoting the time derivatives. The notation O denotes element-wise multiplication of vectors, and squaring, rectification, square-root, and division I
are also performed element-wise. 1 is n-dimensional vector with all entries equal to is the input drive to the circuit and is a weighted sum of the input, x E IR" i.e., z W. The firing rates, y = | and at = v | are rectified ([.)) power functions of the underlying membrane potentials. For the derivation of general model with arbitrary power-law exponents, including the Eq. see Appendix A Note that the term y v serves the purpose of defining mechanism for reconstructing the membrane potential (which can be negative, depending on the sign of the input) from the firing rates y* that are strictly nonnegative. y and y are the firing rates of neurons with complementary receptive fields such that they encode inputs with positive and negative signs, respectively. Note that only one of these neurons fires at a given time. In ORGaNICs, these neurons have single dynamical equation for their membrane potentials, where the sign y
3

---

indicates which neuron is active. Neurons with such complementary (anti-phase) receptive fields are found adjacent each other in the visual cortex [63], and we hypothesize that such complementary neurons are ubiquitous throughout the neocortex. b IRt"" and ba E T are the input gains for the external inputs z and fed to neurons y and a, respectively. R is the set of positive real numbers, {i E RR x > 03. I determines the semisaturation of the responses of neurons y by contributing to the depolarization of neurons a. Ty R and Ta B represent the time constants of and a neurons.
In addition to receiving external inputs, both y and neurons receive recurrent inputs, represented by the last term in both of the equations. W, E R= is the recurrent weight matrix that captures lateral connections between the y neurons. This recurrent input is gated by the neurons, via the term (1 at). Similarly, the nonnegative normalization weight matrix, W R, nx* encapsulates the recurrent inputs received by the a neurons. The differential equations are designed in such way that when W, and = bu bol (i.e., with all elements equal to a constant by), the principal neurons follow the normalization equation exactly (and approximately when W, A I) at steady-state
2
1
y vy.1 0² +
(2)
1-2
[z]² 2 and |-z]² represent the contribution of neurons with complementary receptive fields to the normalization pool, and 1- z² is the contrast energy of the input. Note that the recurrent gain, (1 - ). is particular nonlinear function of the output responses/activation designed to achieve DN, while the input gain, b an input gate that can implement an attention mechanism.
4 Stability analysis of high-dimensional model with identity recurrent weights
We consider the stability of the general nigh-dimensional ORGaNICs (Eq. [
when the recurrent weight matrix is identity, W, = We first simplify the dynamical system by noting that / y vy = y and y + +y = y2 yielding the following equations,
Ty y
a b
(3)
Ta ④ -a by 0² w (y² の [al)
For identity recurrent weights, have unique fixed point, given by,
b
y. b O W ( x²)
a = bj o² W (b²
(4)
Since the normalization weights in the matrix W are nonnegative, at steady-state we have a, > 0, that [[a.] v ag. and the corresponding firing rates at steady-state are,
+b 3 2
y b W (b2 1 x2)
a = bB 0² W (b² Oz?)
(5)
Note that recover the normalization equation, Eq. if bo bg Since the fixed points of y and a neurons are known analytically, to prove that this fixed point is locally asymptotically stable (i.e., the responses converge asymptotically the fixed point), we apply the indirect method of Lyapunov at this fixed point (64). This method allows us to analyze the stability of the nonlinear system in the vicinity of the fixed point by studying the corresponding linearized system. The Jacobian matrix J E I8 x22 about (y,, a,), defining the linearized system, is given by.
D
D
Ya
20VATOT,
J
(6)
D
W WD T
2
D ( (-I+ WD (y2))
where (xx) is diagonal matrix of appropriate size with the elements of the vector the diagonal. necessary and sufficient condition for local stability is that the real parts of all cigenvalues this matrix are negative. We thus proceed by computing the characteristic polynomial for the Jacobian,
Pa (A) = dJt(J- - AI). The roots this polynomial, found by setting pa (A) = O, are the eigenvalues of the system. Consider the block matrix,
An
A1
D
ja
II
-D
20 OT,
Y
J AI=
(7)
Ag
A22
D
2
W (a. Gy.) D (
(-I+ WD

---

Notice that A1 and 12 are diagonal and therefore they commute, i.e., A11A12 A11 we have that det AI) det(As 12 12) which is property the determinant of block matrices with commuting blocks (65] Therefore, the characteristic polynomial of the linearized system after expansion of the terms and simplification is given by,
det(J - XI) det 1 D
+D
+
a Ty
D To
WD
(y:))
+D
a
Ta
Ty Ta
Finding the roots of this polynomial is thus a quadratic eigenvalue problem of the form ((A) = del()²1 + XB K) O. which has been studied extensively (66-69)] L()) be interpreted as the characteristic polynomial associated with a system of linear second-orde differential equations with constant coefficients of the form I& Bx + Kx Therefore, proving the stability of our system (i.e., Re(A) < 0 for {) ((A) = 0}), is equivalent to proving the asymptotic stability of IX Bx Kx 0.
Tisseur et al. [67] and Kirillov et al. [69] list a set of constraints on the damping matrix, B, and stiffness matrix, K, that yield stable system, but they not directly applicable In the context of high-dimensional mechanical system, our system falls under the category of gyroscopically stabilized systems with indefinite damping. Few results are known about the conditions leading to the stability of such systems. By constructing a Lyapunov function, we prove (Appendix ( the following stability theorem that directly applicable our system, following an approach similar to Kliem et al. (70).
system.
Theorem 4.1. For a system of linear differential equations with constant coefficients of the form,
I + Bx Kx=0
where B Tx and K E Raxx a positive diagonal matrix (hence > 0), the dynamical system is globally asymptotically stable if B is Lyapunov diagonally stable.
Since the stiffness matrix,
aa Ty
O
K =
= =
(10)
TYOTA
positive diagonal matrix, sufficient condition for stability of the system is that the damping matrix, B, given by,
B BB B2
B3
a2
=D
D
D
WD (y?)
Ta
Ty
Ta
(11)
D/
0²
W (b2
x2)
1 Ta
ba о W
z (b² t²)
=D
D
WD
Ta
Ty
b
is Lyapunov diagonally stable, i.e., there exists positive definite diagonal matrix T, such that TB -B T positive definite.
Since all of the parameters are positive, and the weights in the matrix W are nonnegative, we can conclude the following: B, and B, are positive diagonal matrices and B3 is matrix with all positive entries (that may or may not be symmetric). Therefore, B is a Z-matrix, meaning that its dia-diago entries nonpositive. Further, Z-matrix is Lyapunov diagonally stable if and only if it is nonsingular M-matrix. Intuitively, M-matrices are matrices with non-positive off-diagonal elements and "large enough" positive diagonal entries. Berman & Plemmons (71) list 50 equivalent definitions of nonsingular M-matrices. We use the one that is best suited for our problem,
Theorem 4.2. (Chapter 6, Theorem 2.3 from (71)) Z-matrix matrix B E RRxx Lyapunov diagonally stable if and only if there exists convergent regular splitting of the matrix, that is, it has a representation of the form
M N, where M- and N have all nonnegative entries, and M ' N has spectral radius smaller than 1.
We now show that, indeed, B has a convergent regular splitting for all combinations of the circuit parameters and for all inputs. We have already shown that B is Z-matrix, therefore, the first
5

---

condition of the theorem is satisfied. Next, we consider the following splitting B = M - N with M B, + B2 and N B3. Since B and B2 are positive diagonal matrices, M is nonnegative, while N is also nonnegative because B3 has all positive entries. Therefore, the only condition left to satisfy is that the spectral radius of M- is smaller than or that the matrix convergent.
The matrix S = M = (Bi + B₂) BB can be written as,
1
ba
Z
==
WD b
(12)
1 (Ta/TTT)) vbj o o² W (b²
o² +
W (b² 3 2z2)
We prove the following theorem (Appendix D which directly applies to S,
Theorem 4.3. A matrix A of the form A D(t) W D (u/ (v+ Wu)) is convergent (i.e., its spectral rudius is less than 1), if W E Raxn and t, u, V R2 satisfy < 4 < 1, u 2 0/ and Wij 2 ffr all j.
Defining t + 1/(1 +
@ O u and v - bj @ o².
can be seen that they satisfy the constraints of the theorem, and thus is convergent. This implies that B has convergent regular splitring and, as a result, the linearized dynamical system is unconditionally globally asymptotically stable for all the values parameters and inputs. Further, the global asymptotic stability of linearization implies the local asymptotic stability of the normalization fixed point for ORGaNICs.
This result holds even when the neurons have different time constants, regardless of their type, as no assumptions were made about the time constants. This finding is significant for machine learning, particularly for designing architectures based on ORGaNICs. It allows neurons/units to integrate information at varying time scales while maintaining a stable circuit that performs normalization dynamically. Moreover, analytical expressions for eigenvalues can be obtained in the following case, Theorem 4.4. Let W, = I, the normalization matrix be given by W = a.E, where E is the all-ones matrix, and the parameters are scalars, i.e., Ty = Ty To = Ta by = by and a = 01. Under these conditions, the eigenvalues of the system admit closed form solutions (detailed in Appendix
This result is particularly useful for neuroscience it elucidates the connection between ORGaNICs parameters and the strength and frequency of oscillatory activity. Since we followed a direct Lyapumov approach to prove Theorem 4.1 as shown in Appendix B we can derive an energy (viz., Lyapunov function) for ORGaNIC's as shown in Appendix H
Theorem 4.5. When W, = I. the energy (Lyapunov function) minimized by ORGaNICs in the vicinity of the normalization fixed point, is given by,
a(y,a) a)
Osi
Tu
uei i (ys vos)2 + (Va.y: VasiMes))
(13)
Yoi
Toi
i=1
Where t; are the diagonal entries T and you (asi) are the steady-state values of neurons ys (a;).
Specifically, for a two-dimensional model (one neuron and one neuron) this expression simplifies to reveal that ORGaNICs behave like a damped harmonic oscillator with energy,
Vly a) Ty 0 xbb2²
bz
+ (Vay - bz)²
(14)
0602
wb222
This result demonstrates that ORGaNICs minimize the residual of the instantancously reconstructed gated input drive (Vay- - bz), while also ensuring that the principal neuron's 's response, y, achieves DN. The balance between these objectives is governed by the parameters and the external input strength. With fixed parameters, weaker inputs, 2, cause the model to prioritize input matching over normalization, whereas stronger inputs inereasingly engage the normalization objective.
5 Stability analysis for arbitrary recurrent weights
Now, we relax the constraint that the recurrent weight matrix must be identity, allowing W,, # I, and see how the stability result changes. This leads to the following set of equations,
T - Oz
Vial) @ (Wry)
(15)
by 2 W (y² [a])
Ta a

---

0 < E 1
Unstable fixed point
Stable fixed point
bao 1/w,
byo 1 1/w,
2-
y y
a
1
y
Figure 1: Phase portraits for 2D ORGaNICs with positive input drive. We plot the phase portraits of 2D ORGaNICs in the vicinity of the stable fixed points for contractive (a, d) and expansive (b, c, e, recurrence scalar w,. A stable fixed point always exists, regardless of the parameter values. (a-c), The main model (Eq. T6). (d-f), The rectified model (Eq, T02). Red stars and black circles indicate stable and unstable fixed points, respectively. The parameters for all plots are: b 0.5, Ta 2 ms, Ty ms, w = 1.0, and 1.0. For (a) & (d), the parameters are Wr = 0.5, bo = 0.5, σ 0.1; for (b) & (e), W. 2.0, 0.5, a 0.1; and for (c) & (f), w, 2.0, bo 1.0, 1.0.
The linear stability analysis becomes intractable for a general W, because we no longer have closed-form analytical expression for the steady states of and a. Additionally, the characteristic polynomial cannot be expressed in way similar Eq 8 Nevertheless, for a two-dimensional system,
T 1リ + bz +
VLal
w.y
(16)
T_ -a + bāo² + wyy
we can prove the following, with a detailed analysis provided in Appendix E
Theorem 5.1. Given that the recurrence contracting, i.e., wy S when 2 0 0) there exists a unique fixed point with > 0 (y. < 0) and > O, and it is asymptotically stable.
Theorem 5.2. Given that the recurrence is expansive, i.e., Wr > 1, there are either 1 or 3 fixed points of which at least one asymptotically stable. When z > 0 (2 0) there exists exactly I fixed point with ys > 0 (y, < 0) and a, > O, and it is asymptotically stable. If byo > 1- 1/w,, there are no additional fixed points. If bbo 1 1/w,, there exist either 0 or additional fixed points with < 0 (ys > 0) and a, >> 0 whose stability cannot be guaranteed.
We plot the phase portraits for these different cases in Fig. The key takeaway is that there is always a fixed point (ye, a.) with a, > O and ys having the same sign as 2. This fixed point is asymptotically stable regardless of the value of w,. Based these results and the proven stability of arbitrary dimensional ORGaNICs when W, = (as shown in Section L. we conjecture that
Conjecture 5.3. Consider high-dimensional ORGaNICs with an arbitrary recurrent weight matrix W, and constraints on the remaining parameters. If the norm of the input drive satisfies ||||| S 1, and the maximum singular value of W, is constrained to be 1, then the system possesses at least one asymptotically stable fixed point.
This conjecture supported by empirical evidence showing consistent stability, as ORGaNICs initialized with random parameters and inputs under these constraints have exhibited stability in 100%
7

---

Algorithm Iterative scheme for the fixed point when the maximum singular value of W, is 1
Input: ORGaNIC's parameters, input (z), tolerance (e), maximum iterations (N) 2: Output: Approximation to the fixed point (y.,
3: <- o² b + W (W, (b @ 2))
I/ initial approximation for neurons initial approximation for neurons
4: y + (W, (b 2 z)) / /va
5: k+ 0
6: while |ly by-bor o 2 (1 - Va) O (W,y) and k < N do
7:
y [
W, + +D (Va) W.) (boz)
y neurons update
8:
a bj c o² W (y² *a)
a neurons update
9:
k 1-
10:
end while
11:
return ay
of trials, see Fig. We further speculate that ORGaNICs may be typically stable beyond this regime as we find that 100% trials yield stable circuit when the constraint on the maximum singular value of W, is increased to 2, but it becomes unstable when it is increased to 3.
6 Experiments
We provide further empirical evidence in support of Conjecture 5.3 that ORGaNIC' is asymptotically stable by showing that stability preserved when training ORGaNIC's using naive BPTT on two different tasks: I) static classification of MNIST, 2) sequential classification of pixel-by-pixel MNIST. Because these ML tasks have no relevance for neurobiological or cognitive processes, we relax one aspect of the biological plausibility of ORGaNICs, specifically, allowing arbitrary (learned) nonnegative values for the intrinsic time constants
6.1 Static input classification task
We first show that we can train ORGaNICs on the MNIST handwritten digit dataset (72 presented to the circuit a static input. This setting corresponds to evolving the responses of the neurons dynamically until they reach a fixed point solution and using the steady-state firing rates of the principal neurons predict the labels, akin to deep equilibrium models [73] While the fixed point of the circuit is known when W, I (given by Eq. 89). we allow W, to be learnable and parameterized it to have a maximum singular value of This constraint allows us to find the fixed point responses of all the neurons without simulation, using a fixed point iteration scheme (Algorithm
Table 1: Test accuracy MNIST dataset
Model
Accuracy
SSN (50:50)
94.9%
SSN (80:20)
95.2%
MIP (50)
98.2%
ORGaNIC' (50:50)
98.1 %
ORGaNICs (80:80)
98.2%
ORGaNICs (two layers)
98.1% %
D that converges with great accuracy in a few (less than 5) steps, see Fig. 4 & 5 We provide an intuition for why this algorithm works with empirical evidence of fast convergence in Appendix G
We trained ORGaNICs on this task (details provided in Appendix I.I and compared its performance SSN (5) trained by dynamics-neutral growth [55). We found that ORGaNIC' perform better than SSN with the same model size, and on par with an MLP (T. We analyzed the eigenvalues of the Jacobian matrix the trained model and consistently found the largest real part to be negative (Fig. 5 indicating stability. Moreover, we found that stability was maintained during training (Fig. 6).
6.2 Time varying input
We trained unconstrained ORGaNICS by naive BPTT on classification task of sequential MNIST (sMNIST), proposed by Le et al. (74). This is challenging task because it involves long-term dependencies and requires the architecture to maintain and integrate information over long timescales. Briefly, the task involves the presentation of pixels of MNIST images sequentially (one pixel for each
Python code for this study is available

---

Table 2: Test accuracy on sequential pixel-by-pixel MNIST and permuted MNIST
Model
sMNIST
PSMNIST
# units
# params
LSTMs 75
97.3%
92.6%
128
68k
AntisymmetricRNN 59
98.0%
95.8%
128
10k
CoRNN [61]
99.3%
96.6%
128
34k
Lipschitz RNN
60
99.4%
96.3%
128
34k
ORGaNICss (fixed time constants)
90.3%
80.3%
64
26k
ORGaNICs (fixed time constants)
94.8%
84.8%
128
100k
ORGaNICs
97.7%
89.9%
64
26k
ORGaNICs
97.8%
90.7%
128
100k
timestep) in scanline order, and at the end of the input the model has predict the digit that was presented. There is more complicated version of this task, permuted sequential MNIST, in which pixels of images are permuted in some random order before being presented sequentially. We train ORGaNICs with different hidden layer sizes (number of neurons) these two tasks by discretizing the rectified ORGaNICs with arbitrary recurrence, Eq, [87] which has all the properties that we have derived for the main model. Since unstable fixed point is undesirable in such task, as it may lead to diverging trajectories, we prefer the rectified model (Appendix (F) over the main model. We proved that the 2D rectified ORGaNICs (Eq. T02) does not exhibit an unstable fixed point for positive inputs, as can also be seen in Fig I The hidden states of the neurons are initialized with uniform random distribution (for more details, Appendix 1.2p. Additionally, we make the input gains and by dynamical with their ODEs given by,
Tb O - + f(W + Wogy W Wb
(17)
The @ bq bo + S(W baz W boyy
W boaa)
We achieved slightly better performance than LSTMs on sMNIST with smaller model size and comparable performance on permuted sMNIST, without hyperparameter optimization and without gradient clipping/scaling (Table 2). We found that the trajectories of are bounded when it is trained on the sequential task (Fig. 7 indicating stability. We also show that the training of ORGaNICs stable and does not require gradient clipping when the intrinsic time constants of the neurons are fixed (Table 2)
Discussion
Summary: While extensive research has been aimed at identifying highly expressive RNN ar- chitectures that can model complex data, there has been little advancement in developing robust, biologically plausible recurrent neural circuits that are easy to train and perform comparably to their artificial counterparts. Regularization techniques such batch, group, and layer normalization have been developed and are implemented as ad hoc add-ons making them biologically implausible. In this work, we bridge these gaps by leveraging the recently proposed ORGaNICs model which implements divisive normalization (DN) dynamically in recurrent circuit. We establish the unconditional stabil- ity of an arbitrary dimensional ORGaNICs circuit with an identity recurrent weight matrix (W,), with all of the other parameters and inputs unconstrained, and provide empirical evidence of stability ORGaNICs with arbitrary W,. Since ORGaNICs remain stable for all parameter values and inputs, we do not need to resort to techniques that are restrictive in parameter space, that require designing unrealistic structures for weight matrices. ORGaNICs' intrinsic stability mitigates the issues of exploding and oscillating gradients, enabling the use of "vanilla" " BPTT without the need for gradient clipping, which is instead required when training LSTMs. Moreover, ORGaNICs effectively address the vanishing gradient problem often encountered when training RNNs. This is achieved by processing information across various timescales, resulting in a blend of lossy and non-lossy neurons, while preserving stability. The model's effectiveness in overcoming vanishing gradients is further evidenced by its competitive performance against architectures specifically designed to address this issue, such LSTMs.
Dynamic normalization: Normalization techniques, such as batch and layer normalization, are fundamental in modern ML architectures significantly enhancing the training and performance of

---

CNNs. However, a principled approach to incorporating normalization into RNNs has remained elusive. While layer normalization is commonly applied RNNs to stabilize training, it does not influence the underlying circuit dynamics since it is applied a-posteriori to the output activations, leaving the stability of RNNs unaffected. Furthermore, DN has been shown to generalize batch and layer normalization [39], leading to improved performance [39] 41 42). ORGaNICs, unlike RNNs with layer normalization, implement DN dynamically within the circuit, marking the first instance of this concept being applied and analyzed in ML. Our work demonstrates that embedding DN within a circuit naturally leads to stability, which is greatly advantageous for trainability. This stability, a consequence of dynamic DN, sets ORGaNICs apart from other RNNs by providing both output normalization and model robustness. As a result, ORGaNICs can be trained using BPTT, achieving performance on par with LSTMs. The key insight is that the dynamic application of DN not only enhances training efficiency but also improves model robustness. This illustrates how the incorporation of neurobiological principles can drive advances in MI..
Interpretability: In the proof of stability, we establish direct connection between ORGaNICs and systems of coupled damped harmonic oscillators, which have long been studied in mechanics and control theory. This analogy not only enables us to derive an interpretable energy function for ORGaNICs (Eq. 13 providing normative principle of what the circuit aims to accomplish, but also sheds light on the link between normalization and dynamical stability of neural circuits. For a relevant ML task, having an analytical expression for the energy function allows to quantify the relative contributions of the individual neurons in the trained model, offering more interpretability than other RNN architectures. For instance, Eq. 13 shows that the ratio of time constants (Ty/Ta) for E-I neuron pairs determines how much weight a neuron assigns to divisive normalization relative to aligning its responses with the input drive z. This insight provides clear functional role for each neuron in the trained model. Moreover, since ORGaNICs are biologically plausible, we can understand how the various components of the dynamical system might be computed within neural circuit (48), bridging the gap between theoretical models and biological implementation, and offering means to generate and test hypotheses about neural computation in real biological systems (which we will be reporting elsewhere).
Although the stability property pertains to continuous-timg time system of nonlinear differential equations, typical implementations for tasks with sequential data involve an Euler dis- cretization of these equations for training purposes. This might lead to stiff dynamical system, potentially causing numerical instabilities and explosive dynamics, highlighting the importance of carefully parameterizing time constants and choosing small enough time step maintain stable dynamics. The proof of unconditional stability is only tractable for the two-dimensional circuit and the high-dimensional circuit with W, , I. Therefore, we can only conjecture the stability of ORGaNICs for arbitrary W,, based on these two limiting cases and on empirical evidence. In the current form, the weight matrices of the input gain modulators, Weg. Wbo. Wbog. and W boa. are each n X n. As result, the number of parameters grows more rapidly with the hidden state size compared to other RNNs. To mitigate this, we plan to explore using compact and/or convolutional weights to prevent a significant increase in the number of parameters as the hidden state size expands.
Limitations:
Attention mechanisms in ORGaNICs: ORGaNIC's have built-in mechanism for attention: mod- ulating the input gain b (e.g., Eq. 17). coupled with DN. This attention mechanism aligns with experimental data on both increases the gain of neural responses and improvements in behavioral performance [19 20 32 76)85). Moreover, this mechanism performs a computation that is analogous to that of attention head in ML systems (including transformers (2) both operate by changing the gain over time. In ORGaNICs, DN replaces the softmax operation typically used in an attention head.
Future work: This study has explored only a single layer of ORGaNICs for the sequential tasks. Future work will examine how stacked layers with feedback connections, similar to those in the cortex, perform on benchmarks for sequential modeling and also on cognitive tasks with long-term dependencies. We have thus far shown that ORGaNIC' can address the problem of long-term dependencies by learning intrinsic time constants. Future investigations will assess the performance of ORGaNIC's for tasks with long-term dependencies by learning to modulate the responses of the and b neurons to control the effective time constant of the recurrent eircuit (without changing the intrinsic time constants) 47). i.e., implementing working memory circuit capable of learning to maintain and manipulate information across various timescales.
10

---

Acknowledgments and Disclosure of Funding
The authors thank anonymous reviewers for their insightful suggestions. The authors also acknowl- edge valuable discussions with Flaviano Morone, Asit Pal, Mathias Casiulis, and Guanming Zhang. This work was supported by the National Institute of Health under award number RO1EY035242. S.M. acknowledges the Simons Center for Computational Physical Chemistry for financial support. This work was supported in part through the NYU IT High Performance Computing resources, services, and staff expertise.
References
[1] Alex Krizhevsky, Ilya Sutskever, and Geoffrey E Hinton. Imagenet classification with deep convolutional neural networks. Advances in neural information processing systems, 25, 2012.
[2] Ashish Vaswani, Noam Shazeer, Niki Parmar, Jakob Uszkoreit, Llion Jones, Aidan N Gomez, Eukasz Kaiser, and Illia Polosukhin. Attention is all you need. Advances in neural information processing systems, 30, 2017.
[3] Sepp Hochreiter and Jürgen Schmidhuber. Long short-term memory. Neural computation. 1735 1997.
[4] Hugh Wilson and Jack Cowan. Excitatory and inhibitory interactions localized popula- tions of model neurons. Biophysical journal, 12(1):1-24, 1972.
[5] Daniel B Rubin, Stephen D) Van Hooser, and Kenneth D Miller. The stabilized supralinear network: a unifying circuit motif underlying multi-input integration in sensory cortex. Neuron, 85(2):02-447, 2015.
[6] Jack D Cowan, Jeremy Neuman, and Wim van Drongelen. Wilson-cowan equations for neocortical dynamics. The Journal of Mathematical Newroscience, 6:1-24, 2016.
[7] Yashar Ahmadian, Daniel B Rubin, and Kenneth Miller. Analysis of the stabilized supralinear network. Neural computation, 25(8):1994 2037, 2013.
[8] James Kotary, Ferdinando Fioretto, Pascal Van Hentenryck, and Bryan Wilder. End-to-end constrained optimization learning: A survey. arXiv preprint arXiv:2103: 16378, 2021.
[9] Priya L Donti, David Rolnick, and Zico Kolter. De3: A learning method for optimization with hard constraints. arXiv preprint arXiv:2104. 12225 2021.
[10] Eldad Haber and Lars Ruthotto. Stable architectures for deep neural networks. Inverse problems,
34( (1):0144004 2017.
(11) Matteo Carandini and David J Heeger. Summation and division by neurons in primate visual cortex. Science, 264(5163): 13316333:333313333 333-1336, 1994.
[12) Duane G Albrecht and Wilson S Geisler. Motion selectivity and the contrast-response function of simple cells in the visual cortex. Visual neuroscience, 7(6):531-546, 1991.
[13] David J Heeger. Normalization of cell responses in cat striate cortex. Visual neuroscience. 92):11111177 1992.
[14] David Heeger. Modeling simple-cell direction selectivity with normalized, half-squared, linear operators. Journal of neurophysiology, 70(5):1885-1888 1993.
[15] Eero P Simoncelli and David Heeger. A model of neuronal responses in visual area mt. Vision research, 3843774,7761, 1998.
[16] John H Reynolds, Leonardo Chelazzi, and Robert Desimone. Competitive mechanisms subserve attention in macaque areas v2 and v4. Journal of Neuroscience, 1736):1736 1999.
[17] Odelia Schwartz and Eero Simoncelli. Natural sound statistics and divisive normalization the auditory system. Advances in neural information processing systems, 13, 2000.
[18] Davide Zoccolan, David D Cox, and James J DiCarlo. Multiple object response normalization monkey inferotemporal cortex. Journal of Neuroscience, (5(36):8150-8168 2005.
[19] Geoffrey M Boynton. framework for describing the effects of attention on visual responses. Vision research, 49(10): 1129-1143, 2009.

---

[20] Joonyeol Lee and John HR Maunsell. A normalization model of attentional modulation of single unit responses. PloS one, 4()): 2009.
[21] Paul M Bays. Noise in neural populations accounts for emrors in working memory. Journal of Neuroscience, 34(10):3632-36)) 36443365,36 2014.
[22] Wei Ma, Masud Husain, and Paul M Bays. Changing concepts working memory. Nature neuroscience, 17(3y:347-356, 2014.
[23] Barbara Zenger-Landolt and David Heeger. Response suppression in v1 agrees with psy- chophysics of surround masking. Journal of Neuroscience, 23(17))6884 (117))6888-683,, 2003.
[24] John M Foley. Human luminance pattern-vision mechanisms: masking experiments require new model. JOSA A, 11(6):1710-177 1711101-9111111 1994.
[25] Matteo Carandini and David Heeger. Normalization canonical neural computation. Nature reviews neuroscience, 13(1):51-62, 2012.
[26] Gijs Joost Brouwer and David Heeger. Cross-orientation suppression in human visual cortex. Journal of neurophysiology, 106(5):2108-211 2011.
[27] James R Cavanaugh, Wyeth Bair, and J Anthony Movshon. Nature and interaction of signals from the receptive field center and surround in macaque v1 neurons. Journal of neurophysiology, 88(5):25302255 22500-546, 2002.
[28] Jing Xing and David Heeger. Serteur interactions in foveal and peripheral vision. Vision research, 40(22):3065-300 -3772, 2000.
[29] Alexander A Petrov, Barbara Anne Dosher, and Zhong-Lin Lu. The dynamics of perceptual learning: incremental reweighting model. Psychological review, 112(4):715, 2005.
[30] Martin J Wainwright, Odelia Schwartz, and Eero P Simoncelli. Natural image statistics and divisive normalization. 2002.
[31] Zachary M Westrick, David Heeger, and Michael S Landy. Pattern adaptation and normaliza- tion reweighting. Journal of Neuroscience, 36(38):9805-9811 2016.
[32] John H Reynolds and David Heeger. The normalization model of attention. Neuron, 61 185, 2009.
[33] David Heeger, Eero P Simoncelli, and J Anthony Movshon. Computational models of cortical visual processing. Proceedings of the National Academy of Sciences, 93(2):623-62 1996.
[34] Siwei Lyu and Eero P Simoncelli. Nonlinear extraction of independent components of natural images using radial gaussianization. Neural computation, (1(6):1485-155) 2009.
[35] Jesús Malo, José Juan Esteve- Taboada, and Marcelo Bertalmío. Cortical divisive normalization from wilson-cowan neural dynamics. Journal of Nonlinear Science, 34(2):1-36, 2024.
[36] Sergey Ioffe and Christian Szegedy. Batch normalization: Accelerating deep network training by reducing internal covariate shift. In International conference on machine learning, pages 448-456 pmlr, 2015.
[37] Jimmy Lei Ba, Jamie Ryan Kiros, and Geoffrey E Hinton. Layer normalization. arXiv preprint 2016.
[38] Lei Huang, Jie Qin, Yi Zhou, Fan Zhu, Li Liu, and Ling Shao. Normalization techniques in training dnns: Methodology. analysis and application. IEEE transactions on puttern analysis and machine intelligence, 1013)1101,, 2023.
[39] Mengye Ren, Renjie Liao, Raquel Urtasun, Fabian H Sinz, and Richard Zemel. Normalizing the normalizers: Comparing and extending network normalization schemes. arXiv preprint 161 1.04220, 2016.
[40] Yuxin Wu and Kaiming He. Group normalization. In Proceedings of the European conference on computer vision (ECCV), pages 3-19, 2018
[41] Michelle Miller, Yue Chung, and Kenneth D Miller. Divisive feature normalization improves image recognition performance in alexnet. In International Conference on Learning Representations, 2021.
[42] Saurabh Singh and Shankar Krishnan. Filter response normalization layer: Eliminating batch dependence in the training of deep neural networks. In Proceedings of the IEEE/CVF conference on computer vision and pattern recognition, pages 1 11237-11246, 2020.
12

---

[43] Matteo Carandini, David Heeger, and J Anthony Movshon. Linearity and normalization in simple cells of the macaque primary visual cortex. Journal of Neuroscience, 17(21):8621-886 1997.
[44] Hirofumi Ozeki, Ian M Finn, Evan S Schaffer, Kenneth D Miller, and David Ferster. Inhibitory stabilization of the cortical network underlies visual surround suppression. Neuron, 62(4):578- 592, 2009.
[45] Matteo Carandini, David Heeger, and Walter Senn. A synaptic explanation of suppression in visual cortex. Journal of Neuroscience, 22(22):10033-1006) 2002.
[46] Tobias Brosch and Heiko Neumann. Interaction of feedforward and feedback streams in visual cortex in firing-rate model of columnar computations. Neural Networks, 54:11-16, 2014.
[47] David J Heeger and Wayne E Mackey. Oscillatory recurrent gated neural integrator circuits (organics), unifying theoretical framework for neural dynamies. Proceedings of the National Academy of Sciences, 1 116(45):22733-2277
2019.
[48] David J Heeger and Klavdia O Zemlianova. A recurrent circuit implements normalization, simulating the dynamics of v1 activity. Proceedings of the National Academy of Sciences, 117036)22444--22505 2020.
[49] Lyndon Duong, Eero Simoncelli, Dmitri Chklovskii, and David Lipshutz. Adaptive whiten- ing with fast gain modulation and slow synaptic plasticity. Advances in Neural Information Processing Systems, 36, 2024.
[50] Junyoung Chung, Caglar Gulcehre, KyungHyun Cho, and Yoshua Bengio. Empirical evaluation of gated recurrent neural networks on sequence modeling. arXiv preprint arXiv:1412.3555
2014.
[51] lly Sutskever, Oriol Vinyals, and Quoc V Le. Sequence to sequence learning with neural networks. Advances in neural information processing systems, 27, 2014.
[52] Kyunghyun Cho, Bart Van Merriënboer, Caglar Gulcehre. Dzmitry Bahdanau, Fethi Bougares, Holger Schwenk, and Yoshua Bengio. Learning phrase representations using rnn encoder- decoder for statistical machine translation. arXiv preprint arXiv: 1078, 2014.
[53] Alex Graves, Abdel-rahman Mohamed, and Geoffrey Hinton. Speech recognition with deep recurrent neural networks. In 2013 IEEE international conference on acoustics, speech and signal processing, pages 6645-6649 Iece, 2013.
[54] Alex Graves. Generating sequences with recurrent neural networks. arXiv preprint arXiv:1308.000 0000000000000..0. 2013.
[55] Wayne Soo and Máté Lengyel. Training stochastic stabilized supralinear networks by dynamics- neutral growth. Advances in Neural Information Processing Systems, 35:29278-2995 292911, 2022.
[56] Francis Song, Guangyu R Yang, and Xiao-Jing Wang. Training secitatory-iniiiiton recurrent neural networks for cognitive tasks: a simple and flexible framework. PLoS computational biology, 1 2016.
[57] Wayne Soo, Vishwa Goudar, and Xiao-Jing Wang. Training biologically plausible recurrent neu- ral networks on cognitive tasks with long-term dependencies. Advances in Neural Information Processing Systems, 36, 2024.
[58] Huaguang Zhang, Zhanshan Wang, and Derong Liu. A comprehensive review of stability analysis of continuous-time recurrent neural networks. IEEE Transactions on Neural Networks and Learning Systems, 1229-122222 2014.
[59] Bo Chang, Minmin Chen, Eldad Haber, and Ed H Chi. Antisymmetricrn: A dynamical system view on recurrent neural networks. arXiv preprint /902 2019.
[60] N Benjamin Erichson, Omri Azencot, Alejandro Queiruga, Liam Hodgkinson, and Michael W Mahoney. Lipschitz recurrent neural networks. arXiv preprint arXiv:2006. 12000000000000 2020.
[61] T Konstantin Rusch and Siddhartha Mishra. Coupled oscillatory recurrent neural network (cornn): An accurate and (gradient) stable architecture for learning long time dependencies. arXiv preprint arXiv:2010.0095/ 2020.
[62] T Konstantin Rusch and Siddhartha Mishra. Unicornn: A recurrent model for learning very long time dependencies. In International Conference on Machine Learning, pages 9168-9178. PMLR, 2021.
13

---

[63] Zheng Liu, James P Gaska, Lowell D Jacobson, and Daniel A Pollen. Interneuronal interaction between members of quadrature phase and anti-phase pairs in the cat's visual cortex. Vision research,
1992.
[64] Hassan K Khalil, Control of nonlinear systems. Prentice Hall, New York, NY, 2002.
[65] John R Silvester. Determinants of block matrices. The Mathematical Gazette, 84(501):460-444 2000.
[66] Peter Lancaster.
and vibrating systems. Courier Corporation, 2002.
[67] Françoise Tisseur and Karl Meerbergen. The quadratic eigenvalue problem. SIAM review, 43(2):235-28) 2001.
[68] Peter Lancaster. Stability of linear gyroscopic systems: review. Linear Algebra and its Applications, 43963)6686-766, 2013.
[69] Oleg Kirillov. Nonconservative stability problems of modern physics, volume 14. Walter de Gruyter GmbH & Co KG, 2021.
[70] Wolfhard Kliem and Christian Pommer. Indefinite damping in mechanical systems and gyro- scopic stabilization. Zeitschrift für angewandte Mathematik und Physik, 60:785-7955 2009.
[71] Abraham Berman and Robert Plemmons. Nonnegative matrices in the mathematical sciences.
SIAM, 1994.
[72] Li Deng. The mnist database of handwritten digit images for machine learning research [best of the web]. IEEE signal processing magazine, 29(6):141-14 2012.
[73] Shaojie Bai, Zico Kolter, and Vladlen Koltun. Deep equilibrium models. Advances in neural information processing systems, 32, 2019.
[74] Quoc V Le, Navdeep Jaitly, and Geoffrey E Hinton. A simple way to initialize recurrent networks of rectified linear units. arXiv preprint /5ff:://00.00944/ 2015.
[75] Martin Arjovsky, Amar Shah, and Yoshua Bengio. Unitary evolution recurrent neural networks.
In International conference on machine learning, pages 1120-1128. PMLR, 2016.
[76] Frederik Beuth and Fred H Hamker. A mechanistic cortical microcircuit of attention for amplification, normalization and suppression. Vision research, 16:2411257, 2015.
[77] Rachel N Denison, Marisa Carrasco, and David Heeger. A dynamic normalization model of temporal attention. Nature Human Behaviour,
2021.
[78] Katrin Hermaann, David Heeger, and Marisa Carrasco. Feature-based attention enhances performance by increasing response gain. Vision research, 74:10-20, 2012.
[79] Katrin Herrmann, Leila Kouhsari,, Marisa Carrasco, and David Heeger. When size matters: attention affects performance by contrast or response gain. Nature neuroscience, 13(12):1554-155) 1559, 2010.
[80] John HR Maunsell. Neuronal mechanisms of visual attention. Annual review of vision science, 1:373-391, 2015.
[81] Amy M Ni and John HR Maunsell. Spatially tuned normalization explains attention modulation variance within neurons. Journal of neurophysiology. 1 118(3):1933-191) 2017.
[82] Amy M Ni and John HR Maunsell. Neuronal effects of spatial and feature attention differ due to normalization. Journal of Neuroscience, (9(2)):5493-55 -5005, 2019.
[83] Philipp Schwedhelm, B Suresh Krishna, and Stefan Treue. An extended normalization model of attention accounts for feature-based attentional enhancement of both response and coherence gain. PLoS computational biology, 12(12):e 1005225, 2016.
[84] Philip L Smith, David K Sewell, and Simon D Lilburn. From shunting inhibition to dynamic normalization: Attentional selection and decision-making brief visual displays. Vision Research, 1 16:2119240, 2015.
[85] Xilin Zhang, Shruti Japee, Zaid Safiullah, Nicole Mlynaryk, and Leslie G Ungerleider. A normalization framework for emotional attention. PLoS biology, 14(1) 10025788 2016.
[86] Matteo Carandini. Amplification of trial-to-trial response variability by neurons in visual cortex. PLoS biology, 299),2264, 2004.
14

---

[87] Jeffrey S Anderson, Ilan Lampl, Deda C Gillespie, and David Ferster. The contribution of noise contrast invariance of orientation tuning visual cortex. Science, (999005498)1168-11999 2000.
[88] Nicholas Priebe and David Ferster. Inhibition, spike threshold, and stimulus selectivity in primary visual cortex. Neuron, 47(4):482-497, 2008.
[89] Dennis S Bernstein. Matrix mathematics: theory, facts, and formulas. Princeton university press, 2009.
[90] Fuzhen Zhang. The Schur complement and applications, volume 4. Springer Science Business Media, 2006.
[91] Roger A Horn and Charles R Johnson. Matrix analysis. Cambridge university press, 2012.
[92] Adam Paszke, Sam Gross, Soumith Chintala, Gregory Chanan, Edward Yang. Zachary De Vito, Zeming Lin, Alban Desmaison, Luca Antiga, and Adam Lerer. Automatic differentiation in pytorch. 2017.
[93] Diederik P Kingma and Jimmy Ba. Adam: A method for stochastic optimization. arXiv preprint arXiv: 141.. 69800 2014.
[94] Kaiming He, Xiangyu Zhang, Shaoqing Ren, and Jian Sun. Delving deep into rectifiers: Surpassing human-level performance on imagenet classification. In Proceedings of the IEEE international conference on computer vision, pages 1026-1034, 2015.
15

---

A Derivation of ORGaNICs
Here, we derive a generalized 2-neuron types (excitatory and inhibitory) ORGaNICs model for a high-dimensional input. The system presented in Eq. Di special case of this generalized model where p 2 and a = v[al. Assuming W is the normalization weight matrix, and z is the input drive, we can write the normalization equations for principal neurons with complementary receptive fields as,
P σ² (L p) Y. o² W ([z)p+
z
(18)
Note that typically the exponent of the input p 2 for cortical neurons. |z| I and |-z| repre- sent the contribution of neurons with complementary receptive fields to the normalization pool. Mathematically, we have, |z] + |-z|" |z|".
Here, we derive, for general p. the dynamical equations that have the fixed point defined by the normalization equation above. First, it is important to distinguish between the membrane potentials and the firing rates of neurons. The coarse (low-pass filtered) membrane potential of a given type of neuron is denoted by the vector, y, with the corresponding firing rates of y (y"))att. The instantaneous firing rates of the neurons are obtained from the corresponding coarse membrane potentials by applying rectification, denoted by [.J. and power law (sub/supra-linear) activation for different types of neurons (8). Therefore, for a set membrane potentials x the instantaneous firing rates are x [x| Specifically for principal neurons, have, yt [y and I-yl" P. Combining the firing of principal neurons with the complementary receptive fields, y and y Eq. 18 can be alternatively written as,
yy y ty.
(19)
Oº "
Now for the principal neuron y receiving an input drive 2j, can rewrite the normalization equation for each neuron as,
12
lngl" of
(20)
In the ORGaNICs paradigm [47). the steady-state activity of the principal neurons (single equation for complementary receptive fields) is a weighted sum of input drive and recurrent drive,
dy Tys dt
(1 az) Σ Wrik (lvt) (yE (1p)
byzj
(21)
Weighted inpot drive
Weighted recurrent drive
Here w, the weights of the recurrent weight matrix W, encoding the recurrent/lateral connections between the principal neurons y; y is the firing rate of the principal neuron k, given by y [yv]?, and y/ = 1-wke P is the firing rate of the complementary principal neuron k. b, is the gain to the input drive 2j which simulates attention (realized via gain modulation). (1 - -af) is the gain to the recurrent drive which controls the recurrent amplification.
Now, we find the dynamics of the inhibitory neurons aj with firing rates at that yield stable dynamics with the fixed point given by Eq. 18 First, we assume that the recurrent weight matrix W, I. Also, note that the following identity holds: (Mit) ) (vp 1/e Lynl L-yk] = yk- Therefore,
21 can be simplified to,
dy; リı byzj + (1 - (1-a)) yi
Ty dt
(22)
At steady-state, the fixed-points (15, az), satisfy the following relationship.
az y3 = b,ا Taking modulus and raising both sides to power p. we get,
(23)
|a uzll P = 1b,
(24)
16

---

or in vector form,
a
o yy P b x
(25)
(at)* o ly ④ x
From Eq. 19 we know that
|z|" = o y.ll + y O [z|")
(26)
Since we can write the element-wise product between two vectors x and X2. X1 ② X2 X2 ④ X1 D(x1)x) Where D(x1) is a diagonal matrix with elements of the vector X1 on the diagonal. Therefore, using this fact we can rewrite the equation above as,
kz| = 8 @ y P + (ly& (W (W
1z|" = o o ys + (ly W |z|"
(27)
D (1y W| |z|" = o" O ys
1z|2= [ D (ly, W ( + l1
Substituting the expression for |z|" into Eq. 25 we get,
yel" b" I D (lya PP) W|
10 1lyl'" ly)l
(28)
We simplify this equation as,
(at) 2 lyel"
I D (l) W|
0 3 yya
bP
(a @ lyal"
(1/ly,, D (ly. w]
(1/ly& ( (o'
(29)
b²
(a 2
(a &
lys
w
bp
b²
Now we assume that all of the entries of b are equal to constant by. therefore, we have,
( (at)²
w ((at)* @ P)
lof"""
(30)
(at)² 60 + W ((at)"oly.l"
Element-wise multiplying both sides by a ) on the left, we get,
0 =
( 1oo" + w ((at)* о ly,l")]
(31)
This equation is true at the steady-state, but it is also in form similar to that of such that we have a weighted input drive and weighted recurrent drive, therefore, we propose the following dynamical equation for a.
180 W o
a
(32)
(at)"
This equation naturally follows Eq, 3T at steady-state and thus also follows the normalization equation
Eq. 19. Note that the equation above is true for any choice of at, but cortical neurons have been experimentally observed to have an exponent close to = 2, i.e., y = [y] and y = [-y] 2 Therefore, we get a particularly simple form of the equations when = 2 and a case, the equation for given by,
|a] In this
To oaa - [800² + W (a (a))²
(33)
We reintroduce the recurrent weight matrix W, and to simulate the effect of attention by gain modulation, we define different gains for y and a neurons. Additionally, replacing y² + yt +y and y + v-+ Vy- yields the dynamical system that we analyze,
oy (1-a+) (w, (Vyt-Vs
(34)
Ta à - by W ((y+ +y 3 @a+2)
17

---

Finally, in the most general form, i.e., any choice of and nonlinearity for a the dynamical system of ORGaNICs is given by,
at) 1 (., (1xt) (5-)
Ty my
(35)
[bj O o² W ((y +y-) (at)*)
Ta
(a+)
B
Stability theorem
Theorem B.1. For system of linear differential equations with constant coefficients of the form,
IX + Bx + Kx
(36)
where B _xxx and K Ruxx a positive diagonal matrix (hence K 0), the dynamical system is globally asymptotically stable if B is Lyapunov diagonally stable.
Proof. The stability of the system is defined by solving the following associated quadratic eigenvalue problem,
((A) det BX B K)
(37)
The spectrum of ((). i... {) C det(L(()) 0} also known as the eigenvalues of the system. The system defined by Eq. is globally asymptotically stable if all of the eigenvalues have negative real parts. We take the direct Lyapunov approach to prove the stability of the linear system. We write Eq. 36 the matrix form as Jz, where,
z
and J K
-B
(38)
We will first prove the Lyapunov stability (Re(As) 0) of this system find the appropriate block diagonal matrices and then we will prove global asymptotic stability (Re(Xa) < 0). To prove Lyapunov stability (Appendix B.2), we propose a Lyapunov function Viz) = z Pz, where P is block positive definite matrix defined as follows,
A
0 T
P =
(39)
where T E Ruxn is positive diagonal matrix such that TB + BTT O (notation for positive definite matrix) and A Roxx flexible symmetric positive definite matrix that we will find using the second Lyapunov criteria. Note that such matrix T exits since B is defined to be Lyapunov diagonally stable. It can easily seen that > using the first criteria Scction B.I since A and T 0. Additionally, since T > 0, is invertible.
Now, for Lyapunov stability, we need (z)) T (PJ+ + JPJ+JPP)2 < 0. Therefore, we find A such that Q = (PJ + J P) is positive semi-defnnit.
Q (PJ
P)
I
0
K
0
T
K
B
I
B
0
T
(40)
KT
TK A
TB
T
We want to define A > 0, such that Q 0. Using the second criteria from Section B.I
need
TB- B T > 0 and - ATB + T) "(TKK - A) > 0. The first condition is satisfied by the definition of T. For the second condition to be satisfied, Note that both and K are positive definite and diagonal, therefore they commute (KT TK) and A is symmetric and positive definite. When A TK, the LHS of the second condition becomes 0 z 0. Therefore, the system is Lyapumov stable.
obvious candidate for A is TK.
Now, we prove the global asymptotic stability of the system by again using the direct Lyapunov approach. We propose the Lyapunov function the same form as before, i.e., V(z) = z T Pz, where P is positive definite matrix. But for global asymptotic stability, we need a more stringent condition
18

---

on the Lyapunov function, V(z) z (PJ + J J'P) z < 0. Drawing inspiration from the previous exercise, we consider the following form of the matrix P.
TK
d
P
(41)
d
T
Here, <> O is scalar whose magnitude is to be determined based on the Lyapunov criteria for asymptotic stability. First, we need P 0. Applying the first criteria from Section B.T we want to satisfy, TK T > Second, for (z) < O, we want Q = - - (PJ+JPP to be positive definite. Q is given by,
TK
c
0
I
0
K
TK d
Q =
d
K
B
d T
(42)
2eK
cB
cB
TB
+
B
T
2c1
Again, apply the first criteria from Section B.1
Q
if and only if TB B 2c1 > 0 and
2cK B (TB + T B positive definite matrix M. Therefore, we have prove that there exists 0 which satisfies the following conditions,
2d) -a -d1)= B O. To simplify notation we replace TB + T with
TK - 22 0
M 2cI > 0
2K cB (M-2el
B >
Assuming t, and k; to be the diagonal values of the positive diagonal matrices T and K, we get the following two conditions, < min ((v),)) and a/2, where a is the smallest eigenvalue of M, or min (AM).
Now, for the third condition, will use number of facts about positive definite matrices which are all listed in [9). We first consider the following matrix inequality, M ž al. This notation is equivalent to saying that M aT is positive semi-definite M aI 0. Therefore, we have,
M 2dI = (a Ic)I
(43)
Assuming, is small enough such that matrices on LHS and RHS are positive definite, we have,
2c Iz -M -M-
(44)
Since B is nonsingular, it is full rank, therefore, we have,
BBT 2 BM (2M-) 221) TB
(45)
a 2c
Multiplying both sides by e, we get,
BB
7 cB (M- 2e1) B
(46)
2c
Notice that BBT is positive definite matrix. Let 3 be the maximum eigenvalue of BB or max (ABBT). Therefore, we can add an upper-bound matrix to the inequality as follows,
cB I
BB
cB (M 2c1) -B
B
(47)
2c
a
Now, we find such that,
cB
2K
(48)
a 2
The range of values for which the above inequality is true is,
2kza
< min
(49)
3 + 4k;
19

---

If e satisfies the inequality above, we have,
cB
2K
BB
- cB (M
2el)" B
(50)
a
a- 2c
Therefore, have, 2K cB (M-
211)-B B and the third condition required for is satisfied.
This implies that there exists a range of 6, given by,
2k:a B + 4h:
min min (liv Ki); 2'
min
(51)
for which
TK d el
P
(52)
T
a valid Lyapunov function for asymptotic stability. Therefore, the dynamical system is globally asymptotically stable.
B.1 Positive definite block matrices (Schur complement)
For a symmetric block matrix of the form,
A
B
P
(53)
B
T
C
with A = A T and = CT, If C invertible the following two properties hold, (90).
* P > 0 if and only if C 0 and A BC-'B" > 0.
IfC > 0, then P ž if and only if A - BC- 2 0.
B.2 Lyapunov stability criteria
Consider non-linear autonomous dynamical system defined as x f(x), with point of equilibrium at x = 0. Where x E D C R" the system state vector and f(x) : - RR" is a continuous vector field on D (contains origin). The dynamical system is called Lyapunov stable if there exists real scalar function V(x)= R + IR, also known as the Lyapunov function, such that it satisfies the following conditions,
V(x) = 0, if and only if = O.
V(x) > if and only if x بf
(x) < O, vx + O. Note that for asymptotic stability, we require the strict incquality
(xx
< 0.
For linear dynamical system of the form x Jx, where Ruxn, with a point of equilibrium
x 0. Consider Lyapunov function Vxx) of the form x Px, such that P > 0. By the definition of a positive definite matrix, it satisfies the first two conditions, namely, V(x) = x Px = when
x 0 and Vx) = x Px > 0 when æ + O. For the third condition, consider V(x).
d
vxx
V(x)
dt
d
Px
(54)
= T + x Px
=x (PJ + PP
For stability, We need V(x) x J x < 0. This is satisfied when the matrix PJ JTP O. summary, linear dynamical system Jx is Lyapunov stable if there exists positive definite matrix P. such that PJ + J P is negative semi-definite.
20

---

C Analytical eigenvalue for the fully normalized circuit
Here we show that when all of the normalization weights in the system are equal, to value a, and the various parameters are scalars, i.c., Ty = Ty.1, Ta = 1. by = byl and a = o1, we can derive closed-form analytical expression for all the eigenvalues. Considering these assumptions, we can break the determinant in Eq. g into diagonal and non non-diagonal part as follows,
D(va.)
WD (y. D (vaa) TyTa
det((J X1) det A Ta
+
X
(55)
Ty
Ta
Consider the non-diagonal part of the matrix in the determinant, (A/Ta) WD (y²). Since W is matrix with all entries equal to positive constant, it is rank 1. Therefore, it can be written as the following outer product.
W
1]
(56)
Therefore, the non-diagonal part of the matrix can be written as,
b²
z² (b2 x2²)
Aa Ta
WD
uvv
(57)
0 02 1
+
where u = [1, 1, and v (b² 3 x²) nant lemma which states that,
1b1 + W (b² o z²)).
We use the matrix determi-
det A yuv ) = (1 y1 A 'u) det(A)
(58)
The matrix A is given by,
D ag
D (va.)
21
Ta
D
2 1 + W (b² o x
D 2Vo22 + W @b²
I
= X X
Ta
Ty
Ta
(59)
D 2 + a b o 2/13)
D Vo + b C 2/
I
-
Ta
To
1
2
b
bb + b z
1
+
A
Ta
Ty
Here, ||x]|| represents the Euclidean norm of Therefore, A = 5I, where d is a quadratic scalar polynomial in Now using Eq. 58, we can write
λa
det Al)
8"
Ta
(60)
X
b z 2
8n-1
Ta
2bg
+ +all b о 2
Solving for the eigenvalues, det(J AI) 0. get 2(n 1) repeated solutions by equation gn-1 0, where each solution is found by solving = 0,
1
1 b 2
2b b z
x² +
(61)
Ta
This gives the following strictly negative eigenvalues,
A Ta
& ==
(62)
T
21

---

The potentially complex eigenvalues are given by solving for zeroes of the first part of the factorized determinant,
b (o + b
0
0
2
+ 2z 0 (63)
x +
TV
Ty
The eigenvalues found by solving this equation have negative real parts (as expected) for all choices of parameters since the coefficient of À and the constant term of the quadratic equation are positive for all choices of parameters and inputs. Therefore, this solution does admit complex roots (thereby oscillations) for certain choices of the parameters which can be found by solving this quadratic equation.
D Convergent theorem
Theorem D.I. A matrix A of the form,
A D(t(
(64)
v Wu
convergent, i.e., its spectral radius is less than one, W R and u, R with the additional constraints < 1; < 1, u, 2 0, oi > 0 and Wii 2 0 for all i.
Proof. Assuming the constraints mentioned in the theorem, first notice that,
u
WD
(65)
Wu
a nonnegative matrix. Further, because D(t) is positive diagonal matrix with all entries less than 1, we notice that the following is true element-wise,
D(t) W
< WD
(66)
v Wu
+ Wu
We denote the spectral radius of A, max À o (A)), by p(A), where a (A) is the spectrum of A and make use of the following inequality for the spectral radius of nonnegative matrices.,
Theorem D.2. (Theorem 8.1.18 from (971) Let X and Y be nonnegative matrices with spectral radius p(X) and p(Y), respectively. If X Y (Iij I Wij, V i,j), then p(X) S p(Y).
Therefore, have,
p (AA) sp
WD
(67)
v + Wu
Now the elements of the matrix are given by,
W11
W12
-, wyyay
WD
W21
W22
0/ +,, way
(68)
v Wu
Let Sty be the j element of this matrix. Upon multiplication of the matrices above, we find that,
Wifdy
8ij
(69)
+ De W3
Now use the following theorem that provides upper bound for the spectral radius based on the entries of the matrix.
Theorem D.3. (Theorem 8.1.26 from (91)) Let S be nonnegative matrix. Then for any positive
vector X E ga with entries Ij, we have,
PSS) < max
Ti
(70)
1s Fi
j=1
22

---

We pick 2j = U/ - Σx u which is a positive vector, and apply the theorem to the matrix. Substituting Sij and 2.j. we get the following bound for the spectral radius,
V We
max IsiSh + Ex WiRlk Uk
Wayug
+
Σ
Wiklk
(71)
2 + Σx W
Upon simplification, get,
V Wu )
E, Wijll
Σx wA
p WD
S
max ISSS Uj + Ex W
max 1S
(72)
U2 + Ex W.
The inequality is true, because, for all i, the denominator is larger than the numerator because of the extra positive term 14. Therefore,
V Wun )
)
(AA) s
WD
<1
(73)
and A is convergent matrix.
Linear stability analysis of the twoddimensional model E
Here consider the dynamical stability of the 2D model containing one neuron each of and when the recurrent scalar, w, can take any positive value. The dynamical system to be analyzed is given by,
TuU = yy bz + 11v[aa wry
(74)
Ta a + bi0² + wy² a
with positive real constraint the following set of parameters, Tyr Tay b, bus a, w. We first notice the symmetry in the dynamical system about y = O. Replacing, ÿ + -y, we get the mirrored dynamical system,
9 = - 1(-2) + (1- [a w-ÿ
(75)
Q = a + 0l0² + wy² [a]
These equations are the same as Eq. 74 up to a sign change in z. Therefore, we can derive analogous conditions for stability for < 0 once the conditions for > are known.
The steady-state of Eq. 74 (yss a.) satisfies the following equations,
ys = bz + (1-v V(a.l) wrye
(76)
0. -0,00 2 + wyy [a.]
(77)
Since the RHS the second equation is always positive, if a root exists, we have a, >> 0. Therefore, remove the rectification around a, and look for positive solutions for a.. Upon rearranging the terms, we get,
(1 w, + Wr / /as)as)))))))))))) V = ba
(78)
(1 - wy a. = 0/
(79)
Substituting ys from Eq. 78 in Eq. 79 we get the following quartic equation in m = Vas.
2 2(1 w. w., 3 ((1 w wb² 2² 2/ 2) m2
(80)
2(1 w. ).,,), b302 2 (1 w, ) 0
For a valid fixed point (ye, ag). i..., y/ E and a2 E R a necessary condition that a, must be positive and real, which in turn implies m must be positive and real.
Theorem E.1. A fixed point, (y., a.), is valid if and only if it satisfies va. 0.
23

---

Proof. Due to Lemma E.2 we have: fixed point, (ya, az), is valid if and only if it satisfies m > boo.. Further due Lemma E.3 we have that the conditions m > 0 and m boo are equivalent. Combining these two we get the statement of the theorem.
Lemma E.2. fixed point, (ys, (s), is valid and only if satisfies vu, >> bno.
Given, m = va > byo, we have as > byoz, therefore, Og E R and ye Vm²² 02 / (m) therefore, 1/s E P. = Given as R and y/ E R, implies Rt Now, Eq. 79 posits that, m byo / wy? therefore, m byo.
Proof.
Lemma E.3. For a positive real root, m = Vug to the quartic equation (Eq, 80 the condition is equivalent the condition >> byo. Further, fixed point satisfies 0 < buo.
Proof. Given byd, have O because by > 0 and > 0.
Since is root of given wy 1, > boo. Therefore, > 0 is equivalent to m >> boo.
the quartic, it also satisfies Eq. 79 Therefore, have buo / / wy Since
m > 0, this implies that
This also implies that there exists no fixed point that satisfies < m byo.
Further, the Jacobian matrix at the fixed point of the dynamical system, in terms of the parameters and the fixed point is given by,
2r, wh Vaa
Ty
J
(81)
2wa Us
1+wp?
From the linear stability theory, we know that fixed point (Yar a.) asymptotically stable when the real part of the cigenvalues of J are less than 0, i.e. Re(Ar) < 0. For 2D system, this is equivalent the conditions: Tr(J) and det(J) 0. The trace of the Jacobian matrix is given by,
Tr(J)= W,
Wr
a2
wy
Ty
Wr
+
Wy V02
1 wwy
(82)
Ty
Wy
+
Wry
/
02
bzo²
Ty
ag
The determinant of the Jacobian matrix is given by,
Wr - We /ry
1 + wy?
WyHss
2waayys
dJt(J)
+
Ty
Ta
2 VaanTyy
Ta
(1 w, w,)(1 wy?)
wr w a
(83)
TaTy
TaTy
(1 wrr
+ wr wr
TaTa
0g
Now we consider the different conditions on and the input drive = for stability and state the various cases as theorems along with their proofs,
E.1
Contractive constraint on recurrence (0 < WT S 1)
2 >>00 There exists unique fixed point with ys > 0 and a2 > O and it is asymptotically stable.
Proof. The existence and uniqueness of the fixed point, along with its stability properties, are established by Lemma E.4 Additionally, the positivity of a, (i.e., 0) ensures that the expression wy + Wy 0 O. Consequently, given that 0, follows from Eq. 78
that y, > 0.

---

2 < 0:: : There exists unique fixed point with y < 0 and o₂ > 0 and it is asymptotically stable.
Proof. Since this condition becomes equivalent to = > 0, up to a sign change in y (Eq. 75 it is straightforward see why this is true.
Lemma E.4. Given 0 wr and I, there exists unique fixed point and it is asymptotically stable.
Proof. We observe that given 0 < w, S 1, if there exists valid fixed point, which satisfies ag Tr(J) in Eq. [2 is less than (since Wr + wr > 0 for any ya, > 0), and dct(J) in Eq. 83 greater than for all combinations of parameters and z R. Therefore, constraints on the parameters that allow for at least one fixed point. The fixed point, m = v satisfies the quartic polynomial in Eq 80 Due to Theorem E. for valid fixed point, we need find positive real roots that satisfy m > 0.
need to find the
The sequence of signs of the coefficients of the polynomial when 0 S wr < 1 is given by (,.+ -.. We use Descartes' Rule of Signs, which states the following: The number of positive real roots of a polynomial p(x) is either equal to the number of sign changes (omitting the zero coefficients) between consecutive non-zero coefficients of p(z)., or is less than this by a multiple of 2. Since there is exactly one sign change from left to right in the sequence, regardless of the sign of the coefficient of m² we know that the equation above has exactly one real positive root for m. This further implies that a, has exactly one positive root and the corresponding fixed point (Me, ),), is locally dynamically stable for all of the combinations of parameters and 2 E IR.
Note that the result also holds for w, = 1. There is exactly one positive fixed point (ya, a,) and it is given by, i.e. y. bz/ /e/Dbooz + wt? z2 and a, lifo + wb2=2.
E.2 Expansive constraint recurrence (w, 1)
Now we consider the case when w, > 1. We first find an alternative form of the determinant of the Jacobian 83) that is more amenable to the case, W, > 1.
1 ((1 wr)(1 wy²) + Wr wrrvas)
det(J)
Ta
(84)
(1 Wr + Wr ag + (wr - 1) 1)wy?)
Ta
Rewriting the trace,
Tr(J) (1
- Wir + Wr
02002 Og To
(85)
Ty
The properties of stability can be summarized in the following two cases,
= >>> 0: There exists exactly one fixed point with ys > and a, > 0 and it is asymptotically stable. Further, if byo > 1 - 1/W,, there exist no additional fixed points. But if byo < - 1/w,, then there exist either two fixed points with ya < 0 and ag 0 and they
may or may not be stable.
Proof. Since the fixed point satisfies Eq. 78 there two possibilities, either 1 1- W.-
Wy 0a > 0 or Wy + Wy / 0a, < O. We consider them separately,
1-w, w w. aa > 0: : If= > 0, we have ys > 0 because ys satisfies Eq. 78 Further due
to Lemma B.5 this fixed point is unique and asymptotically stable for all combinations of parameters.
1- Wr + Wr 0 If > 0, we have ya < 0 because y, satisfics Eq. 78 Further, are given buo > 1 - 1/w,, no root exists. We prove this by contradiction, assuming
a root exists that satisfies ya,- < 1- 1/ w,,, since byo > - 1/w,, this root satisfies
0 < Va,- <
byo, but due to Lemma E.3 such root exists, therefore we have
contradiction.
Now if we have bog < 1 - 1/w,, the root for m = Va, must satisfy bng < m < 1 1/w,, Since either one or three roots satisfy m 0 (equivalently m > boo) in
25

---

Eq. 80 and exactly one root satisfies m > 1 1/20, (Lemma E.5) we conclude that the number of roots for in the interval (by0r, 1/w,) are either two or none.
< < There exists exactly one fixed point with ys < O and > 0 and it is asymptotically stable. Further, if boo > -1-1/ 1/w,, there exist no additional fixed points. But if boo < 1/w,, then there exist either two no fixed points with ys > 0 and 02 >> 0 and they may or may not stable.
Proof. Since this condition becomes equivalent to z >> O, up to sign change in y (Eq. 75 it is straightforward see why this is true.
Lemma E.5. Given w, > I. there exists unique fixed point, (ye, az), that satisfies 1 1- + w Va O and it is asymptotically stable.
Proof. The fixed point, m = Vag. satisfies the quartic polynomial in Eq 80 for valid fixed point, need to find positive real roots that satisfy >> U.
Due to Theorem E.
The sequence of signs of the coefficients of the polynomial is given by (+,,, +, -). Regardless of the sign of the coefficient of m², using Descartes' Rule of Signs, we find that since there are 3 sign changes, there are either 1- W, / > O, these roots are valid only if they satisfy 1 w, + w.m > O, We make a transformation m + f + (1 l/w,) which gives us,
positive real roots for m. Since, we are given w, > and
> 1 1- 1/wr-
w 2 2(1 - w, ww, m 3 + ((1 - w,)2 wb²z 22 61j020. w. 2) m² 2wb² 2²(1 1/ w/) )m wb²²² 2² (1 1/ w,, )= =0
(86)
The number of valid roots is given by the number positive real roots of this polynomial. The signs of the coefficients given by, (-, =, - -). There is exactly one sign change, hence there is unique fixed point that satisfies 1 w, wr > O.
Now we prove that this fixed point is asymptotically stable. It can be easily seen that Tr(J) < in Eq. (85 and (Jt(J) > Eq. 84 when 1 wy + wryaa, > and wr, > 1. Therefore, this unique fixed point is asymptotically stable.
L
F Analysis of the Rectified model
Here, we present the stability results for the model with only the positive part of the complementary receptive fields present. The model is given by the following dynamical equations,
Ty + - b + (1 a*) [W,y]
(87)
Ta ④ á
bgg の 0² W
@
We will follow the same procedure for stability analysis as we did for the main model and state the key steps in the various derivations for stability.
F1 Stability of the high-dimensional system
We analyze the stability of the system when W, = I. Upon substituting y+ [y| and a -
a., we get the following dynamical system,
aa Lyl
(88)
To 2 à - - 22
(Ly] [a))
The fixed point is given by.
b G Z
y. V b 0²
l-boz
W b
(89)
a, = by o 3b2
26

---

Also, [y&l is given by,
b o
ye
(90)
o
The Jacobian matrix, about this fixed point is given by,
va.0lootteee
-D D 2 WD (yJ) D (-I WD (D(Ly.J3)
28 Baya Bar
=
(91)
Here 1220 is an indicator function that returns vector the size of with locations where 2; and elsewhere, Similarly. 1z<0 returns a vector the size of z with 1 at locations where < 0 and 0 elsewhere. Further J AI is given by,
A11
A
A21
A22
D
Va,0l
AI
D
yy
(92)
Ty
200 (ay
D
(a @ ly.))
W
D Ta
(-1
W (Lv, AV
Since A₁ and A12 commute, can write, det([ simplification, we get,
AI) det det(A22A11 12)
Upon
1 + X D 1
Va. ④ 1220 + 1250
det(J - AI) =det
+D
Ta
Ty
(93)
D
WD (Lyx13)
+D
aas
O 1x20 + 1xc0
Ty
④ Ta
Note that, here we used the fact that y. has the same sign system is stable
seen from Eq.
89 The dynamical
eigenvalues this characteristic polynomial have negative real parts. Just like
main model, map this quadratic eigenvalue problem of the form [(X) det X Now, we see if the conditions of Theorem are met. The stiffness matrix is given by,
Va 2 1220 12250 =
b o
W b の 2/2 の 1220 + 1xc0
K K-
(94)
TTO
3 Ta
Clearly, K a positive diagonal matrix for choices of parameters and input. Now we check the B is Lyapunov diagonally stable matrix, which is implied when B admits regular convergent splitring, i... has a representation of the form = M - N, where M and have all nonnegative
entries and M N has spectral radius smaller than 1.
Va, C 1220 1z50
B D Ta
D
D
WD (lv.12)
(95)
Ty
Ta
Since M is positive diagonal matrix, all the entries of M are nonnegative. Also, since all the matrices involved in the definition of N are nonnegative, therefore their product, N, is nonnegative. We are only left to prove that M -NN-- S has spectral radius smaller than 1. Consider the matrix S.
==
WD
(96)
1+ (Ta/T)) 2 (va, 2 1220 + l2co)
by 00 W Wboz
Theorem puts a bound of on the spectral radius of this matrix. Define t 1/ (1 (Ta/TT)TT)) (va, ② 1220 l2<u)). + bb 0z z/2 and v + by 2 0² we notice that they follow the constraints of the theorem, therefore, is convergent. This implies that has convergent regular splitting, therefore, the linearized dynamical system is unconditionally globally asymptotically stable (and nonlinear dynamical system is locally asymptotically stable) across all the values of parameters and inputs.
27

---

F.2 Analytical eigenvalue for fully normalized circuit
Following procedure similar to that in Appendix C when all of the normalization weights in the system are equal, to value and the various parameters are scalars, i.c. Ty = Ty Ta Tal. bo by 1 and o1, we can write the analytical expressions for eigenvalues. The characteristic polynomial is given by,
Aa D Ta
81 1220 + 821 )
det(J - A1) =
81'07"
(97)
where, u = [1, 1J" v= = b o / (o² W 2|²); nı and n2 are the number of nonnegative and negative values, respectively, in the input drive ô1 and d2 are given by,
81 = X2 +1
2
b
z
a z
Tv
Ta
(98)
1
62 = X2 +\
+
Ta
Ty
Ta
Simplification of Eq. 97 gives us,
all[b : 2
det - AI) = 81 Ta by + + 2
61
(99)
Since the characteristic polynomial is product of quadratic polynomials, can solve them analytically. The strictly negative eigenvalues are given by,
& A= Vo + + 2 2
A A=
(100)
Ty
Ty
The potentially complex eigenvalues are given by the solution to the following quadratic equation,
Ta all a bo 22)22 + Vo + a 2 12
o²t 2²bg
+ a b o 2 1
+
+
0
Ty
Ty
(101)
F33 Linear stability analysis of the two-dimensional model
The dynamical system to consider is,
Tyú y ba (I v[al] [w,yl
(102)
T a + 1100² + w Lu! LL]
Since valid fixed point must have a, 0, the fixed point (y/. az) satisfies,
(1 - Wp + wr wz Lyal - [-yal = bz
(103)
(1 w Uu = b200
(104)
We divide this into two cases,
Ya > 0: The fixed point is given by the equations,
(1 wr + wy va,) Us bz
(105)
wy/²) 1/02
(106)
A fixed point, (ya, as), is valid only if it satisfies yo E IR and a2 E IR.. The Jacobian matrix is given by,
W. waya.
Wr
Ty
2/a. Ty
(107)
2100 .No
14 Ta
Ta
Note that this is equivalent to the main model, with the additional constraint of ys > whose stability analysis is presented in Appendix E
28

---

< : Eq. 103 & 104 yield us unique fixed point ys = bz and a, = =8/02. Since ys < 0 The Jacobian matrix about (ys. a,) given
and y, 62, this is only possible when by,
1 4y
(108)
Since the eigenvalues of J are X -1 1/TT -1/Ta are real and both negative, this fixed point always stable.
Combining the cases above and results already established in Appendix characterize the stability of the system as follows,
E3311
Contractive constraint recurrence (0 wr 1)
2 There exists unique fixed point with ys and o, > 0 and is asymptotically stable
z < <0: There exists a unique fixed point with ys < 0 and ag > 0, given by, (bz, b8 and
it asymptoticallyy stable.
F33.2 Expansive constraint on recurrence (w, 1)
0: There exists unique fixed point with ys and 0 > 0 and it is asymptotically
stable.
2 <00- There exists exactly one fixed point with y/ < and u, given by, (bz, biga"). and it is usymplotically stuble. Further, if boo 1- 1/Wr, there exist additional fixed points. But if byo < 11 1/Wr, then there exist either two or no fixed points with ys and 0g 0 and they may may not be stable.
G Iterative algorithm
In this section, we present an iterative approach to finding the fixed point for ORGaNICs with an arbitrary recurrent weight matrix. We show that this algorithm converges few steps (2-10) with great accuracy. We consider the system given by Eq. 15
Ty - v (1 a @ (W,y)
(109)
④ á a byg o² + W (y²
C [a])
The fixed point of this system (y, and a.) is found by solving the following simultaneous equations,
y. = h z +
(W,y.)
(110)
a. b
W (y: aa)
(111)
These equations do not admit a closed-form analytical solution when W, I. We first find good approximation for the initialization of y, and a, and then define the iterative algorithm. The equation for can be written in terms of
y. (I+ (D I)) -)),- (b
(112)
Now applying the Woodbury matrix identity, which states that
(A + -CVV = A-1 A A "C + VA-1
(113)
the inverse in Eq. 112 with I,U I- (Va,)
I lad W,, we get,
I D (v - -1) 1 + w,) W (b (boz)
y. =
(114)
We approximate the above equation by assuming that W, is a symmetric matrix with the eigende- composition given by QAQ with Q QQ = I and A (A) is diagonal matrix with
29

---

the eigenvalues as its diagonal entries. This gives us the following approximation,
@ (D (va,)= I) 1- + A) -
Q)) QQQq)) (b
1- QD
t
Q QQQQq (boz)
va,- -
A * Va, x
Q (boz)
1 QD
A X va.
(115)
1
100)) (b 22)
-- Q
I D-D
X X
a
1
QD
Q ((boz)
y
1
QD
AQ (boz)
A X
va.
We approximate the cigenvalues, X, by the maximum eigenvalue of the W, and assume the entries of ag are identical. This gives us the following initial guess for ys.
1
y: D
QAQ (boz)
Am
X2 m
+
12 v2
a
(116)
W, (b ④ z)
X2 + X2
For the initial guess of a.. we use Eq. the corresponding y found by using Eq. 116 This gives us the following,
111 and plug in the following on the RHS a, - by o o² and
W, Ob
a. = o² ob + w
Obi
(117)
Am - Ag, + X2
Next we update the y. and a, by performing the following iterations derived using Eq. For instance, (y' a are given by,
110
y W, + D a big
a W.)
(bOz)
(118)
W (y!) a?
This procedure is summarized in Algorithm Substituting Am = 1 in Eq. 116 & 117 yields simpler initial conditions and gives us Algorithm Even though we had assumed that W, should be symmetric, in practice we find that this algorithm leads to fast convergence even for non-symmetric matrices. The fast convergence owed to the fact that we have a good initial approximation of the solution. We also find that this iteration scheme works only for recurrent weight matrices with maximum singular value I.
H Energy of ORGaNICs
Here, we find the energy (Lyapunov function) that is minimized by the dynamics of the ORGaNICs in the vicinity of the normalization fixed point. We consider the dynamical system with W, I, which is given by Eq. Upon linearizing about the fixed point get the following linear dynamical system,
D
/aa
-D
Y
20 20
OTOOT
y y, a
(119)
à
D
W WD (a, Oy.)
I
WD (y?))
30

---

Algorithm 2 Iterative scheme for finding the fixed point
1: Input: ORGaNIC's parameters and input (z), Tolerance e, maximum iterations N
2: Output: Approximation to the fixed point (y,,
((...) (boz)) + bjood
3: a by +W
① (b) 3 02)
initial approximation for
4:
(b(b
// initial approximation for y
Am + Am va
5: k
6: while ||y - b oz (1- va) 3 (W,y) /|>))>>> and l < N do
7:
y -1 [W,,, +D (vā)) W.) (boz)
// update
8:
b bj @ 0² W (y² *a)
|/ a update
9:
k+
10:
end while
11: return a)
This system is dynamically equivalent (admits the same cigenvalues) to system of coupled harmonic oscillators with the following equations,
x + D
1
+
as
D
1 Ta
WD (v?) x D
a Ty 3 To
= = 0.
(120)
Ty
We can rewrite the linear system in terms of the position, x, and the velocity, v,
I
x
(121)
D
D
WD (y3)
v
TOOTa
Tv
Tu
Since this system is of the form I* given by Vzz) = z Pz, or.
Bx Kx
Eq.
36
the energy of this dynamical system is
TK c
x v
V(x, v) = x
T
=
x (TK)x + Tv T
(122)
Here, T is any positive diagonal matrix such that TB + B 0 and K is also positive diagonal matrix given by D /a./ (Ty O TT)). Now, for valid Lyapunov function, we can take e to be arbitrarily small, Eq. 51 Therefore, the energy minimized by the dynamical system is given by V(x, v) x (TK) x + Tv. This is a high-dimensional version of the energy of a damped harmonic oscillator. For single oscillator, have V(x, 1(kr² + u2) which proportional to the total energy (kinetic + potential) of the oscillator.
We now express this energy in terms of the variables relevant ORGaNICs, i.e., find V a). First, we denote the Jacobian matrices in RHS of Eq. 119 & 121 by & B, respectively. We note the simple fact that A & B related by similarity transformation change of basis). This means that there exists an invertible matrix U, such that = U BU and the corresponding transform is given by xx v] Uly y, a. IT. Assuming that U is invertible, UA BU. To solve this, we consider a block matrix representation U and find the following solution,
write this equation
D
OTOT D
(123)
This change of basis gives us the transformation,
D
y Varor
x
ys
(124)
ax
1
a
2
or,
a
O T
(y -a)
y
(125)
a
(a
a.
V
3 (y y.)
2
31

---

Substituting these expressions into V(x.v) x (TK)x + + Tv, and assuming the diagonal entries of the matrix T to be t; and substituting the diagonal entries of K, (Tsi/(Ty, Tas) get the energy in terms of y and a.
Ui agi 3/2 Tai
V a)
(y yos)
Osi Yoi
(aasi i (y/ Yos)
Yod 2 2
(a: agg)
Usa
(126)
We notice that Taylor expanding the term v Vall about yasiMst and ignoring the second order terms, we get,
Valyi Valiles + Vagi (ye Vos) + Yoi (a, aee)
(127)
2 2
Therefore the energy function, V(y, a), given by,
V =
t
Osi
Tys
aas (yi - Uoi)² + (Vasy) - vastbes)
(128)
Yes
Tas
Notice that Voillot bizi. This gives us the following expression for the energy function,
a(y,a) =
Σe
4
Ous
Tui
Usi (y:
ysi)² + (Vaiyi - b,
(129)
Uss
Tai
Further, for an ORGaNIC's model containing one and one a neuron, after removing the proportion- ality constants, the energy function is given by,
Viy, a) aa (y- un)² 2 (Vay bz)²
(130)
After plugging in the steady-state values, we get.
2
Vly, a) Ty Ta
/b802 + wb?
bz
(Vay bz)²
(131)
12122 2 + wb2
For this system, it is easy to verify that this is a valid Lyapunov function and is minimized by the dynamics of the circuit. We now demonstrate that it has the properties of Lyapunov function. First, Vyy,, a.) = 0 and Vyy. a,) > A. y + Ys & a F a,, this can be easily seen from Eq. 130 Second, we need to show that, Vyy, a) < 4 y F ys & F as- Using Eq. 119 & 131 V (y, a), we can write the total time derivative of the energy be,
dV 0V dy dt
0V da Da dl
dy dt
(132)
(2ya. 3a,ys + ays) 2
( (VasTa wyf + Ty)
20,
Since as 0 and,
(Va,To wy?T + Ty)) (VasTa + Ty 1i/02 + wb? ))
1z02
>
(133)
for all the choices of parameters and y F ys: a # as we have Vyy, a) O, therefore, it is valid Lyapunov function and can be interpreted as the energy that decreases with time via the dynamics of ORGaNICs.
I Training details
The code (written in PyTorch 922 to produce all the results can be found at http://gthub.com//mmattaaaa.aaaa lwb/dynamic- nore For both the static input and the sequential input, we train ORGaNIC' on the MNIST handwritten digit dataset [72), and to the best of our knowledge, does not pose any privacy concern and has been used widely by the ML community
freely. The simulations were performed on HPC cluster. All of the models were trained on single 110 (80GB) GPU. We use Adam optimizer with default parameters for minimizing the loss function.
32

---

I.1 Static MNIST input
We performed a random split of 57,000 training samples and 3,000 validation samples and picked the model with the largest validation accuracy for testing. To make direct comparison to SSN [55], we
use the same architecture structure as theirs. First, we train an autoencoder (Table 3 to reduce the dimensionality of MNIST images to 40 by using mean -aaaand loss function. Then, 40-dimensional vector as input to ORGaNICs and train it using the cross-entropy loss function. We additionally make the input gain b dependent on the input x, b f(Wixx), where is sigmoid. A layer of ORGaNICs is given by Eq. 15
this
W xx + (1 Vlal) (W.y)
Ty yy = + bxX O
(134)
- Ta 9 -
W (y 0 al)
The "output" of layer is the steady-state state firing rate of the neuron with the positive receptive field, ie., y = = We parameterize W, to have a maximum singular value and instead of simulating the dynamical system to find the fixed point, use the iterative Algorithm [ with maximum number of steps 10. More details about the parameters are given in Table 4 kaiming
uniform initialization is used from (94). Additional hyperparameters are given in Table 7 We train
ORGaNICs a single-layer setting with the number of y neurons encoding the input, N₁ = 50, 80.
We also train two-layer ORGaNICs (Table with N1 120 and N₂ 60 neurons in each layer. The model is trained using backpropagation and takes approximately 10 min to fully train.
Table 3: Autoencoder architecture
Layer
Shape
Nonlinearity
Input - encoder (layer-1)
784 x 360
ReI.U
encoder (layer-1) >> encoder (layer-2)
360 120
ReL.U
encoder (layer-2)- - embedding
120 x 40
sigmoid
embedding - decoder (layer-l)
40 x 120
ReLU
decoder (layer-1) - decoder (layer-2)
120 x 360
ReLU
decoder (layer-2) 2) -> output
360 x 784
sigmoid
Table 4: ORGaNIC' parametrization for static MNIST classification
Parameter
Shape
Learned
Initialization
W
N x M
yes
kaiming uniform
W br
N X M
yes
kaiming uniform
W
N X N
yes
identity
W
N N
yes
ones
bo
N
yes
random normal
N
no
ones
Table 5: ORGaNIC's architecture for static MNIST classification
Layer
Shape
Nonlincarity
Input + ORGaNICs (layer-1) ORGaNICs (layer-I) + ORGaNICs (layer-2)
40 x N
None
N1 X N2
None
ORGaNIC' (layer-2) + fully-connected
N₂ 10
None
1.2 Permuted and Unpermuted sequential MNIST
We performed a random split of 57,000 training samples and 3,000 validation samples and picked the model with the largest validation accuracy for testing. The unpermuted sequential MNIST task defined as follows: for a given 28 28 image, we flatten it to get one-dimensional, 784 timestep-long input. Then these pixels are presented as an input (X,, one pixel each time-step i) to the Fuler discretized rectified ORGaNICs model (Eq. 87 with rectified input drive, given by the following
33

---

equations,
Af
Yi+1 = Yi +
(-y; + b, [ x + (1 a o [wyil)
At
a;+1 = a; +
(-a, b 2 + W (y, oa.²))
Ta
(135)
1 = b; +
C (-b; f Wi X W Wy a))):
At
bo,441 bo.
(-bo, f(W bn + Who Woyyy + W Whaa:))
The
When we are done presenting the pixels we use the last hidden state, i... Y784, to make the predictions. To make this more challenging also train ORGaNIC' on permuted sMNIST where we first permute the pixels of all the images in some random order and the rest of the task is the same. Instead of parametrizing т. we parametrize At/Ty 0.05 f(py). At/T 0.01 * f(Pa) and 0.1 ( so we can control the dimensionless relative time
At/T = 0.1 * f (Ps) and Át/Tb,, constants. In practice, we find better to make the a neurons sluggish compared to y. This based on the intuition given by the two-dimensional phase portrait for different relative time constants Fig- 3 All the parameters (including W,) are unconstrained for this task with initialization specified stable, we did not need to gradient clipping for training, which is commonly used for LSTMs. Additionally, we train the model using a StepLR learning rate scheduler The model is trained using backpropagation through time (BPTT) and takes approximately 30 hours to fully train.
in Table 6 Since ORGaNICs
with parameters given in Table 7
Table ORGaNICs parametrization for sequential MNIST classification
Parameter
Shape
Learned
Initialization
W
N
yes
kaiming uniform
W ba
N X
yes
kaiming uniform
W by
N
yes
kaiming uniform
W ba
N X
N
yes
kaiming uniform
byu
N X
yes
kaiming uniform
W buy
N X N
yes
kaiming uniform
Wiya bya
X
yes
kaiming uniform
W.
N X N
yes
identity
W
N X
yes
ones
N
ones
Table 7: Hyperparameters
Hyperparameter
Static MNIST
Sequential MNIST
Batch size
256
256
Initial Learning rate
0.001
0.01
Weight decay
10-5
10-3
Step size (StepLR)
None
30 epochs
Gamma (StepLR)
None
0.8
34

---

Supplementary figures
0 < uy S 1
Unstable fixed point
Wr > 1
Stable fixed point
byar 1//,,
byo 1/us,
a
9
0.
Figure 2: Phase portraits for 2D OR GaGII's with negative input drive. We plot the phase portraits of 2D ORGaNICs the vicinity of the stable fixed point for contractive (a, d) and expansive (b, c, e, f) recurrence scalar w,. A stable fixed point always exists, regardless of the parameter values. (a-c), The main model (Eq. [6). (d-f), The rectified model (Eq, [02). Red stars and black circles indicate stable and unstable fixed points, respectively. The parameters for all plots are: b 0.5, Ta 2 ms, Ty ms, w = 1.0, and 2 1.0. For (a) & (d), the parameters are w, = 0.5, by 0.5, o 0.1; for (b) & (e), Wr 2.0, 60 0.5, 0.1; and for (c) & (f), w, 2.0, bo 1.0, 1.0.
a
Ty Ta
b
Ty
Ty
2
0
19
Figure 3: Phase portraits for 2D rectified ORGaNICs for different time constants. Red stars indicate stable fixed points. The parameters for all plots are: w, 1.0, by 0.5, b = 0.5, 0.1. w 1.0, and z = 1.0. For (a), the time constants are Ta ms, ms; for (b), Ta = 10 Ty ms: for (c), 2 ms, Ty 10 ms.
35

---

b
muste
masied
122
Error (c)
+
02
C
Numberr of Itenstons
0.5
M4
1,5
G
sunnaaiion
5.0
iteration
8
0 1-3
y 2011
a (3p
-00
0.25
-00
Ttenution 100
Iteration 100
Neuron iodex
Nouron index
Figure 4: Fast convergence of the iterative algorithm. The results for 20-dimensional ORGaN- ICs y and 10 a neurons) with random parameters and inputs with the additional constraint the maximum singular value of W, equal to 1 and ||||| < 1-sigma S.D.) and maximum errors (e) as a function of number of iterations. is calculated as the norm of the difference between the solution (found by simulation starting with random initial- ization) and the iteration solution. (b), An example of a randomly sampled W,, (c), Steady-state approximation function of iteration number. Different lines represent different neurons. (d), Overlap between the iteration solution (after 15 iterations) and the true solution.
(a), Mean (with error bars representing
b
240
150
Frequency
160
100
50
80
- -0.285
-0270
-0.255
-0.240
0.200
-0.175
-0.150
-0.125
-0.100
max(Re(A))
max(Re(A))
Figure 5: Histogram for the eigenvalue with the largest real part. We train two-layer ORGaNICs (To Ty 2 ms) with static MNIST input where W, is constrained to have maximum singular value 1. We plot the histogram of cigenvalues of the Jacobian matrix with the largest real part, for inputs from the test set. We find that all the eigenvalues of the Jacobian have negative real parts, implying asymptotic stability. (a), histogram for the first layer. (b), histogram for the second layer. Note that since this is implemented in feedforward manner, this is cascading system with feedback, hence we can perform the stability analysis of the two layers independently.
36

---

0.0
first layer
sccond Inyer
-0.1
max(X)
-0.2
01.3
3000
6000
9000
Training steps
Figure 6: Eigenvalue with the largest real part while training on static input (MNIST) classifi- cation task. This plot shows the largest real part of eigenvalues across all test samples as training progresses. The fact that the largest real part consistently remains below zero indicates that the system maintains stability throughout training.
a
b
y
-2
-2
200
400
600
800
200
400
600
800
Time
Time
Figure 7: Trajectories the hidden states (y). This plot shows the dynamics of the hidden state as the input is being presented sequentially. We train ORGaNICs (128 units) RNN (a), unpermuted sequential MNIST and (b), permuted sequential MNIST. The inputs are picked randomly from the test set. The hidden state trajectory remains bounded, indicating stability.
37

---

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: We provide proof, empirical evidence, and citations of any claims made in the abstract and the introduction.
2. Limitations
Question: Does the paper discuss the limitations of the work performed by the authors?
Answer: [Yes]
Justification: We discuss the limitations of the model in the Discussions section.
3. Theory Assumptions and Proofs
Question: For each theoretical result, does the paper provide the full set of assumptions and a complete (and correct) proof?
Answer: [Yes]
Justification: Detailed proofs all the theorems are provided in the Appendix.
4. Experimental Result Reproducibility
Question: Does the paper fully disclose all the information needed to reproduce the main ex- perimental results of the paper to the extent that it affects the main claims and/or conclusions of the paper (regardless whether the code and data are provided or not)?
Answer: [Yes]
Justification: We describe the architecture in detail and provide code to implement it in PyTorch.
5. Open access to data and code
Question: Does the paper provide open access the data and code, with sufficient instruc- tions to faithfully reproduce the main experimental results, as described in supplemental material?
Answer: [Yes]
Justification: We provide URL to the GitHub repository with the code used reproduce all the results the paper.
6. Experimental Setting/eetails
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: This information can be found in the supplementary.
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 in the plot dealing with the random initialization of the model.
8. Experiments Compute Resources
Question: For each experiment, does the paper provide sufficient information on the com- puter resources (type of compute workers, memory, time of exccution) needed to reproduce the experiments?
Answer: [Yes]
38

---

Justification: We provide all of the information in the supplementary.
9. Code Of Ethies
Question: Does the research conducted the paper conform, in every respect, with the NeurlPS Code of Ethics
Answer: [Yes]
Justification: We do not violate any code of ethics listed in the guidelines.
Broader Impacts
Question: Does the paper discuss both potential positive societal impacts and negative societal impacts of the work performed?
Answer: [NA]
Justification:
There is no societal impact of the work performed.
11. Safeguards
Does the paper deseribe safeguards that have been put in place for responsible release of data or models that have a high risk for misuse (c.g., pretrained language models, image generators, scraped datasets)?
Question:
Answer: [NA]
Justification: The model used in the paper have been trained on toy datasets and therefore pose risk.
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: cited creators the publicly available MNIST dataset.
13. New Assets
Question: Are assets introduced in the paper well documented and is the documentation provided alongside the assets?
Answer: [Yes]
Justification: We provided an anonymized URL in the supplementary for a well-documented GitHub repository to reproduce all the results in the manuscript.
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 as details about compensation (if any)?
Answer: [NA]
Justification: The paper does not involve crowdsoureing nor research with human subjects. 15. Institutional Review Board (IRB) Approvals or Equivalent for Research with Human Subjects
Question: Does the paper describe potential risks incurred by study participants, whether such risks were disclosed to the subjects. and whether Institutional Review Board (IRB) approvals (or equivalent approval/review based the requirements of your country or institution) were obtained?
Answer: [NA]
Justification: The paper does not involve crowdsourcing nor research with human subjects.
39


Reviews:
[{"review_id": "LzzVJudoog", "reviewer": "Reviewer_NGE1", "summary": "Recurrent neural networks are widely applied in both solving machine learning tasks and modeling neural recordings.  However, they can be difficult to train due to exploding/vanishing gradients and other sources of instability, often requiring hyperparameter optimization.  In this paper, the authors argue that the ORGaNICs model is capable of solving sequential tasks without the need for hyperprameter optimization due to the conjectured property of ``unconditional stability,'' which guarantees the existence of a stable fixed point in the dynamical system for any model parameters.  The biological plausibility of the ORGaNICs dynamics also suggests that a similar dynamical system could guarantee stability in biological neural circuits for processing sensory information and generating motor sequences.\n\nThe stability property is proven rigorously in the case where the recurrent weights are equal to the identity and in the case where only a single recurrent and inhibitory unit are included in the model.  The authors provide new insights into the behavior of ORGaNICs by providing a Lyapunov function for the dynamics near the fixed point in the case where $\\mathbf{W}_r = \\mathbf{I}$.  An algorithm for quickly finding a stable fixed point of the ORGaNICs dynamical system in general is provided.  This algorithm is used to probe empirically for the fixed points of the ORGaNICs dynamics with general recurrent weights $\\mathbf{W}_r$.\n\nThe authors test their prediction empirically in numerical experiments where they train an ORGaNICs network on two tasks: Static and sequential MNIST classification.  In both cases, they show that ORGaNICs obtains competitive performance without hyperparameter optimization.", "questions": "What hyperparameters are present in the ORGaNICs model as trained on the static and sequential MNIST tasks?  While the time constants are learned, how is the learning rate for naive BPTT selected?  Would this hyperparameter plausibly need to be optimized to achieve performance competitive with a well-tuned LSTM?\n\nIn the verification of the iterative fixed-point finding algorithm, how is the ``true'' fixed point determined?  By running the dynamics from a random starting point to the fixed point, or by running the dynamics from an initial estimate obtained using the iterative fixed point finding algorithm?", "limitations": "Even for static inputs, there are scenarios in which a biological circuit may want to produce a non-static output.  For example, in central pattern generators of motor sequences. If indeed ORGaNICs always converges to a stable fixed point, this may represent too constrained of a circuit to explain many interesting neural phenomena.", "rating": 6, "confidence": 2, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "This paper includes many interesting insights and clever mathematical arguments.  The proof of unconditional stability through the indirect method of Lyapunov and by relating the problem to a second-order dynamical system is mathematically sound and provides interesting insights through the analogy to a damped harmonic oscillator.  The finding that ORGaNICs achieves competitive performance without the need for hyperparameter optimization  on the static and sequential MNIST tasks is also very promising.", "weaknesses": "The main theorem makes the very restrictive assumption that the recurrent weights $\\mathbf{W}_r = \\mathbf{I}$, such that they exactly cancel the intrinsic decay rate of the recurrent neurons.  While there are neural circuits where all recurrent interactions take place through inhibitory interneurons (e.g. cerebellar granule cells interacting through inhibitory Golgi cells), the intrinsic decay must also be taken into account.  If the assumption can be relaxed to say that $\\mathbf{W}_r = \\alpha \\mathbf{I}$ for some $\\alpha \\in [0, 1]$ then I would be inclined to raise my score.\n\nAlso, the argument for unconditional stability of the ORGaNICs model is not well-supported in the paper's current form.  Beyond the case $\\mathbf{W}_r = \\mathbf{I}$ and the two-dimensional case, there this is tested only for the MNIST task with constant inputs.  In these tests, it also appears that the fixed-point-finding algorithm is used to quickly discover the stable fixed point.  However, it would be stronger to show that the dynamics naturally converge to the fixed point in a reasonable amount of time when started from randomly generated initial configuration.  In fact, it is possible that the dynamics naturally lead the system to a stable limit cycle and never converge to the stable fixed point, even if it is guaranteed to exist.  We ask the authors to clarify the methods used to locate the fixed point of the ORGaNICs network in their experiments and show that the system's dynamics tend toward this fixed point and not another limit cycle in a reasonable amount of time.  The numerical training also enforces a maximum singular value of $\\mathbf{W}_r$, which is a strong condition required for stability.  ``Towards Unconditional Stability...'' would be a more appropriate title given the strong constraints required for stability guarantees in both the analytical proofs and numerical tests.\n\nFor the sequential MNIST task, the authors claim that they are able to achieve LSTM-like performance without hyperparameter optimization or gradient clipping/scaling.  During these experiments, the neural time constants are treated as parameters that can be learned.  This result is different than showing that for any fixed values of the time constants, ORGaNICs can be trained without issue.  The authors should verify that fixing the time constants and then training also leads stable training and no need for gradient clipping.  This is important to maintain a claim of biological plausibility where time constants are intrinsic properties of the neurons and not trained.\n\nThe static and sequential MNIST tasks also represent a rather limited scope in which to test for unconditional stability in practice.  It's possible (though I agree, unlikely) that for other more unwieldy tasks unconditional stability will not hold.  Either a proof of unconditional stability for general recurrent weight matrices or a wider range of numerical tests are necessary to make the case for unconditional stability of the ORGaNICs dynamics.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "**Weaknesses:**\n- “The main … raise my score”. **Answer:** Our proof for the stability of the high-dimensional model (using the indirect method of Lyapunov) relies on the existence of analytical expressions for the fixed point. Such an expression exists when $\\mathbf{W}_r=\\mathbf{I}$, but does not exist when $\\mathbf{W}_r \\neq \\mathbf{I}$. However, we were able to prove unconditional stability for the 2-dimensional model (see Fig.1, Theorem.5.1 and Theorem.5.2), for any positive value of $\\alpha (w_r) \\in (0,\\infty)$ and any combination of choices for other parameters. Building on these two limiting cases, for which a rigorous proof of stability is tractable (i.e., for any $w_r$ in mean field, when d=2, and for $\\mathbf{W}_r=\\mathbf{I}$ when d>2, with no restrictions on other parameters), we conjecture at the end of Section 5 that the stability property holds in higher dimensions for systems where the largest singular value of $\\mathbf{W}_r$​ is less than 1 (therefore covering the case when $\\alpha \\in [0, 1]$). We provide empirical evidence in support of this conjecture in Fig. 4 and Fig. 5. To satisfy the reviewer’s question directly, we simulated 10,000 random networks with  $\\mathbf{W}_r=\\alpha \\mathbf{I}$ with $\\alpha \\in [0,1]$ and plot the distribution of the maximum real part of eigenvalues, all of which are found to be less than 0 (Fig.1 of pdf), indicating stability. We do not exclude that a statistical approach could be used to prove the stability of a general network “in expectation”, but 1) it would fall short of a rigorous proof of stability like ours; 2) it is sufficiently ambitious that it should be the subject of a subsequent paper.\n- “Also, the  … amount of time”. **Answer:** We would like to point out to the reviewer an additional experiment at the end of Section 5 where we mention empirical evidence of stability for a more general case of the recurrent matrix. where the maximum singular value of $\\mathbf{W}_r$ is constrained to 1. This conjecture is supported by empirical evidence showing consistent stability, as ORGaNICs initialized with random parameters and inputs under these constraints have exhibited stability in 100\\% of trials Fig.4.” In Fig. 4 the “true”  fixed point was found by simulating the network from a zero initial condition. In all of the trials where the conditions in the conjecture were met, the dynamics always converged to the stable fixed point found by the iterative scheme. To show that in all simulations there is no limit cycle and that we can start from a random initial condition, we re-performed this analysis (10,000 networks), starting from random initial conditions \\in [0, 1] and fixed $\\mathbf{W}_r=\\alpha \\mathbf{I}$ for some $\\alpha \\in [0,1]$. We find that the system is stable (Fig.1 of pdf) and the simulations always converge to the same fixed point as found by the iterative scheme. So, while we cannot prove the conjecture by the indirect method of Lyapunov, numerical evidence suggests that the conjecture is true almost surely.\n- “For the … and not trained”. **Answer:** This is an excellent point. We chose the time constants to be learnable to demonstrate that ORGaNICs is competitive with LSTMs on sequential MNIST tasks without hyperparameter optimization and without using ad hoc techniques for training, such as gradient clipping/scaling. As per the referee's suggestion, we have now conducted additional experiments where ORGaNICs is trained with fixed time constants (Table 1 of pdf). Our findings indicate that even with fixed time constants, ORGaNICs maintain stable training without the need for gradient clipping/scaling.\n\n- “The static and … ORGaNICs dynamics”. **Answer:** We have provided empirical evidence for unconditional stability for a general recurrent weight matrix, particularly when the maximum singular value of $\\mathbf{W}_r$ ​ is less than 1 (Section 5, Fig. 4, and Fig. 5). We prove unconditional stability, using a non-trivial approach, in two specific cases of $\\mathbf{W}_r$ ​: first, when $\\mathbf{W}_r = \\mathbf{I}$, and second, for the 2D model. While a proof for unconditional stability for a general $\\mathbf{W}_r$​ would be most desirable, it is currently beyond reach as detailed in our response above. Therefore, we must rely on empirical evidence for the stability of a general $\\mathbf{W}_r$. In a subsequent version of the paper, we will provide results for the formal language parity task to further corroborate our point.\n\n**Questions:**\n\n- “What hyperparameters … tasks?” **Answer:** Learning rate = 0.001 (static MNIST), 0.01 (sMNIST); Batch size = 256, Weight decay = $10^{-5}$ for both. This information has been added in the Appendix.\n\n- “While the … LSTM?” **Answer:** We do not have a specific scheme for selecting an optimal learning rate. Since the point of this paper is to showcase greater stability, trainability, and interpretability compared to other RNNs, and not to achieve SOTA on the benchmarks, we did not tune the hyperparameters. This will be explored in future work and should lead to better performance on ML benchmarks.\n\n- “In the verification … algorithm?” **Answer:** The “true” fixed point is found by simulating the network from a zero initial condition. However, in response to the reviewer, we have verified that for a random initialization as well, the simulation converges to the same fixed point. \n\n**Limitations:**\n\n- “Even for … phenomena”. **Answer:** This class of models has already been shown to produce non-static oscillatory activity as well as to reproduce a wide range of neural phenomena, including sustained activity, sequential activity, motor preparation and motor control (Heeger, D. J., & Mackey, W. E. (2019)); simulate the dynamics of V1 activity (Heeger, D. J., & Zemlianova, K. O. (2020), S. Rawat, D.J. Heeger, and S. Martiniani. Cosyne Abstracts 2024); predict the emergence of communication subspaces in interareal communication (S. Rawat, D.J. Heeger, and S. Martiniani. Cosyne Abstracts 2023).", "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": "0frc7g6k9X", "reviewer": "Reviewer_2CUH", "summary": "The paper discusses the development and analysis of \"Oscillatory Recurrent Gated Neural Integrator Circuits\" (ORGaNICs), a biologically plausible model of recurrent neural networks that implements divisive normalization. The authors prove the unconditional local stability of ORGaNICs with an identity recurrent weight matrix using the indirect method of Lyapunov function. They also demonstrate empirical stability for higher-dimensional circuits. The model's performance is evaluated on static and dynamic classification tasks, showing comparable results to LSTMs without specialized training strategies, thanks to its intrinsic stability properties.", "questions": "1. Can this model provide more interpretability in relevant machine learning tasks? \n2. Besides normalization and stability, can this E-I model offer more insights into neuroscience compared to machine learning models?\n3. Many existing machine learning models have already implemented various forms of divisive like normalization such as layer normalization. What are advantages or innovations of the divisive normalization implemented in this model.\n4. In comparison to standard machine learning models, biophysically plausible models typically offer superior interpretability. However, the paper offers limited discussion on this topic.", "limitations": "1. In the sequence modeling benchmarks, only 2 tasks are evaluated. It is recommended to include a broader range of sequence modeling tasks to allow for more comprehensive comparisons, such as those where SSM excels in long sequence modeling or where LSTM performs well in formal language tasks like the parity task.", "rating": 6, "confidence": 4, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "1. Biological Plausibility: ORGaNICs are designed with a structure and dynamics that are more aligned with biological neural circuits than traditional artificial neural networks.\n2. Biological divisive normalization: Implemented divisive normalization using a more bioplausible network.\n3. Trainability: ORGaNICs are a biophysically plausible extensions of LSTM and can be directly trained by BPTT.", "weaknesses": "1. Generalization: The paper does not extensively discuss how well the results might generalize to more complex or different types of tasks beyond the tested benchmarks.\n2. Performance on Benchmarks: The proposed model performance still has a gap compared to machine learning models.\n3. Scalability: parameter size scales too fast with model size and there is no obvious way for impovement.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "**Weaknesses:**\n\n- “Generalization: The paper ... the tested benchmarks”. **Answer:** We expect ORGaNICs to perform well on other ML benchmarks, especially those concerning sequential data. This will be done in a future study, therefore we refrain from making any specific claims about performance on other benchmarks. In this paper, we sought to demonstrate the stability, trainability, and interpretability of ORGaNICs compared to other RNNs, and present proofs and evidence supporting these claims. Specifically, we were able to rigorously prove for the two limiting cases that ORGaNICs is absolutely stable and demonstrated empirically that this stability holds broadly. We achieved good performance on ML tasks using naive BPTT without resorting to gradient clipping/scaling or other techniques required for alternative neurobiological models like SSN, and without hyperparameter tuning. We derived a Lyapunov function for our model offering an interpretable optimization principle explaining the dynamics of ORGaNICs.\n\n- “Performance on Benchmarks: … machine learning models. **Answer:** There was no hyperparameter tuning nor an attempt to get SOTA performance in this paper. The top-performing ML models are designed with properties of low-dissipation and stability to solve the problem of vanishing and exploding gradients. In contrast, ORGaNICs is a biologically plausible circuit designed to implement divisive normalization (DN), a computational principle experimentally found in a wide range of cortical areas (i.e., different brain regions) across different species, to simulate cortical activity. Unlike ML models, stability in ORGaNICs is not engineered but it emerges from its neurobiological plausible design. Despite not being optimized for the tasks considered, ORGaNICs performs competitively with other RNN models. To our knowledge, no other neurodynamical biologically plausible model achieves this.\n\n- “Scalability: parameter size scales … way for impovement”. **Answer:** In future work, we will design the gains $\\mathbf{b}$ and $\\mathbf{b_0}$ to be dependent locally on $\\mathbf{y}$ and $\\mathbf{a}$, This will significantly reduce the number of parameters compared to the current model where the matrices $\\mathbf{W_{by}}$, $ \\mathbf{W_{ba}}$, $ \\mathbf{W_{b_0 y}}$ and $ \\mathbf{W_{b_0 a}}$ are $n \\times n$.\n\n**Questions:**\n\n- “Can this model provide … machine learning models?” **Answer:** We draw a direct connection between ORGaNICs and systems of coupled damped harmonic oscillators, which have been studied in mechanics and control theory for decades. This connection allows us to derive an interpretable energy function for a high-dimensional ORGaNICs circuit (Eq. 128), providing a normative principle of what the circuit aims to accomplish (see Eq. 13 and the subsequent paragraph). For a relevant ML task, having analytical expressions for the energy function (which is minimized by the dynamics of ORGaNICs) allows us to quantify the relative contributions of the individual neurons in the trained model, offering more interpretability than other RNN architectures. For instance, Eq. 128 reveals that the ratio $\\tau_y / \\tau_a$ of a neuron pair ('y' and its corresponding 'a') indicates the \"weight\" a neuron assigns to normalization relative to aligning its responses to the input. This provides a clear functional role for neurons in the trained model. Furthermore, since ORGaNICs is biologically plausible, we understand how the different terms of the dynamical system may be computed biologically within a neural circuit (Heeger & Mackey, 2019). This bridges the gap between theoretical models and biological implementation, providing a framework to test hypotheses about neural computation in real biological systems.\n\n\n- “Many existing … implemented in this model”. **Answer:** Divisive normalization (DN) was introduced as a model of the steady-state response of neurons, functioning as a static nonlinearity similar to batch and layer normalization. The DN equation has also been shown to generalize batch and layer normalization, as detailed in the work by Ren, M., Liao, R., Urtasun, R., Sinz, F. H., & Zemel, R. S. (2016). However, in the brain, there are no static nonlinearities; it has thus been proposed that DN is achieved via a recurrent circuit. ORGaNICs is such a recurrent circuit designed so that the responses of neurons at steady state follow the DN equation. While batch and layer normalization are ad hoc implementations that do not affect the dynamics (they are applied to the output layer), ORGaNICs implements DN dynamically. Additionally, whereas batch and layer normalization do not inherently affect the stability of the model (because they do not influence the dynamics), our paper demonstrates that a model implementing DN naturally exhibits stability, which is greatly advantageous for trainability. This stability, derived from the dynamic implementation of DN, sets ORGaNICs apart by providing both output normalization and model robustness.\n\n- “In comparison… limited discussion on this topic”. **Answer:** We thank the reviewer for this suggestion. We will add a paragraph based on our answer to the first question. \n\n**Limitations:**\n\n- “In the sequence modeling …  the parity task”. **Answer:**  We thank the reviewer for the suggestions regarding additional tasks. To satisfy the reviewer, in a subsequent version of the paper, we will include ORGaNICs’ performance on the parity task. However, we stress that the primary focus of this paper is to demonstrate the stability, trainability, and interpretability of ORGaNICs compared to other RNNs, which we have achieved through a thorough analytical and numerical analysis, and not to provide an extensive benchmark of the model on ML tasks. A more comprehensive test of ORGaNICs on additional long sequence modeling tasks, along with hyperparameter optimization, will be pursued in future work.", "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": "jc0nESpYlH", "reviewer": "Reviewer_FwjG", "summary": "This paper studies the stability properties of a model of cortical circuits which was introduced in 2019 (and I wasn't yet aware of): the ORGaNIC model by Heeger & Mackey. This LSTM-like model uses a simple set of differential equations that unify several phenomena observed in cortex, including normalization (Carandini & Heeger, 2012), and in some sense map LSTMs onto plausible cortical circuitry. This paper presents two main theoretical results: (i) that the model is locally stable around its unique fixed point when the recurrent matrix is the identity, and that (ii) existence and uniqueness of a stable fixed point can be shown for a general 2D model (really just one \"principal neuron\", and a gate variable) under certain conditions on the parameters. The paper wraps up with a few training experiments on (regular/sequential/permuted) MNIST classification, showing that ORGaNICs perform well.", "questions": "- The abstract (and main text) says you trained your models using \"simple backpropagation through time\" but appendix I says you trained all models using Adam -- in principle there is no contradiction between these two claims, as backprop is also involved in Adam, but I wonder why you wrote \"simple backpropagation\" (which could easily be misunderstood as \"simple gradient descent\", i.e. not Adam). What's non-simple backpropagation?\n\n- The formatting of references is annoying :D please use parentheses to clearly separate references from the main text.\n\n- typo on l104: \"dynamics-neural growth\" → \"dynamics-neutral growth\"", "limitations": "Some technical limitations are stated in the Discussion section.\n\n*EDIT*: following the rebuttal, I am raising my score to a 5.", "rating": 5, "confidence": 2, "soundness": 3, "presentation": 2, "contribution": 2, "strengths": "This paper is very strong on a technical level. The maths are sound as far as I could tell from going through the paper carefully (mostly main text; took only a cursory look at the appendices) but not re-deriving the equations myself. The proof of local stability in ORGaNICs when $W_r = I$ is elegant and actually rather sophisticated; perhaps some of the most useful collaterals of this proof are (i) the explicit derivation of the eigenvalues of the Jacobian when additionally all normalization weights are equal to the same (positive) value, and (ii) the energy function of Eq 13 that provides further insight into the dynamics of ORGaNICs. Both of these add to the interpretability of this model class. The iterative fixed-point-finding algorithm (Algorithm 1) seems novel and useful, too.", "weaknesses": "Overall this is a highly esoteric paper for which I have difficulty assessing potential impact (hence my high uncertainty rating). The theoretical results are impressively detailed but seem fairly limited in scope (how are the W_r = I and 2D cases relevant for either ML or neuroscience?). I am also not super convinced by the utility of the broader conjecture on stability at the end of section 5. The experiments are fairly limited in scope and breadth, too (e.g. no mention of hyperparameter tuning for the SSN comparison).\n\nThe abstract claims that it is “thanks to its intrinsic stability property” that ORGaNICs perform comparably to LSTMs; while this seems like a sensible hypothesis I don't think the paper really shows that, yet I believe the paper would be much stronger if that was shown to be true. In this respect, I wonder in what way the comparison to SSNs (e.g. Soo & Lengyel) helps make this case; is it because SSNs also perform a form of normalization at steady-state, much like ORGaNICs, but have otherwise no stability guarantees? Is it even true that SSN training is brittle / often fails because the network loses stability, thus giving rise to a stronger-than-usual tradeoff between training stability and learning rate magnitude? How could the author rule out that there might be another fundamental difference between these 2 models that has little to do with stability and yet underlies the performance difference? After all, no formal stability guarantees exist for standard LSTMs as far as I am aware, and that doesn't prevent them from training very well on most tasks (was gradient clipping even necessary for the LSTMs in those particular experiments?). In summary, I suppose the paper lacks a couple of convincing ablations to make the point that stability guarantees increase training performance.\n\nOn clarity of exposition: this is a pretty hard paper to follow, it's very dense and doesn't really offer the hierarchy of exposition that a reader would want to see -- the theoretical results are presented in a very flat way, with a permanent back and forth between details and bigger picture that I found hard to navigate.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "**Weaknesses:**\n\n- “Overall this is a … uncertainty rating)”. **Answer:** \nWe argue that the potential impact is high. LSTMs and GRUs have had a huge impact on ML/AI applications even though they are not always stable and hence require ad hoc techniques for training (e.g., gradient clipping/scaling). Furthermore, although normalization has been shown to improve the training and performance of CNNs, there has been no principled way of adding normalization to RNNs (putting a normalization layer between RNN layers is not the same as normalizing the activations within each RNN layer). ORGaNICs overcomes these limitations, with comparable capabilities to LSTMs and GRUs, with built-in recurrent normalization (motivated by neurobiology) that, as we show, is sufficient to guarantee dynamical stability. Moreover, we note also that normalization is an essential architectural motif in modern ML, and that divisive normalization has been shown to generalize batch and layer normalization, as shown in the work by Ren, M., et. al., (2016). As such, the ability to integrate normalization in the dynamics of an RNN, and to analytically demonstrate the impact of normalization through the derivation of an interpretable normative principle (Lyapunov function) is far from esoteric and very much needed in the theory of deep learning. Finally, our work establishes an important precedent for how the incorporation of neurobiological principles can drive advances in ML.\n\n- “The theoretical … end of section 5”. **Answer:** Considering the cases of $\\mathbf{W}_r = \\mathbf{I}$ and the two-dimensional model, we explore two different limits of arbitrary $\\mathbf{W}_r$. The first limit involves relaxing constraints on all parameters in a high-dimensional model while assuming $\\mathbf{W}_r = \\mathbf{I}$. The second limit involves relaxing constraints on all parameters, including $w_r$, for a 2D model. These two limits provide a foundation for understanding the empirical results regarding stability in models with arbitrary recurrence, which is intractable with the approach used in the paper. There is also a rich history of studying neural mean field theories (Wilson, H. R., & Cowan, J. D. (1972); Kraynyukova, N., & Tchumatchenko, T. (2018)) which yield two-dimensional E-I circuits used for modeling the cortex. The conjecture is crucially important because we want to learn an arbitrary recurrent weight matrix during ML and neuroscience tasks. We have empirical evidence in support of this conjecture (see Fig.4,5 and Fig.1,2 of the attached pdf). We also use this conjecture in the paper itself (for the static input task). \n\n- “The experiments … SSN comparison)”. **Answer:** We thank the reviewer for the suggestions regarding additional tasks. To satisfy the reviewer, in a subsequent version of the paper, we will include ORGaNICs’ performance on the parity task for comparison. However, we stress that the primary focus of this paper is to demonstrate the stability, trainability, and interpretability of ORGaNICs compared to other RNNs. There was no hyperparameter tuning for any of the tasks – ORGaNICs simply outperforms SSNs on the first attempt at the task we considered. \n\n- “The abstract claims… to be true”. **Answer:** By that statement, we mean to convey that ORGaNICs mitigates the problem of exploding gradients: through its inherent stability; and vanishing gradients: by processing information across various timescales while maintaining stability, resulting in a blend of lossy and non-lossy neurons. This leads to a trainable RNN circuit, with a performance comparable to LSTMs without the need for specialized techniques like gradient clipping/scaling. We have revised this sentence for better clarity.\n\n- “In this respect, ...  training performance”. **Answer:** Like SSN, ORGaNICs is a neurodynamical model originally designed to explain cortical activity. The comparison with SSN is simply meant to demonstrate that just like SSN, ORGaNICs can be successfully trained on a static MNIST task. So we are not trying to make the point that ORGaNICs performs better due to its stability guarantees, but rather that ORGaNICs can be trained by naive BPTT and perform comparably to SSN trained by the non-trivial method of dynamics-neutral growth (DNG). The result of the comparison is that ORGaNICs performs slightly better (despite no hyperparameter tuning) while being much easier to train than SSN. If we were to speculate on the performance gap between ORGaNICs and SSN, our best guess would be that DNG constrains the range of possible models (viz. parameters) thus reducing the expressivity of SSN. It is beyond the scope of this paper to identify why SSN trained by DNG does not do as well, but what is certain is that it is much harder to train SSN than ORGaNICs. The main takeaway is that ORGaNICs can be trained naively because it is stable, and it performs as well or better than alternative neurodynamical models (e.g., SSN) trained by sophisticated techniques.\nRegarding the LSTMs, they were trained using gradient clipping in Arjovsky, M., Shah, A., & Bengio, Y. (2016).\n\n- “On clarity of exposition: … navigate”. **Answer:** We thank the reviewer for the suggestion. We will revise the introduction of the paper to provide a clearer roadmap of the main results of the paper and how they are organized.\n\n**Questions:**\n\n- “The abstract … backpropagation?” **Answer:**   We mean to convey the fact that we train ORGaNICs by naive backpropagation (not simple gradient descent), meaning that 1) we do not use gradient clipping/scaling that is typically adopted during BPTT when training RNNs such as LSTMs; 2) we do not need to resort to specialized techniques such as the sophisticated DNG method used to train SNN, an alternative neurodynamical model. We have resolved this ambiguity by revising the statement.\n\n- “The formatting … text”. **Answer:** We thank the reviewer for the suggestion, we have updated the references using the suggested format.", "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": "tzv9fE202F", "reviewer": "Reviewer_8Tz2", "summary": "This work is based on ORGaNICs, a particular type of RNN architecture that implements divisive normalization (from neuroscience). The authors explore whether ORGaNICs are stable enough to be meaningfully (and stably) trained by gradient descent. Due to their stability, which the authors have proven theoretically for specific cases and empirically in general, the authors claim that ORGaNICs can learn long-term dependencies and solve vanishing and exploding gradients found in other neurodynamical models.", "questions": "ORGaNIC begins with the word \"oscillatory\". Can the authors comment on whether the network must always be oscillatory? I see that the analysis looks into stable and unstable fixed points, which means that the answer is no, which makes the name confusing to me.\n\nLikewise, can the authors elaborate on the \"gate\" in ORGaNIC and how it is the \"biophysically plausible extension\" (line 86) of gating in machine learning? Specifically, why it is biophysically plausible, and why it is an extension.", "limitations": "In line with my issue of overselling this work, there really is not any honest discussion of genuine limitations in the discussion section, which the authors claim they have done so in the checklist. Discretization problems (lines 350-354) and using single-layers (line 355) are not specific to this work at all. I would encourage the authors to talk about genuine limitations. For example, I suspect but cannot prove, that there would be some increased computational complexity compared to vanilla RNNs, since having multiple differential equation for a single model means that more things need to be done in a single time step.\n\nThis work, focusing on a biologically plausible model, shows the model performing machine learning experiments (in contrast, in Heeger and Mackey [2019] the model is performing neuroscience experiments). \n\nAlso, while this model implements divisive normalization, as I mentioned above there are many biological phenomena that the brain exhibits. It would be unreasonable to ask the authors to state what ORGaNICs cannot do, but at least I would like to ask the authors if ORGaNICs can do everything that SSNs (which they cited) can do? \n\nThese are the limitations that I believe readers would like to see.", "rating": 5, "confidence": 4, "soundness": 3, "presentation": 4, "contribution": 1, "strengths": "Trained RNNs in neuroscience have always had a standard vanilla architecture or have minimal modifications. This work introduces the possibility of training a more complex model motivated by biology while having no stability issues.", "weaknesses": "Subjectively, I believe that the scope of this work is too limited for the audience at NeurIPS, but I will not consider this in my final score. \n\nThe general theme of my objective issues with this work is the slight overselling of this work:\n\n1. Fundamentally, ORGaNIC is an architecture. Song et al. [2016] implemented Dale's law by manipulating the recurrent weight matrix, which is a **method** applicable to any architecture. Soo et al. [2024] trained vanilla RNNs on long-term dependencies using skip connections through time, which once again is a **method** applicable to any architecture. So the claim that these \"**models** do not accurately reflect cortical dynamics\" in lines 105-106 do not make sense to me. To illustrate my point, can the authors comment on whether there are any issues if Soo et al. [2024] is applied to ORGaNICs (beyond the fact that it is not needed)? \n\n2. This work seems to be strongly driving the point that divisive normalization is a key criteria for biological plausibility. There are a variety of biological properties that the brain exhibits. For example, there are 22 different properties of V1 identified in Marques et al. [2020], of which specific versions of DN make up some of them. It is not convincing to select a model based on one particular effect found in V1. \n\nMarques et al. [2020]: https://www.biorxiv.org/content/10.1101/2021.03.01.433495v2 \n\n3. The authors claim that ORGaNICs is stable and therefore solves exploding and vanishing gradients, which does not make sense. The way the authors used the word \"stable\" here is in the language of dynamical systems. \"Stability\" in the context of exploding and vanishing gradients has a different meaning. The fact that the model can learn sequential MNIST (gradient stability) does not mean that it is due to its (dynamical) stability. A completely (dynamically) stable architecture can still have vanishing gradients. Lines 115-116 suggest to me that the authors are somewhat aware of this but not completely.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "**Weaknesses:**\n\n\n- “Fundamentally, ORGaNIC … one particular effect found in V1”.  **Answer:** DN is observed in numerous cortical areas beyond V1 and across different species (Carandini & Heeger, 2012). It explains a wide range of experimental phenomena in various neural systems and cognitive processes (citations in lines 49-53). DN can also be shown to be the appropriate nonlinearity for removing dependencies in sensory inputs, resulting in a set of neural responses that are statistically independent (Simoncelli and Olshausen, 2001). Consequently, DN is not just an effect observed in V1, but a computational principle (Carandini & Heeger, 2012) that can be derived from a statistical optimization criterion. It is thus unsurprising that DN has been shown to dramatically improve the performance of ML models in several applications (e.g., Balleé, Laparra, Simoncelli, 2015; Balleé, Laparra, Simoncelli, 2017). Moreover, batch and layer normalization (of which DN is a generalization, see Ren, M., Liao, R., Urtasun, R., Sinz, F. H., & Zemel, R. S. (2016)) are critically important architectural motifs in modern deep learning architectures. There is thus overwhelming evidence that DN is an essential feature of both natural and artificial neural systems. In this paper we demonstrate that DN can be integrated into the dynamics of an RNN, leading to a provably stable circuit for which an interpretable normative principle (Lyapunov function) can be derived. This work does not just propose an architecture, but provides 1) unprecedented insight into the impact of normalization on the dynamics of RNNs, through an interpretable normative principle; 2) establishes an important precedent for how the incorporation of neurobiological principles can drive advances in ML.\n\n\n- “The authors claim … aware of this but not completely”. **Answer:** In RNN architectures derived from discretized ODEs, such as ORGaNICs, the stability of the dynamical system is intrinsically linked to the vanishing and exploding gradients (VG and EG) problem. Mathematical analysis and details can be found in, Haber and Ruthotto (2017); Chang, Chen, Haber, and Chi (2019); and Erichson et al. (2020). ORGaNICs effectively addresses both problems: EG: Mitigated through the architecture's inherent stability. VG: Addressed by processing information across various timescales (and via the long effective time constant when the membrane time constants are fixed), resulting in a blend of lossy and non-lossy neurons. The effectiveness of ORGaNICs in tackling VG is evidenced by its competitive performance against architectures specifically designed to address this issue, such as LSTMs.\n\n\n**Questions:**\n\n\n- “ORGaNIC begins … name confusing to me”. **Answer:** ORGaNICs is oscillatory in the sense that it can be mapped to a damped harmonic “oscillator” as we show in the paper. This means that for the right parameters and inputs, stochastically driven ORGaNICs exhibit a peak in the power spectrum, which is consistent with the  LFP (Local Field Potential) oscillatory activity observed in neural recordings, even though there are no sustained oscillations.\n\n- “Likewise, can the authors elaborate … it is an extension”. **Answer:** Gating in ORgaNICs is performed by modulating the input gain ($\\mathbf{b}$) and the recurrent gain ($\\mathbf{1} - \\mathbf{a}^+$). Possible mechanisms of how the responses of these neurons may be computed by loops through higher visual cortical areas and/or thalamocortical loops are discussed in the Discussion/Mechanisms of Heeger, D. J., & Zemlianova, K. O. (2020) and its SI Appendix. The gate in ORGaNICs is an extension of the gates in LSTMs/GRUs because the recurrent gain/gate in ORGaNICs (unlike LSTMs/GRUs) is a particular nonlinear function of the output responses/activation, designed to achieve normalization.\n\n\n\n**Limitations:**\n\n\n- “In line with my issue … in a single time step”. **Answer:** The computational cost of ORGaNICs is similar to LSTMs since like LSTMs we have three extra sets of variables ($\\mathbf{a}, \\mathbf{b}, \\mathbf{b}_0$) which are not directly used for prediction. Therefore, ORGaNICs is more computationally expensive than vanilla RNNs. We thank the reviewer for pointing this fact out. We have included this fact in the discussion. However, we do not consider this a limitation, as the increased computational cost results in improved stability, trainability, and interpretability. \n\n- “This work, … neuroscience experiments)”. **Answer:** Since training biologically plausible neurodynamical models using backpropagation is a challenging task (Soo, W., & Lengyel, M. (2022)), we demonstrate the trainability of ORGaNICs by naive backpropagation and backpropagation through time, without gradient clipping/scaling, on well studied ML tasks. Training ORGaNICs on neuroscience experiments will be done in future work.\n\n\n- “Also, while this model implements … SSNs (which they cited) can do?” **Answer:** ORGaNICs has been shown to model a wide range of neural phenomena, including sustained activity, sequential activity, motor preparation, and motor control (Heeger, D. J., & Mackey, W. E. (2019)), as well as simulate the dynamics of V1 activity (Heeger, D. J., & Zemlianova, K. O. (2020); S. Rawat, D.J. Heeger, and S. Martiniani. Cosyne Abstracts 2024) and the emergence of communication subspaces (S. Rawat, D.J. Heeger, and S. Martiniani. Cosyne Abstracts 2023). To address the reviewer’s question, ORGaNICs can do everything that an LSTM can do (while being stable, easily trainable, and interpretable). While a direct comparison of ORGaNICs and SSN across many tasks has not been done, on the one task we consider (static MNIST, following Soo, W., & Lengyel, M. (2022)) ORGaNICs outperforms SSN on the first try (i.e., without hyperparameter optimization). Consistent with our claims, ORGaNICs also performs comparably to LSTMs on sequential tasks such as sequential  (and permuted) MNIST.", "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": "XxxRW8fIaf", "reviewer": "Reviewer_ma5B", "summary": "Stability is a critical notion in the understanding of dynamical systems as well as for learning them. This paper studies the stability of a biologically plausible recurrent cortical model, ORGaNICs, that implements divisive normalization. More precisely:\n- it demonstrates the local stability of two specific subclasses of the model.\n- it shows that the model performs competitively on standard machine learning tasks.", "questions": "- The model considered has the same number of E and I neurons, which is, to the best of my knowledge, far from being the case in the brain. Is the system still stable when the number of I neurons is much smaller than the number of E neurons?", "limitations": "Limitations are properly addressed.", "rating": 7, "confidence": 3, "soundness": 4, "presentation": 3, "contribution": 3, "strengths": "Understanding how properties of biological neural networks lead to better learning properties is an important area of research and this paper nicely contributes to it by highlighting the potential importance of divisive normalization in stabilizing neural dynamics and facilitating training. The paper is overall well-written and nicely connected to both theoretical and empirical neuroscience results. The theoretical stability analysis is non-trivial and provides important insights into the model. The empirical results are solid.", "weaknesses": "The main weakness of the paper is that it might not be fully suited to the NeurIPS community in its current form:\n- the introduction of the dynamics lacks some intuition on what the dynamics achieves. Equation 2 partially does so, but is likely not sufficient to fully convey the intuition. If possible, connecting to existing RNN architectures would help. Introducing the Lypanuov energy earlier could also help.\n- many notions used in Section 4 might deserve a better introduction as NeurIPS is not a dynamical systems conference. This holds for \"Lyapunov diagonally stable\", \"Z-matrix\", \"M-matrix\" and, to some extent, to \"stability\" and \"indirect method of Lyapunov\". To make their results more accessible, the authors may want to provide intuitive definitions of these terms and why they are important.\n\nAdditionally, the connections made to the machine learning literature are sometimes imprecise:\n- in L43, it is written that stability helps mitigate the vanishing gradient problem, which is inaccurate (e.g. a system with very fast time constants is stable but suffers from vanishing gradients). On top of that, stable recurrent networks having slow time constants can still suffer from some form of exploding gradients [Zucchet and Orvieto, 2024](https://arxiv.org/abs/2405.21064).\n- in L61, the authors write that divisive normalization generalizes batch and layer normalization. While there are definitely some links, this is not true: divisive normalization affects the recurrent dynamics whereas both batch and layer normalizations are usually applied after the recurrent layer. In particular, those normalization schemes do not affect the stability of the system.\n- in L349, divisive normalization is directly compared to the form of attention used in Transformers. Transformers normalize over the time axis, which is highly implausible and different from the type of attention mentioned here. This point deserves to be made more precise.\n\nFinally, there are some additional links to existing ML literature that would nicely fit within the paper:\n- the kind of neural network studied in Section 6.1 is known as [Deep equilibrium models, DEQ](https://arxiv.org/abs/1909.01377). One of the main difficulties that come with training such models is to [keep the dynamics stable](https://arxiv.org/abs/2106.14342) throughout learning. The studied architecture may be less prone to these behaviors: adding the DEQ baseline to the experiments and monitoring the stability of the dynamics during training would be interesting additions to Section 6.1.\n- in Section 6.2, the baseline architectures are rather \"old\". The state-of-the-art networks on tasks requiring modeling long-range dependencies are state-space models (e.g. [S4](https://arxiv.org/abs/2111.00396) or [LRU](https://arxiv.org/abs/2303.06349)). Additionally, sequential models are known to be extremely sensitive to initializations on tasks such as sMNIST, see [Amos et al. 2023](https://arxiv.org/abs/2310.02980). It is therefore important to describe the initialization scheme of the recurrence in the main text, as it is likely a main driver of the performance. \n\nI hope these remarks can help the authors to improve their paper.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "**Weaknesses**:\n\n\n- “the introduction of the dynamics … why they are important”. **Answer:** We thank the reviewer for the insightful feedback. We will add more intuition about the dynamics in the model description and compare ORGaNICs to LSTM, an RNN architecture similar to ORGaNICs. We will also provide intuitive definitions of the terms used in the stability proof.\n\n\n- “in L43, it is written that stability … Zucchet and Orvieto, 2024”. **Answer:** We appreciate the referee's observation. While stability indeed addresses the issue of exploding gradients, excessive stability—where the real parts of the eigenvalues of the Jacobians are significantly smaller than 0—can lead to a lossy system prone to vanishing gradients. ORGaNICs effectively mitigates both issues: Exploding gradients: through its inherent stability. Vanishing gradients: by processing information across various timescales (and via the long effective time constant when the membrane time constants are fixed) while maintaining stability, resulting in a blend of lossy and non-lossy neurons. The efficacy of ORGaNICs in mitigating vanishing gradients is demonstrated by its competitive performance against architectures specifically designed to address this issue, such as LSTMs. We will provide a more detailed explanation of this in our revised manuscript. Moreover, we posit that the built-in normalization feature of ORGaNICs may alleviate the \"curse of memory\" described by Zucchet and Orvieto (2024), as normalization is proposed as a potential solution to this problem.\n\n\n- “in L61, … stability of the system”. **Answer:** Divisive normalization (DN) was introduced as a model of steady-state response of neurons, functioning as a static nonlinearity similar to batch and layer normalization. However, in the brain, there are no static nonlinearities; it is proposed that DN is achieved via a recurrent circuit. ORGaNICs is such a recurrent circuit designed so that the responses of neurons at steady state follow the DN equation, with stability being an emergent property of the circuit. The DN equation indeed generalizes batch and layer normalization, as shown in the work by Ren, M., Liao, R., Urtasun, R., Sinz, F. H., & Zemel, R. S. (2016).\n\n\n- “in L349, … made more precise”. **Answer:** The key similarity is that both operate via changing the input gain, which also models a wide range of results from neurophysiological and psychophysical experiments. In particular, it has been shown that normalization explains experimental measurements of temporal attention (Denison, Carrasco, Heeger, Nature Human Behavior, 2021), i.e., over the time axis, analogous to transformers.\n\n- “the kind of neural network … additions to Section 6.1”. **Answer:** We thank the reviewer for pointing this out. We checked if the models dynamics are stable throughout training and found that this is indeed true. Please see Fig.2 of the attached pdf.\n\n- “in Section 6.2, … of the performance”. **Answer:** We are aware that state-space models perform better at these tasks and we will mention this (along with citations) in the paper, but our goal is to compare ORGaNICs with other RNNs designed with the property of less-dissipation and stability to mitigate the problem of vanishing and exploding gradients. We thank the reviewer for pointing out the sensitivity of initialization on sequential modeling tasks. We have added that detail (uniform random initialization) in the revised manuscript.\n\n**Questions:**\n\n\n- “The model … E neurons?” **Answer:** Experimental data suggest that the ratio of E and I neurons is 4:1. Our results are also applicable for a different E/I ratio than 1.", "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:
We thank all the reviewers for dedicating their time and providing valuable feedback on our submission. Your insights and comments are greatly appreciated and will help improve the quality of our work.

In response to the reviewers' comments, we have conducted additional experiments and analyses. Please find the following new results in the attached PDF:
- $\textbf{Table 1:}$ Based on the suggestion from Reviewer 1, we demonstrate that ORGaNICs can be trained with fixed time constants on both sequential and permuted MNIST datasets without requiring gradient clipping or scaling.
- $\textbf{Figure 1:}$ In response to Reviewer 1's inquiry about ORGaNICs' stability when the recurrent weight matrix defined as $\mathbf{W}_r=\alpha \mathbf{I}$ for $\alpha \in [0,1]$, we simulated 10,000 networks with random weights, inputs and values of $\alpha$, and with random initial values of $\mathbf{y}$ and $\mathbf{a}$. We found that all simulations converge to the fixed point identified by the iterative algorithm. These fixed points are stable, as the largest real part of the eigenvalues, for all the networks, is strictly negative.
- $\textbf{Figure 2:}$ Based on the insight from Reviewer 5 that ORGaNICs, when trained on static input classification task (Section 6.1), acts as a deep equilibrium model (DEQ) (Bai, S., Kolter, J. Z., & Koltun, V. (2019)). DEQs are known to be prone to instabilities during training. We demonstrate that since ORGaNICs is intrinsically stable for all parameters and input in the conditions assumed in Section 6.1, it stays stable throughout the training. Specifically, the fixed points remain attractive as all of the eigenvalues, across all the test samples, have negative real parts throughout the training.

Meta-review: This paper proves the unconditional stability of the ORGaNICs model, a biologically plausible RNN that implements divisive normalization. The authors provide rigorous mathematical analysis, particularly under specific conditions like the identity recurrent weight matrix, and demonstrate empirical stability across various tasks, including static and sequential classification benchmarks. Reviewers raised concerns about the broader applicability of the stability proof and the reliance on empirical evidence for more general cases. However, the authors addressed these concerns with additional experiments. While all reviewers voted for acceptance, the paper could benefit from clearer exposition and more comprehensive comparisons.
Decision: Accept (poster)