Conductance-based reversal potentials in spiking recurrent neural networks enhance energy efficiency and task performance

1 Introduction

Spiking recurrent neural networks (SRNNs) are a direct analog of biological neural networks (Maass, 1997) and can achieve relatively good performance on complex tasks as compared to classic RNN architectures (gated recurrent units, long short-term memory) or feedforward artificial neural networks (ANNs) (Bellec et al., 2018; Yin et al., 2021; Fang et al., 2021). SRNNs also offer energy-efficient solutions for neuromorphic computing (Roy et al., 2019). The performance of SRNNs has been improved due to several developments, including: the use of surrogate gradients to implement back-propagation-through-time (BPTT) (Werbos, 2002); the use of spiking neurons with adaptation mechanisms (Bellec et al., 2018); and increasing or optimizing the heterogeneity across neurons (Perez-Nieves et al., 2021; Bittar and Garner, 2022). However, many properties of biological neurons have yet to be implemented in SRNNs, and the ways in which these properties affect their performance and energy efficiency remain partially unexplored.

Several studies indicate that the performance of SRNNs can be significantly boosted by using neurons with adaptive thresholds or currents, as compared to the standard LIF models (Bellec et al., 2020, 2018; Baronig et al., 2025; Bittar and Garner, 2022). Networks with adaptive Linear Integrate Fire neurons (adapLIFs) are able to reproduce both subthreshold behaviors of so-called excitability class 1 and class 2 neurons (Baronig et al., 2025; Hodgkin, 1948; Izhikevich, 2000), while LIF neurons are restricted to class 1 neurons. Furthermore, these neurons exhibit resonant properties that give rise to subthreshold oscillatory dynamics, and they exhibit longer memory time constants, allowing them to solve tasks with greater memory requirements. The addition of such neurons can be seen as an example of boosting performance by introducing heterogeneity (Perez-Nieves et al., 2021).

However, the implementation of adapLIF neurons in a task-solving SRNN (Bittar and Garner, 2022; Baronig et al., 2025) still leads to voltage dynamics with major fluctuations that are not biologically plausible. A recent study suggests that clipping a neuron's voltage distribution can improve the performance of a feedforward SNN (Guo et al., 2022). In biological neurons, voltage dynamics are constrained by specific mechanisms, namely the existence of reversal potentials. Reversal potentials give rise to voltage-dependent synaptic currents and can, in fact, reverse the sign of GABAergic (pertaining to the neurotransmitter gamma-aminobutyric acid [GABA]) and glutamatergic synaptic inputs. For example, GABAergic inputs become excitatory when the receiving neuron's voltage is below the chloride reversal potential. Here, we define a new neuron model in which GABAergic and glutamatergic neurons give rise to recurrent synaptic inputs that are governed by a reversal potential. The SRNN thus contains two classes of neurons, each with positive and negative reversal potentials, mimicking Dale's law (i.e., neurons that use a single neurotransmitter), and is either inhibitory or excitatory depending on the receiving neuron's membrane potential.

We performed simulations on the benchmark Spiking Heidelberg Digits (SHD) (Cramer et al., 2022) classification task and evaluated our biologically inspired model, comparing it to the state-of-the-art neuron model adapLIF. We provide a theoretical analysis of the model's dynamics and an analysis of the subthreshold resonance of adapLIF neurons. We show that introducing a reversal potential boosts the SRNN's task performance and energy efficiency.

2 Results

Our starting point is the well-known adapLIF neuron model (Baronig et al., 2025; Bittar and Garner, 2022; Bellec et al., 2018), which we defined in our baseline model. We used a symplectic-Euler (SE) discretization adapted from Baronig et al. (2025), which is described in “Methods” section,

ht=αht-1+(1-α)(-wt-1+It)-VthresholdSt-1St=H(ht-Vthreshold)wt=βwt-1+aht+bSTIjt=∑inWjiinxit+∑j≠iWjiSit-1

Note that using the bounds on the parameters α, β, aandb as discussed in Bittar and Garner (2022), the trained SRNN contained neurons whose membrane potentials were hypersensitive to the first external input value and diverge to either positive or negative infinities. We prevented this effect by using the bound a≥0 as discussed in Baronig et al. (2025); note that this bound is also required for neurons to develop subthreshold resonance (see “Methods” section). We compared the baseline model with two other models in which aspects of biological networks/neurons were introduced:

• Dale model, which uses the same equations above but with a restriction on the recurrence matrices to obey Dale's law;
• AM+GA models, which use Dale's law and in addition implement the voltage reversal potentials, such that synaptic inputs are modulated based on the post-synaptic neuron's membrane potential (Figure 1).

The voltage reversal potentials were mathematically described by changing the input current, *I*, in the baseline model to:

Ijt=∑inWjiinxit+∑i≠jnexcWji(θexc-hjt-1)Sit-1 +∑i=nexc+1≠jninhWji(hjt-1-θinh)Sit-1

Here, θ_*exc* and θ_*inh* are the reversal potentials analogous to glutamatergic (AMPA) and GABAergic inputs, respectively. To determine the relative contributions of conductance-like gain that these “glutamatergic” and “GABAergic” reversal potentials cause, we further defined two additional models: AMPA, GABA, where we implemented the voltage reversal potential either for glutamatergic or GABAergic inputs (i.e., only the green/red line in Figure 1). We have set these threshold values at −1 and +2 as biological estimates based on the AMPA and GABA reversal potentials at positive and negative voltages (Rogers and Rogers, 1987); for the fixed versions of the models the values were kept constant whilst for the trainable versions these values were learnable parameters in the network, and equal inside each layer.

Next, we tested the models on the standard Spiking Heidelberg Digits (SHD) classification task. The spiking RNNs were trained offline with supervised learning for two different-sized networks: Architecture 1 (Arch-1L), which had one hidden layer (and a readout layer); and Architecture 2 (Arch-2L), which contained two hidden layers (and a readout layer). The SRNNs were trained via BPTT (see “Methods” section). We evaluated task performance on the test set through accuracy; and energy consumption through two measures, which describe practically useful energy measures when implementing these networks: Scaled Median Absolute Deviation (SMAD) and Spike energy consumption. These energy measures estimate the average energy consumption per trial by weighing the mean network activity, st¯, and the mean synaptic transmission of the network caused by said spiking activity, τt¯. This energy measure reflects the fact that energy consumption in biological neural networks is both driven by the spike itself as well as the synaptic activity caused by the spiking (Levy and Baxter, 1996; Lennie, 2003); note that we are not concerned here with energy consumption in neuromorphic chips. In more detail, the latter energy measure can be described for a particular sample m and time point t:

Etm=c1stm¯+c2τtm¯Et¯=∑m=1Mc1stm¯+c2τtm¯Etotal=∑tEt¯

where

τtm=|stm|⊗|W||W|=∑j|Wji|

Here, c₁, c₂ are constants that weigh the impact of activity and synaptic transmission (following Ali et al. (2022) c1=13,c2=23). We can inversely correlate this measure with network sparseness and thus with the network's energy efficiency (see “Methods” section for details).

Table 1 summarizes the accuracy and energy efficiency of the different models.

Table 1. Results for different models for two different size networks.

The number of layers on the left refers to the number of hidden layers in the trained networks. Accuracy and Voltage fluctuation relate to performance on the test set averaged across 5 iterations of the model; Voltage fluctuation is measured using the SMAD. The train (trainable) models had learnable reversal thresholds whilst the fixed models kept the values at our biological estimates.

In comparison to the baseline model, the Dale model showed a reduction in spike energy and voltage fluctuations for both the Arch-1L and Arch-2L. Introducing the voltage reversal potential led to a substantial and systematic reduction in both energy consumption and voltage ranges compared to the baseline and Dale models. Also, the differences in the mean value between the fixed and their dynamic counterparts were non-significant across the models (except for the voltage fluctuations in the AM+GA and GABA models).

In Arch-1L, we observed an approximately 3-fold reduction in spiking energy, and a 1.5-fold reduction in the SMAD of the membrane potential fluctuations for the AM+GA models as compared to the baseline and Dale's model. In Arch-2L, we observed similar improvements in energy efficiency. Thus, neurons fired more sparsely (see Tables 2, 3 and Supplementary Figures S1, S2) and had a narrower distribution of voltages in the model with reversal potentials. Analysis of the separate effects of AMPA and GABA suggests that both tend to reduce spiking energy and voltage fluctuations, but act additively when combined (Table 1). To illustrate these effects, we show some of the neuron responses (Figure 2) for the baseline network with large membrane potential fluctuations, compared to sparser firing and constrained membrane potentials with the voltage reversal potential for the AM+GA model.

Table 2. Parameters across neurons for different models after training.

For two-layer architectures, values are averaged across layers. Values are mean ± std across iterations. For the train AM+GA and fixed AM+GA models, |λ|, frequency, and fraction oscillatory were computed from the conductance-dependent dynamics.

Table 3. Parameters for 2-layer architectures split per layer.

Values are mean ± std across iterations. For the train AM+GA and fixed AM+GA models, |λ|, frequency and fraction oscillatory were computed from the conductance-dependent dynamics.

Next, we examined the effect on task performance.

Introducing the Dale model for Arch-1L significantly reduced the accuracy as compared to the baseline model, with a reduction of close to 1.5%, typical of introducing Dalean constraints on neural networks. For Arch-2L, the accuracy did not significantly change from the baseline to the Dale model.

However, the accuracy drop in the Arch-1L Dale network as compared to the base model was recovered by introducing the reversal potentials; strikingly for Arch-2L there was a significant increase in accuracy using the AM+GA models compared to the baseline and Dale models.

Analysis of the performance suggests that including either an AMPA or GABA reversal potential is enough to get accuracy close to the full train AM+GA model; however, best performance and energy reduction are obtained by including both.

Notably, even though there is no significant difference in the mean values of the fixed/train/baseline(for Arch-1L), the std between iterations is significantly smaller across all measures for the train AM+GA model, suggesting more consistent performance across model iterations.

The learned values for the reversal thresholds are aligned with our biological priors and are consistent between AM+GA models and the individual AMPA/GABA models (Table 4).

Table 4. Excitatory and inhibitory thresholds θ across architectures and model types.

Rows with N/A in θ_inh correspond to AMPA-only models; rows with N/A in θ_exc correspond to GABA-only models; rows with both values correspond to AM+GA models.

2.1 Analysis of unit properties

To understand how reversal potentials affect neuronal activity in the network, we conducted further theoretical analyses of adapLIF neurons with reversal potentials. As previously shown (Baronig et al., 2025), adapLIF neurons can exhibit resonances at different frequencies, a result of their being a system of two coupled differential equations that can exhibit subthreshold resonance. The performance gain observed with adapLIF neurons likely relates to increased heterogeneity and time constants (Bellec et al., 2018). However, subthreshold resonance may also play a role because RNNs with units configured as damped harmonic oscillators coupled oscillator RNN (CoRNNs) have potential theoretical benefits and enhance performance on various tasks (Rusch and Mishra, 2021), in particular when units show resonance at the frequency of the external input (Effenberger et al., 2025).

By analyzing neural properties (see “Methods” section and Supplementary Figures S1, S2), we found that neurons exhibited a broad distribution of oscillatory frequencies in the range 0–16 Hz (see Supplementary Figures S3, S4), suggesting that the network learns substantially heterogeneous resonance responses. The AM+GA models with reversal potentials had substantially fewer oscillatory units than the Dale's model and the baseline model, with lower resonance frequencies (Table 2). Thus, the improved energy efficiency and performance were not associated with increased oscillatory activity but rather with greater heterogeneity between oscillatory and non-oscillatory activity.

Furthermore, in Arch-2L, the AM+GA models also had spike patterns closer to Poisson processes (CV closer to 1, Tables 2, 3) which translates to spike train statistics that match biological neurons (Dayan and Abbott, 2005). Notably, the AM+GA model exhibited more adaptation in response to subthreshold activity, with larger values of a.

It can be shown that the inclusion of a reversal potential makes the neuron's properties dynamic and dependent on the time-varying input, rather than static. To see this, we note that our system of equations is written as

ht=αht-1+(1-α)(-wt-1+It)-VthresholdSt-1wt=βwt-1+aht+bStIjt=∑inWjiinxit+∑i≠jnexcWji(θexc-hjt-1)Sit-1 +∑i=nexc+1≠jninhWji(hjt-1-θinh)Sit-1

If we consider only one neuron j and it is under subthreshold dynamics, we can simplify these equations into:

ht=αht-1+(1-α)(-wt-1+It)wt=βwt-1+aht+bStIt=Iind-J(St-1,W)ht-1

Here, Iind=∑inWjiinxit+∑i≠jnexcWjiθexcSit-1-∑i=nexc+1≠jninhWjiθinhSit-1 has the feedforward component and the recurrent component independent of h_*j* and J(St-1,W)=∑i≠jnexcWjiSit-1-∑i=nexc+1≠jninhWjiSit-1 represents the recurrent gain dependent on the membrane potential. Note that J > 0 and that J represents the sum of the magnitudes of excitatory and inhibitory conductances. We can then write the state transition matrix by plugging in h_*t*→w_*t*:

Here, we can see that the inputs themselves determine the state transition matrix A, making the response property of the neuron input-dependent, where

A=[α-(1-α)J -(1-α)a(α-(1-α)J) β-a(1-α)]

To understand the impact of J (the amount of conductance), we perform an eigenvalue analysis in the subtreshold regime (see the “Methods” section for exact derivations). This analysis shows oscillation magnitudes and frequencies that are dependent on recurrent network activity itself, as it determines J as a function of time for each neuron.

The expression for the eigenvalue magnitudes equals (see “Methods” section)

Thus, increased conductance J strictly decreases the magnitude of subthreshold oscillations, because β>0 and ᾱ>0. This reduction in amplitude is related to the classic shunting effect of synaptic conductance in biophysics (Koch, 2004).

The eigenvalues are complex as long as

ΔA=Tr(A)2-4det(A)<0(2)

The discriminant is an upward parabola as a function of J, reaching a minimum at J^*; that is, the system is most likely to oscillate at this point. Oscillatory behavior requires J to lie in the interval {J₋, J₊}, where

J±=J*±2βaᾱ(3)

and

J*=α-β-aᾱᾱ.(4)

We find that all of the neurons in the trained model have J^* < 0.

Hence, increased conductance J invariantly decreases the oscillation magnitude; the eigenvalue magnitudes keep decreasing until the system is no longer oscillatory.

Thus, the best performing model train AM+GA, which has more oscillatory units for J = 0 (see Supplementary Tables S4, S5), shows that the tendency to oscillate decreases dramatically with increased J. Note that it is required that a ≥ 0 for the system to exhibit oscillations (see “Methods” section). As a → 0, the window of conductances J in which oscillations occur becomes increasingly narrow.

Finally, for complex discrete eigenvalues we have a frequency f=ω2πfs with ωc=arccos(12Tr(A)detA) (note that deriving the frequency from continuous dynamics as done in other study is inaccurate). We obtain that the maximum frequency is obtained for

Jc=α+ᾱa-βᾱ.

As a consequence, the frequencies may initially increase slightly as a function of J up to J_*c*, then decrease until the system stops oscillating at J⁺, at which point the frequencies are zero. For our range of parameters, on average J_*c* ≤ 0 (see Table 2), such that frequencies will on average decrease with increased conductance J. In fact, for AM+GA model we compute the mean neuron parameters with our dynamic values of J and averaged them for Arch-1L and Arch-2L, in Tables 2, 3, which yield significantly lower results when compared to the static (no J taken into consideration) results (see Supplementary Table S4).

In sum, the AM+GA models with reversal potentials have substantially fewer oscillatory units with weaker oscillations. The reversal potential further leads to weaker oscillations and a non-linear dependence of frequency on the conductance J.

3 Discussion

The performance of SRNNs has been improved due to the introduction of several features, including the use of surrogate gradients to implement back-propagation-through-time (BPTT) (Werbos, 2002); units showing adaptation mechanisms (Bellec et al., 2018); increased heterogeneity of single-neuron properties (Perez-Nieves et al., 2021; Bittar and Garner, 2022). However, many properties of biological neurons have yet to be implemented in SRNNs, and the effects of these properties on their performance and energy efficiency remain partially unexplored. The implementation of this adapLIF in a task-solving RSNN (Bittar and Garner, 2022; Baronig et al., 2025) still results in non-biologically plausible voltage dynamics, as subthreshold values fluctuate outside the biological range. In biological neurons, voltage dynamics are constrained by specific mechanisms, namely the existence of reversal potentials (Koch, 2004). Here, we define a new neuron model in which GABAergic and glutamatergic neurons give rise to recurrent synaptic inputs that are governed by a reversal potential. The SRNN thus contains two classes of neurons that mimic Dale's law (i.e., neurons that use a single neurotransmitter), which are inhibitory or excitatory depending on the receiving neuron's membrane potential.

We trained SRNNs on the Spiking Heidelberg Digits and show that SRNNs with reversal potentials cut spike energy by up to 3 × , while maintaining or increasing task accuracy. We provide a theoretical analysis of the model's dynamics and examine the subthreshold resonance of adapLIF neurons. Improved model performance with reversal potential is concomitant with a reduction in the number of oscillatory units and a decrease in oscillation frequency compared with the baseline model. The reversal potential makes subthreshold resonance conductance dependent, such that increased conductance reduces oscillatory behavior and decreases frequency.

Further study is needed to precisely identify the reasons why including the reversal potential increases task performance and energy efficiency. There may be several reasons: (1) The reversal potential naturally constrains the voltage range of the neurons, which is explained by the fact that the sign of excitation and inhibition changes around the reversal potential, thereby improving the energy efficiency. This may improve task performance by keeping voltages within the range where the surrogate gradients do not vanish, thereby mitigating the vanishing-gradient problem. In agreement, a previous study suggests that clipping a neuron's voltage distribution can improve the performance of a feedforward SNN (Guo et al., 2022). (2) Mathematically, the reversal potential allows, for example, inhibitory input to have divisive effects on excitatory inputs rather than just subtractive effects, thus enhancing the computational capacity of the network (Koch, 2004). That is, the reversal potential adds a nonlinear integration of inputs inside the neuron, which may enhance the computational capacity of neurons. (3) We showed that the reversal potential allows for the resonance frequencies of the individual neurons to become input-dependent. Hence, a neuron's resonance function may change in response to input, thereby creating more complex, time-dependent feature-selective responses tuned to specific input patterns.

Reversal potentials are, biologically, an unavoidable consequence of the electrochemical gradient of ion channels defined by intracellular and extracellular concentrations. For example, GABA is inhibitory because the intracellular chloride concentration is much lower than the outside of the cell, such that the chemical force is stronger than the electric force when the voltage is above the reversal potential (Koch, 2004). In fact, this is a property that depends on the KCC2 chloride receptor and changes through development, but can also vary between areas (e.g., dentate gyrus) (Ben-Ari, 2007; Sauer et al., 2012). Thus, it stands to reason that the precise GABAergic reversal is a property that is evolutionarily optimized. In this study, we showed both a fixed biological estimate approach, as well as optimal value convergence, which brought more consistent results across different weight initializations and can be an indicator that biological tuning around these thresholds brings consistent energy efficiency/performance across individual areas.

In this study, we report an increase in performance by incorporating Dale's law with reversal potentials; however, there are other methods that may improve task performance that we did not examine here. Baronig et al. (2025) includes a biologically unrealistic q parameter in the adaptation equation (see “Methods” section) that needs to be adjusted for each architecture. They also perform additional preprocessing on the inputs, strongly reducing the number of input channels.

Recently, there has been some interest in the ability of RNNs with oscillatory units (coupled oscillator RNN) to improve task performance (Rusch and Mishra, 2021; Keller and Welling, 2023; Effenberger et al., 2025; Baronig et al., 2025; Lanthaler et al., 2023). The subthreshold dynamics of adLIF neurons can produce oscillatory dynamics (Baronig et al., 2025), with subthreshold resonance around relatively low frequencies around 5–8 Hz (i.e., in the delta range). (Note that we improved the frequency calculation of these units, taking into account the discrete nature of the dynamics.) In general, it is unclear whether the tendency to oscillate or rather the introduction of an additional current/self-influence accounts for possible increases in performance in RNNs with oscillatory units compared to RNNs without oscillatory units (Rusch et al., 2021); noting that a damped harmonic oscillator can be broken up in two first order equations, one with a fast and one with a slow component that allows learning long temporal dependencies. Interestingly, we showed that reversal potentials strongly reduced the fraction of oscillatory units (e.g., from 0.85 to 0.26 in a one-layer neural network) and the magnitudes of the eigenvalues, even though task performance and energy efficiency improved compared to the Dale network. The fact that better performance was obtained with far fewer oscillatory units but similar time constants suggests that the propensity to oscillate may not, per se, boost task performance or energy efficiency. In the networks trained here, the propensity to oscillate was strongly reduced by the input stimulus, leading to a marked increase in conductance. Interestingly, this shows some similarity to biological neurons and networks of neurons, which tend to show more low-frequency activity in the absence of sensory stimulation (i.e., during rest), while stimulus-evoked activity classically shows desynchronization of neural responses (Buzsáki, 2006).

In sum, the findings here are another step toward greater biological realism in spiking neural networks, enabling efficient, high-performance training with Dale's law. These findings likely have important applications for neuromorphic computing technologies, in which energy efficiency and network sparseness play an important role.

4 Methods

4.1 Neuron model

We ran a discretized version of the adapLIF model, using the symplectic-Euler (SE) discretization (Baronig et al., 2025). Some small changes to the dynamics were made compared to those in the studies by Bittar and Garner (2022); Baronig et al. (2025), as we will discuss below. The discrete update equations according to the SE-discretization method are:

ĥt=αht-1+(1-α)(-wt-1+It)(5)

St=H(ĥt-Vthreshold)(6)

ht=ĥt(1-St)(7)

wt=βwt-1+(1-β)(aht+bSt)(8)

with h the hidden membrane potential, S_*t* a binary variable for the spike output of the neuron, w the adaptation variable, and I the input currents to the neuron at that time step. The neuron parameters α, β, a, and b correspond to the membrane decay rate, adaptation decay rate, and sub- and supra-threshold gain constants, respectively; these parameters remain fixed post-learning. We initialize them over a biologically plausible range based on studies by Bittar and Garner (2022); Baronig et al. (2025) and allow them to train individually per neuron to add heterogeneity to the network (Perez-Nieves et al., 2021).

The discrete update equations we used for the simulations are:

ht=αht-1+(1-α)(-wt-1+It)-VthresholdSt-1(9)

St=H(ht-Vthreshold)(10)

wt=βwt-1+aht+bSt(11)

4.1.1 Differences from previous works

There are two major differences from Equations 5–11. First, we avoid the hard reset, Equation 7 used in the study by Baronig et al. (2025), as this represented in preliminary tests a 7% drop off in accuracy, and use instead a soft reset similar to the one used in the study by Bittar and Garner (2022), by adding “−V_*threshold**S_t*−1” in Equation 9.

Second, we drop the (1−β) term in Equation 8 as in preliminary tests the network did not learn efficiently: accuracy peaked [55%, 60%] [compared to (85%, 95%) we get otherwise]. This is due to (1−β) being very small (in our range of values), effectively removing the dynamical component of our adaptation variable. The scaling of (1−β) was solved in the study by Bittar and Garner (2022) identically to us (seen in the source code, as in the original paper they have that scaling factor present); and in the study by Baronig et al. (2025) by adding another scaling factor q, reshaping Equation 8 to:

wt=βwt-1+q(1-β)(aht+bSt)(12)

We chose to remove the term from our main results since Baronig et al. (2025) suggests q needs to be tuned for each task as a global gain on the neuron's adaptation constants. In fact, removing this factor can be seen as absorbing the (1−β) factor into the arbitrary parameters a and b, which rescales the values of the initialization intervals for a and b, and whose bounds we have discussed previously, and in more detail below.

4.1.2 Additions to the neuron model

To the baseline equations, Equations 9–11, we added two biologically realistic features: (1) Dale's law; and (2) a voltage reversal potential.

(1) Dale's law restricts the communication of a neuron to a single neurotransmitter for all postsynaptic connections. This means, for our computation, that we divide our nodes into two subpopulations, with their recurrent weights either positive (excitatory neurons) or negative (inhibitory neurons).

The two populations comprise 50% of the population in this case, implemented to introduce an inductive bias toward balanced E/I.

(2) By using a reversal potential, dendritic inputs to neurons are scaled by their local voltage potential. Its existence is well established in the neuroscientific literature (Brown et al., 1971; Mofidi et al., 2020; Eisenberg et al., 2015; Koch, 2004).

To model the reversal potential, we changed our input current equation to:

Ijt=∑inWjiinxit+∑i≠jnexcWjirecurrent(θexc-hj)Sit-1+∑i=nexc+1≠jninhWjirecurrent(hj-θinh)Sit-1

4.2 Theoretical derivations of the neuron model

The system of equations in the subthreshold regime is given by

ht=αht-1+ᾱ(-wt-1+It)wt=βwt-1+aht+bStIt=Iind-J(St-1,W)ht-1

where ᾱ=(1-α). Here, Iind=∑inWjiinxit+∑i≠jnexcWjiθexcSit-1-∑i=nexc+1≠jninhWjiθinhSit-1, is the feedforward component plus the recurrent component independent of h_*j*.

J(St-1,W)=∑i≠jnexcWjiSit-1-∑i=nexc+1≠jninhWjiSit-1, is the recurrent gain dependent on the membrane potential.

Setting S_*t* = 0 (subthreshold regime), the system of equations can be written as:

ht=(α-ᾱJ)ht-1-ᾱwt-1+ᾱIindwt=a(α-ᾱJ)ht-1+(β-aᾱ)wt-1+aᾱIind

We can then define the subtreshold state transition matrix A,

A=[α-ᾱJ-ᾱa(α-ᾱJ)β-aᾱ]

We have

λ=12Tr(A)±12Tr(A)2-4det(A)(13)

and

Tr(A)=α+β-ᾱ(a+J)(14)

det(A)=β[α-ᾱJ](15)

4.2.1 Eigenvalue magnitude

We first consider the case that there is a pair of complex eigenvalues, that is, the neuron exhibits subthreshold resonance (which occurs in more than half of the measured neurons). In this case, we have the eigenvalue magnitude

Hence, increased recurrent input conductance J tends to decrease the magnitude of subthreshold oscillations as β > 0 and ᾱ>0.

4.2.2 Effect of input on tendency to oscillate

The eigenvalues are complex as long as

ΔA=Tr(A)2-4det(A)<0(17)

The discriminant is an upward parabola, reaching a minimum at J^*, that is, the system is most likely to oscillate at this point, with complex eigenvalues in the interval

J±=J*±2βaᾱ(18)

where

J*=α-β-aᾱᾱ.(19)

We can see that J^* < 0 typically over our range of parameters (Table 4). As the conductance J increases, there is a conductance value J⁺ at which point the neuron stops oscillating. Hence, an increase in J leads to a decrease in eigenvalue magnitudes until the neuron falls outside the range of oscillatory behavior.

Increased adaptation a widens the window in which oscillations occur; increased membrane potential time constant (α) also increases the window. If a > 0, there exists an interval J ∈ (J₋, J₊) in which the eigenvalues are complex, since ΔA(J)=ᾱ2(J-J*)2-4βaᾱ. In the absence of an input, J = 0, we can see that Δ_*A* < 0 if

a-<a<a+, a-,a+>0(20)

Here

a±=α+β±2αβᾱ.(21)

4.2.3 Bounds on the stability of the system

Stability is given by the Schur inequalities.

1-TrA+detA>0⇒aᾱ>-[ᾱ-β+αβ+ᾱ(1-β)J].(22)

amin(J)=-1-α-β+αβ1-α-(1-β)J=-(1-β)(1+J)(23)

If a > 0 and J > 0 then a > a_min(J). This motivates the bound a > 0 to ensure stability.

We have from |detA| < 1

βα-1βᾱ<J<βα+1βᾱ.(24)

Finally we have

amax(J)=1+α+β+αβᾱ-(1+β)J, a<amax(J).(25)

For our ranges of values, this implies that J_*max* is very large (>10) and thus the system, given our parameters, is stable. It is further evident that as the network grows larger in size, the network weights should be scaled as 1N to prevent J from growing too large.

4.2.3.1 J's critical point

For complex discrete eigenvalues, we have a frequency f=ω2πfs. We will derive the dependency of ω on J. We shall use the correct expression for discrete-time dynamics (note that the continuous-time expression would not necessarily be accurate). We define the numerator/denominator of ωc=arccos(12Tr(A)detA) as g, p, and the full argument of the arccos as f. We are going to be searching for the critical point with respect to J.

∂Jωc=∂Jarccos(gp︸f)=-∂Jf1-f2

Expanding the terms and simplifying, we obtain:

∂Jωc=ddJarccos(f(J))=-f′(J)1-f(J)2

Hence, the critical point occurs when f′(J) = 0. We have

f′(J)=-ᾱβ(α-β+ᾱ(a-J))4(β(α-ᾱJ))3/2=-ᾱβ(α-β+ᾱ(a-J))4p3

Thus, the critical point occurs when Jc=α+ᾱa-βᾱ.

For the second derivative evaluated at the critical point, we have

∂J2ωc|Jc=-f″(Jc)1-f(Jc)2

Note that

f″(Jc)=ᾱ24β(β-ᾱa)3/2=ᾱ2β4p(Jc)3, p(Jc)=β(β-ᾱa)

Thus, f″(Jc)>0, and the critical point J_*c* is always a maximum.

Note that we require that Jc≤J+, that is, we require that it lies in the domain of complex eigenvalues. Thus, if J_*c* > 0, initially frequencies increase as a function of J, until J = J_*c*, after which the frequency decreases; if J_*c* < 0, frequencies always decrease. At J⁺, we have Tr(A)² = 4det(A) and hence ω_*c* = 0. Hence, frequencies decrease until ω_*c* = 0 after which subthreshold oscillations disappear.

4.3 Data

We trained networks on the Spiking Heidelberg Digits (Cramer et al., 2022). This dataset consists of close to 10,000 aligned audio recordings of spoken digits (0–9) in English and German from 12 speakers, with two exclusive ones to the test set. These recordings are labeled according to digit and language, yielding a total of 20 classes. We generate the neural network inputs by binning the spike times into 4 ms intervals going from 0 until 1 s, generating a 250 step long time series. There are 700 input channels, the values of which are integers representing the number of times that a neuron fired in that time bin.

4.4 Task and neural network architectures

We examined the performance of different neural networks in learning to classify the SHD data in a supervised setting. We used a connected internal feedback spiking recurrent neural network (SRNN) with a final readout layer, similar to the one discussed in Bittar and Garner (2022). We used two different architectures, which had the same type of readout layer:

4.4.1 Architecture 1L

Here we used a two-layer deep network (1 hidden layer + readout layer) with 256 neurons in the hidden layer, with a total number on the order of 250k trainable parameters. The input currents for each neuron Ijt are going to be of the form:

Ijt=∑inWjiinxit+∑i≠jnWjirecurrentSit-1(26)

with xit the value for the ith input channel. In vector notation, we have

It=Winxt+WrecurrentSt-1(27)

where W_*in* is the 700x256 input weight matrix and W_*recurrent* is the 256x256 recurrent weight matrix with diagonal filled with 0s.

4.4.2 Architecture 2L

Here, we trained a three-layer deep network (two hidden layers + one readout layer) with 360 neurons in each hidden layer. The total number of trainable parameters on the order of 650k parameters.

We can again write the input current for each neuron Itl where l refers to the hidden layer of the network:

Itl=Winlxtl+WrecurrentlSt-1l(28)

with dim(Win2)=128x128 and xt2=St1.

4.4.3 Readout layer and loss function

The readout layer of our network consisted of a leaky integrator neuron (LI), achieved by setting a = b = 0 and removing the reset mechanism. From this output, we would compute the loss function defined in the study by Bittar and Garner (2022). For each time step, we apply a softmax function to the hidden potentials of the readout layer across neurons, then sum these normalized values across time steps for all neurons and compute a cross-entropy loss with the target class for each trial. For each batch, n, the loss was computed as:

oc=∑t=0Tsoftmax(hct)Ln=-∑c=019logexp(oc)∑i=019exp(oi)

For the reported results, we used the same regularization loss function as discussed in Bittar and Garner (2022). Additionally, for Arch-1L, we obtained similar performance using the regularization function discussed in (Zenke and Vogels, 2021) and computing oc=max(hct). The settings used for the Readout layer are compiled in Table 6.

4.4.4 Parameter initialization

We briefly explored the hyperparameter (α, β, a and b) space but since this space is quite large and we are focused on comparability, we decided to base it on previous work with similar models (Bittar and Garner, 2022; Baronig et al., 2025). The hyperparameters were initialized uniformly across the accepted range and, during learning, were clamped to their bounds to ensure realistic dynamics.

In the hidden layers, the feedforward weights were initialized uniformly across |w|<1fanin, and the recurrent weights were initialized using PyTorch's orthogonal method discussed in the study by (Saxe et al., 2013); as in state-of-the-art SRNN models (Baronig et al., 2025). For the readout layer, the weights were initialized from a normal distribution with a specified mean and standard deviation (SD).

When using Dale's law, we additionally split the recurrence matrix in half for the excitatory/inhibitory subpopulations. For the excitatory afferent weights, any negative initialized weights were swapped to their symmetric; for the inhibitory weights, the symmetric happened. When clamping was applied to this model, it was set to the minimum absolute value of the initialized weights.

4.4.5 Training details

The SHD does not have a predefined validation set, so we randomly split the training set into a training and validation set (80/20%), and all our results came from the test set.

We trained all the models for 200 epochs on the SHD training set. We trained each model for 5 iterations using the same initialization. Every setting range used was fixed for each architecture, and shared between architectures (Table 5).

Table 5. Hidden layer settings for both architectures.

Table 6. Readout layer settings.

Similar to the study by Baronig et al. (2025), we define α=e-Δtτm, and the membrane time constant as:

τm=τmmin+ωm(τmmax-τmmin)

The trainable parameter is ω_*m* ∈ [0, 1]. This is done so that all parameters are on the same order of magnitude for the ADAM optimizer. The analogous is done for the synaptic rate (β) and time constant (τ_*s*).

The different settings only came from the number of neurons, already specified above.

The setting names are slightly different from those used in the dictionary available in the code for readability.

For the surrogate gradient we used the SLAYER function with α = 5 and c = 0.4 for all layers (Shrestha and Orchard, 2018).

4.5 Analyses

For all analyses, we computed these values for each iteration of model X and then computed the mean and standard deviation (SD) across the five iterations. Per iteration, we ran through the full test set via our sampling function (in batches), computing the value and saving it into an array; we would perform this loop 10 times and then average the result as our final per-iteration value.

For Accuracy, we took y*=argmaxc∑tsoftmax(ht)c (c denotes the classification neuron that goes from 0 to 19, and y^* is our guess for the label) and compared it to the label, y, of the class for each batch.

For Test Spike Energy, we computed the average neuron energy consumption at each time point and summed them for the duration of the stimulus. Our measure was inspired by the study by Ali et al. (2022), Equations 7, 8. We defined it for a particular sample m and time point t in a batch:

Etm=c1stm¯+c2τtm¯(29)

Et¯=∑m=1Mc1stm¯+c2τtm¯(30)

Etotal=∑tEt¯(31)

where

τtm=|stm|⊗|W|(32)

The operator¯ is the average across neurons for a batch until Et¯. From this point on, it also includes the average across batches. Indeed, one can take the average across batches or neurons interchangeably, which is how we do it in practice. We chose this notation to maintain consistency with the original notation discussed in the study by Ali et al. (2022).

This metric has two components for each time step: 1) st¯, the mean network activity, and 2) τt¯, the mean network synaptic transmission caused by said spike activity. Both of these components were averaged across layers to ensure comparability between architectures. This measure is inversely correlated with the network's sparseness. Since the weight matrices did not differ dramatically between models, we can strongly correlate how sparse a network was with how low its energy value was (hence, how energy-efficient it was).

For Voltage fluctuations, we used the Scaled Median Absolute Deviation (SMAD) from 0 for the resulting membrane potentials of a batch:

SMAD=median(|h|)0.67

The sampling function has a built-in shuffling property. Given the way we loop through the test set, for each iteration, we compute a Monte Carlo estimate of the expected measures. However, due to the order of magnitude of the samples, our value is statistically valid.

For |λ|, Frequency and Fraction Oscillatory, from Tables 2, 3, in the AM+GA models we computed all the values of J per trial and used them to compute the respective |λ|, Frequency and Fraction Oscillatory at each time step where neurons had spiked, and hence J could be computed. We 1)then averaged over those time steps, 2)then across trials for each neuron, and 3)finally averaged the values of each neuron for the per-iteration “dynamic” value. We did these computations for each layer separately. We then did the weighted average of these “dynamic” values with the “static” values (that is, independent of J), where the weights were defined by the average number of time steps where the network was active/silent out of the total number of time steps, respectively. The final results are then the mean and standard deviation of the per-iteration comparison. For Table 2 we also averaged the per-layer result.

We used the conventional coefficient of variation (CV), where ISI is the interspike interval:

CV=ΔISIISI¯

4.6 Bound on a and reversal potentials

The bounds on the subthreshold adaptation constant, a, have a strong impact on the membrane dynamics Baronig et al., (2025). We bounded a ≥ 0 based on the study by (Baronig et al. 2025) and our own derivations (see “Methods” section). This regime, a ≥ 0, represents the adaptation regime, with higher voltages leading to greater adaptation.

The initial values for the reversal potentials were based on biological estimates, approximately matching reversal potentials in biological neurons, and were not subject to optimization. We trained 2 types of models using reversal potentials: fixed and train (trainable). As the name suggests, for the fixed models we maintained the reversal threshold at our biological estimate; for the train models we initialized the reversal thresholds at our biological estimates but treated them as learnable parameters of the network. For each layer, all neurons shared the same reversal threshold.

Generative AI statement

The author(s) declared that Generative AI was not used in the creation of this manuscript.

Any alternative text (alt text) provided alongside figures in this article has been generated by Frontiers with the support of artificial intelligence and reasonable efforts have been made to ensure accuracy, including review by the authors wherever possible. If you identify any issues, please contact us.

Publisher's note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Author disclaimer

The views and opinions expressed are however those of the author(s) only and do not necessarily reflect those of the European Union or the European Research Council. Neither the European Union nor the granting authority can be held responsible for them.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fncom.2026.1870593/full#supplementary-material

Supplementary material Supplementary material at the source

Supplementary material Supplementary material at the source

Supplementary material Supplementary material at the source

Supplementary material Supplementary material at the source

Supplementary material Supplementary material at the source

Footnotes

Edited by: Vignayanandam Ravindernath Muddapu, Azim Premji University, India

Reviewed by: Timothée Masquelier, Centre National de la Recherche Scientifique (CNRS), France

Syeda Shamaila Zareen, Shenzhen Technology University, China