Manuscript

Abstract

We compared reciprocal excitatory–inhibitory circuits with loop-off controls, then tested trained MNIST classifiers under activity constraints, changes to inhibitory timing and spike perturbations. Reciprocal feedback supported rhythmic low-rate activity, while trained classifiers retained useful accuracy at lower excitatory rates. At 25 Hz maximum pixel input, increasing stream presentation duration from 25 to 200 ms raised mean accuracy from 75.2% to 90.3%. The experiments do not isolate a causal advantage of gamma timing, measure energy use or test autonomous stream segmentation.

AI use disclosure

AI agents were used to draft, rewrite, restructure and condense this manuscript, including passages of scientific exposition. The lead author specified the intended meaning, directed revisions, corrected wording and interpretation, and repeatedly reviewed the manuscript, including paragraph-by-paragraph review. The lead author exercised editorial judgment over the retained text. Following review and approval by all authors, the final manuscript reflects the scientific account and interpretations they endorse, and all authors take responsibility for its accuracy, claims, citations and conclusions.

Results

Reciprocal coupling creates a gamma-rhythmic low-rate regime

We first asked whether reciprocal excitatory–inhibitory coupling could produce rhythmic population activity while maintaining low firing rates in individual excitatory neurons. We compared untrained conductance-based spiking networks driven by Poisson input, with reciprocal coupling between excitatory (E) and inhibitory (I) populations enabled (loop-on) or disabled (loop-off; Fig. 1A,B). Illustrative 400-ms responses used separately chosen Poisson drives of 5 Hz per channel for loop-off and 45 Hz for loop-on (Fig. 1C,D).

Loop-off excitatory spikes were dispersed, with silent inhibitory neurons (Fig. 1C); loop-on instead produced recurring E/I volleys and a 56.0-Hz spectral peak with higher-frequency harmonics, which loop-off lacked (Fig. 1D,E,G).

In matched 2–100-Hz drive sweeps on the same 1,024-input networks as the illustrative rasters, loop-off excitatory firing increased from 1.7 to 484.6 Hz while inhibition remained silent. Loop-on excitatory rates were constrained, staying between 1.7 and 11.6 Hz, while inhibitory firing reached 127.5 Hz (Fig. 1F,H).

We next varied the strengths of both directions of the excitatory–inhibitory loop to examine how coupling affected firing rates and rhythmic activity. Across an 11×11 reciprocal-coupling grid, absent E→I or I→E coupling left excitatory firing near 169 Hz and the rhythmicity metric (autocorrelation lobe–trough contrast; Appendix B2) near zero (Fig. 2A–C). With both E→I and I→E pathways present, stronger coupling reduced excitatory firing to single-digit rates, sustained inhibitory firing and increased the rhythmicity metric towards one across a broad region.

Selected rasters illustrate the same progression (Fig. 2D–F): the rhythmicity metric increased from 0.00030 without coupling to 0.268 at intermediate coupling and 0.989 at strong coupling, as excitatory firing became sparse and inhibitory volleys increasingly regular.

We used a separate four-variable mean-field model to examine how oscillations emerged as external drive increased (Appendix C). Leading complex-conjugate eigenvalues crossed zero real part at $I_{\text{ext}}^{\ast} = 0.594$ nA, corresponding to 27.6 Hz (Fig. 2G). Above onset, excitatory-rate amplitude increased continuously, with nearly coincident upward/downward sweeps and no resolved hysteresis (Fig. 2H). A positive amplitude-squared slope and $R^{2} = 0.999$ met the predefined numerical criteria for supercritical onset; criticality was not established analytically.

We then compared how slowing inhibition affected oscillation frequency in the mean-field model and the loop-on classifiers in inhibitory-decay training group TR03. Labels TR01–TR06 identify groups of training regimes, not individual networks or executions (Table 1). Increasing inhibitory decay from 4.5 to 27 ms reduced mean-field onset frequency from 30.2 to 17.9 Hz and median spiking-classifier spectral frequency from 67.3 to 12.2 Hz (Fig. 2I). The mean-field and spiking models showed the same qualitative behaviour—frequency decreased as inhibitory decay slowed. The comparison was intended to assess qualitative behaviour rather than establish quantitative agreement between mean-field onset eigenfrequencies and finite-drive spectral peaks.

The fixed reciprocal loop preserves accuracy at lower excitatory rates

We next trained loop-on and loop-off networks to classify handwritten digits from MNIST using backpropagation through time with surrogate gradients, asking whether reciprocal inhibition could reduce excitatory firing while preserving classification accuracy. In full-MNIST reference group TR01, loop-off and loop-on reached mean final validation accuracies of 95.3% and 91.4%, respectively, with dense loop-off firing and recurring loop-on volleys (Fig. 3A–C). In the separate reduced-pool group TR02, loop-on achieved 89.8% mean test accuracy at a mean excitatory firing rate of 16.6 Hz, versus 90.9% at 88.2 Hz for loop-off: a 5.3-fold firing-rate reduction and 1.1-percentage-point accuracy reduction (Fig. 3D). By adding a spike penalty with progressively lower firing-rate ceilings to both networks, we reduced excitatory firing and classification accuracy, but the networks responded differently. Loop-off firing fell throughout the sweep, accompanied by a steeper loss of accuracy, whereas loop-on firing decreased less at the tightest ceilings while retaining higher accuracy. At the 1-Hz ceiling, loop-on averaged 4.7 Hz and 81.3% accuracy, compared with 2.1 Hz and 63.0% for loop-off. Thus, equal training ceilings did not generally produce matched firing rates. At the 10-Hz ceiling, however, both averaged approximately 9 Hz, with higher accuracy for loop-on (88.6% versus 83.5%).

We next asked whether adding reciprocal inhibition after training could produce sparse rhythmic activity while preserving classification accuracy. Adding reciprocal inhibition to the unpenalised TR02 loop-off classifiers without retraining replaced dense excitatory activity with recurring E/I volleys (Fig. 4A,B). Across the coupling sweep, mean excitatory firing fell from 85.9 to 7.2 Hz (11.9-fold), inhibitory firing reached 40.0 Hz, and accuracy fell from 90.4% to 43.2% (Fig. 4C,D). The activity pattern therefore did not require learning with inhibition, but inserting the loop after training incurred a substantial classification cost.

Using TR05, we next asked whether the rhythmic activity observed with fixed reciprocal coupling persisted when the recurrent excitatory–inhibitory weights were trained alongside the input and readout weights. Training recurrent weights reduced the rhythmicity metric relative to fixed loop-on recurrence. The final rhythmicity metric was 0.999 with fixed recurrence, versus 0.094, 0.095 and 0.018 when recurrence was trained from standard, one-tenth-standard and zero initialization (Fig. 5C). These means describe three networks per condition responding to one fixed encoding of a reference digit, not test-set variability.

We also compared the recurrent connection weights before and after training. In all three networks in each condition, standard and one-tenth-standard training eliminated most E→I connections; most I→E connections remained positive and their mean strength across all entries increased (Fig. 5D–G). Mean E→I strength fell from standard but rose from one-tenth-standard initialization; zero-initialized recurrence remained zero. Low values of the rhythmicity metric accompanied asymmetric reorganization in these three-replicate comparisons, without identifying a causal weight change.

Excitatory firing is organised by population-cycle participation

We next trained the separate TR03 loop-on networks at each inhibitory decay time to examine how the timescale of inhibition shaped population rhythms, excitatory firing and classification accuracy. Across inhibitory decay times of 4.5–27 ms, mean spectral frequency fell from 67.5 to 11.7 Hz and excitatory firing from 18.3 to 2.8 Hz (Fig. 6A). The affine fit to six condition means had slope 0.285 Hz/Hz, intercept −0.70 Hz and $R^{2} = 0.997$. Mean test accuracy fell from 91.2% to 81.9% (Fig. 6B). Slower rhythms thus accompanied lower activity and a classification cost. The slowest condition lay below conventional gamma frequencies, so we describe the full sweep as population cycles rather than gamma cycles.[1]

To examine how individual neurons participated in these population rhythms, we counted each excitatory neuron’s spikes in TR03 classifiers within cycles defined by successive inhibitory bursts. Across the 18 equally weighted networks, the mean fractions were 76.35% with no spikes, 22.09% with one spike and 1.56% with two or more: 98.44% contained at most one spike (Fig. 6C–H). These distributions support predominantly one-spike participation among active neurons alongside widespread within-cycle silence, without establishing a strict firing ceiling or constant participating fraction.

The operating regime has asymmetric perturbation sensitivity

We next tested how the unpenalised TR02 loop-on and loop-off classifiers responded to disrupted spike transmission by deleting naturally generated spikes or inserting additional spike events during inference. Both networks maintained high accuracy when 80% of excitatory and inhibitory spike transmissions were randomly deleted: 88.9% for loop-off and 89.2% for loop-on. Deleting all transmissions reduced both to 10.6% (Fig. 7A). Adding extra spikes had a larger effect on loop-on networks. When spikes were added to both populations at each network’s baseline excitatory firing rate, accuracy fell to 76.4% for loop-off and 38.5% for loop-on. Doubling the rate of added spikes reduced accuracy further, to 53.1% and 11.4%, respectively (Fig. 7B). Loop-on therefore tolerated substantial deletion but was more sensitive to insertion at matched baseline-relative doses.

We next tested how inhibitory spike timing affected the unpenalised TR02 loop-on classifiers by artificially shifting and replaying recorded inhibitory spikes while preserving each neuron’s spike count within every presentation. Shifting each inhibitory spike separately reduced excitatory firing. Shifting all inhibitory spikes from the same fixed 22.8-ms time window together increased it (Fig. 8). At a proposed jitter standard deviation of 14 ms, independent shifts reduced mean excitatory firing from 16.6 to 0.0075 Hz and accuracy from 89.8% to 11.9%; group shifts increased firing to 68.3 Hz while accuracy declined to 82.5%. Replayed inhibitory firing remained 108.3 Hz. Thus, timing altered excitatory recruitment despite preserved inhibitory counts.

Accuracy across simulation timesteps

We next examined whether classification performance depended on the simulation timestep by comparing the TR04 loop-on classifiers trained separately at each tested timestep.

Mean test accuracy ranged from 88.37% to 89.63% across the twelvefold timestep range (Fig. 9). Mean excitatory firing was 14.27 Hz at the finest timestep and 17.02 Hz at the coarsest. Classification performance therefore persisted across the tested resolutions, while excitatory activity remained slightly timestep-dependent. Because each network was trained and evaluated at the same timestep, this comparison combines timestep effects during training and evaluation; it does not establish numerical convergence with fixed weights.

Loop-on networks classify continuously presented inputs

We next tested whether loop-on classifiers could classify successive images while carrying their hidden state from one presentation to the next. A TR06 loop-on classifier trained across variable input rates correctly classified five successive digits presented at changing durations and input rates (Fig. 10A). Sparse excitatory firing (Fig. 10B) and inhibitory volleys (Fig. 10C) persisted throughout the stream. The readout accumulated spike-count evidence during each presentation (Fig. 10D). We selected an example in which all five digits were classified correctly; the first candidate met this criterion. The example therefore demonstrates feasibility, not reliability across arbitrary streams. Output state and counts reset at supplied boundaries, so this result does not establish autonomous boundary detection. We then evaluated streaming accuracy across different presentation durations and input strengths, using repeated trials to assess performance beyond the illustrative example. At 25-Hz maximum-pixel input, increasing duration from 25 to 200 ms raised mean accuracy from 75.2% to 90.3% (Fig. 10E). At 200 ms, increasing input from 0.5 to 5 Hz raised accuracy from 26.3% to 82.7% (Fig. 10F). Mean accuracy then increased across the tested rates from 5 to 25 Hz, reaching 90.3%. Brief presentations and weak drive constrained performance. The grid paired image identity and order across conditions; its uncertainty remains conditional on one sampled image bank.

Methods

Experimental design

We combined untrained circuit experiments, a mean-field calculation and trained MNIST classifiers to examine how reciprocal excitatory–inhibitory coupling shapes population rhythms, excitatory firing and classification. In untrained circuits, we compared loop-off and loop-on responses across input strengths and mapped the effects of reciprocal coupling strength; the mean-field model examined oscillatory onset and its dependence on inhibitory decay time. Classifier experiments compared accuracy and firing under activity penalties, tested the effects of enabling reciprocal coupling during inference without retraining and of recurrent-weight training, and examined inhibitory timescales and cycle participation. We also measured sensitivity to spike deletion, spike insertion and inhibitory timing perturbations, assessed timestep dependence, and evaluated continuous image streams. Table 1 summarizes the designs, defines training groups TR01–TR06 and identifies classifier reuse.

Spiking networks

Hidden excitatory and inhibitory neurons followed conductance-based leaky integrate-and-fire dynamics:

Here, $V$ is membrane voltage and $t$ is time. Membrane capacitance $C_{m}$ was 1 nF for excitatory neurons and 0.5 nF for inhibitory neurons; leak conductance $g_{L}$ was 0.05 and 0.1 µS, respectively. The leak reversal potential $E_{L}$ was −65 mV, while excitatory and inhibitory reversal potentials $E_{E}$ and $E_{I}$ were 0 and −80 mV. Synaptic conductances $g_{E}$ and $g_{I}$ varied with incoming spikes, as described below. At the spike threshold $V_{\text{th}} = - 50$ mV, neurons emitted a spike and returned to the reset potential $V_{\text{reset}} = - 65$ mV. Voltage remained at reset for an absolute refractory period $\tau_{\text{ref}}$ of 1.2 ms in excitatory neurons and 0.6 ms in inhibitory neurons.

Synaptic conductances decayed exponentially and increased upon incoming spikes:

Here, $g_{x [ k ]}$ is conductance at step $k$, $w_{x j}$ the conductance increment from source $j$, and $s_{j [ k ]} \in \{ 0 , 1 \}$ its delivered spike indicator. Decay times were $\tau_{E} = \tau_{\text{AMPA}} = 2$ ms and ordinarily $\tau_{I} = \tau_{\text{GABA}} = 6$ ms. Exponential-Euler voltage integration used $\Delta t_{\text{sim}} = 0.1$ ms, except in the timestep study (Fig. 9). The inhibitory-timescale grid is specified below; exact updates and transmission ordering are in Appendix A1.

Classifier networks in TR01–TR06 contained input→E, E→I, I→E and E→output projections, with E→E and I→I connections disabled; loop-off additionally disabled reciprocal E–I coupling. Input and recurrent weights used lower-clamped Gaussian draws with fan-in normalization (Appendix A4 Table A1). Here, $N_{\text{pre}}$ is the presynaptic population size; input weights were initially zeroed with probability $q_{\text{zero}} = 0.95$ but remained trainable. Readout weights were dimensionless and initialized directly without fan-in normalization.

Untrained circuit experiments

We compared untrained loop-off and loop-on circuits containing 1,024 excitatory and 256 inhibitory neurons (Fig. 1). Reciprocal E→I and I→E coupling was absent in loop-off and enabled in loop-on, using the initialization parameters specified in Table A1. Each simulation lasted 400 ms at a 0.1-ms timestep, beginning with membrane voltages of −65 mV and zero synaptic conductances; no burn-in interval was discarded. Illustrative rasters used 1,024 independent Poisson input channels at 5 Hz per channel for loop-off and 45 Hz for loop-on. The matched input–output sweeps reused the same 1,024-input networks and fixed weights within each loop condition, testing both architectures at 2, 5, 10, 20, 40, 70 and 100 Hz per channel. Each architecture–input condition comprised one trial with seed 42. Population firing rates were measured over the full trial, and spectra were calculated from the illustrative raster trials as described below.

We mapped reciprocal coupling in untrained networks containing 1,024 excitatory and 256 inhibitory neurons (Fig. 2A–F). Each excitatory neuron received its own independent 100-Hz input channel through a fixed input weight of 0.5 µS. The E→I initialization parent mean $\mu_{E \to I}$ ranged from 0 to 3 µS in 0.3-µS increments, and the I→E parent mean $\mu_{I \to E}$ ranged from 0 to 6 µS in 0.6-µS increments, forming an $11 \times 11$ grid. Recurrent weights used Gaussian draws with parent standard deviations equal to 10% of their respective parent means, followed by lower clamping at zero and division by presynaptic population size. Each grid condition comprised one network and one trial with seed 42. Simulations lasted 1 s at a 0.1-ms timestep, with an inhibitory decay time of 6 ms. We discarded the first 100 ms and measured population firing rates and the excitatory-population rhythmicity metric (autocorrelation lobe–trough contrast) over the remaining 900 ms. Illustrative rasters show the first 200 ms after burn-in at three selected grid conditions.

Activity measurements

Mean per-neuron firing rates in Fig. 1F,H, Fig. 2A,B, Fig. 3D, Fig. 4C, Fig. 5B, Fig. 6A, Fig. 9 and Fig. 8C,D were calculated as

where $r_{P}$ is the firing rate of population $P$, $n_{\text{spike} , P}$ is its total spike count, $N_{P}$ is its neuron count and $T_{\text{obs}}$ is the observation duration in seconds, summed across presentations when pooling trials.

For spectral analysis, we formed each presentation’s excitatory-population trace by averaging binary spike indicators across neurons at each simulation step. We subtracted the trace mean and estimated its power spectral density using Welch’s method with a single Hann window spanning the full presentation. For the inhibitory-timescale classifiers in Fig. 2I and Fig. 6A, spectra were averaged across the 1,000 test presentations within each network before locating the largest peak between 5 and 150 Hz. The spectral-peak frequency $f_{\text{peak}}$ was refined by parabolic interpolation (Appendix B1). The 200-ms presentations gave a frequency-bin spacing of 5 Hz.

The illustrative spectra in Fig. 1E,G instead used individual 400-ms trials, giving 2.5-Hz bin spacing. Their peak annotation additionally required inhibitory spiking; absence of a marker was not a statistical test for rhythmicity. Interpolation formulas and boundary cases are specified in Appendix B.

We quantified excitatory-population temporal structure using the rhythmicity metric (Fig. 2C and Fig. 5C). Spikes were pooled across excitatory neurons into 1-ms count bins. We normalized and smoothed the uncentered autocorrelogram over lags up to 100 ms as specified in Appendix B2.

We identified the first local minimum at a lag of at least 2 ms and the maximum at positive lags preceding that minimum. The rhythmicity metric was defined as

Here, $R_{\text{contrast}}$ is the dimensionless rhythmicity metric, and $A_{\text{lobe}}$ and $A_{\text{trough}}$ are the preceding maximum and first minimum of the smoothed autocorrelation. The coupling-grid rhythmicity metric used the 900-ms post-burn-in interval; the classifier rhythmicity metric used the fixed reference-image presentation below. Appendix B specifies the feature-selection and boundary rules.

Mean-field calculation

We examined oscillatory onset using a separate four-variable mean-field model (Fig. 2G–I). Its state comprised excitatory and inhibitory population rates $r_{E}$ and $r_{I}$, excitatory conductance onto inhibitory neurons $g_{E}^{I}$, and inhibitory conductance onto excitatory neurons $g_{I}^{E}$:

Here, dots denote derivatives with respect to time in milliseconds; rates were expressed in ${\text{ms}}^{- 1}$ and conductances in µS. The external current $I_{\text{ext}}$, measured in nA, drove only the excitatory population. The functions $\Phi_{E}$ and $\Phi_{I}$ mapped mean input current to steady-state firing rate using noisy leaky integrate-and-fire gains.

The rate-relaxation constants $\tau_{r , E} = 20$ ms and $\tau_{r , I} = 5$ ms matched the respective membrane time constants. Synaptic decay constants were $\tau_{\text{AMPA}} = 2$ ms and, unless varied, $\tau_{\text{GABA}} = 6$ ms. Fixed driving-force magnitudes were $\Delta V_{\text{exc}} = 65$ mV and $\Delta V_{\text{inh}} = 15$ mV; summed conductance increments were $G_{E \to I} = 1$ µS and $G_{I \to E} = 2$ µS. The gains used the leak, reset and threshold parameters specified above, refractory periods $\tau_{\text{ref} , E} = 1.2$ ms and $\tau_{\text{ref} , I} = 0.6$ ms, and an effective voltage-noise scale $\sigma_{V} = 4$ mV.

This phenomenological closure fixed synaptic driving forces at rest, omitted conductance-dependent shunting and assumed a noise scale that was neither measured nor fitted to spiking activity. The population equations were deterministic; noise entered through the gain functions. The gain integral and closure assumptions are specified in Appendix C.

We continued fixed points across external currents $I_{\text{ext}} = 0$–4 nA and located onset at $I_{\text{ext}}^{\ast}$, where the leading complex-conjugate eigenvalue pair first crossed from negative to nonnegative real part (Fig. 2G). Eigenvalues came from the continuous-time Jacobian $J_{\text{flow}}$ of the four-variable model; continuation and refinement details are given in Appendix C2. The onset frequency was

Here, $\lambda_{J}^{\ast}$ is either eigenvalue of the refined crossing pair and $\operatorname{Im}$ denotes its imaginary part in radians per millisecond; $f_{\text{Hopf}}$ is in hertz. The factor 1000 converts cycles per millisecond to hertz. Numerical settings are specified in Appendix C.

We assessed criticality with 2-s upward/downward drive integrations, measuring peak-to-peak E-rate amplitude $A_{\text{pp}}$ in ${\text{ms}}^{- 1}$ over the final 500 ms at 1,001 fixed sample times (Fig. 2H). Numerical supercriticality required no resolved hysteresis and a positive amplitude-squared slope, using Appendix C’s predefined thresholds. No first Lyapunov coefficient was calculated, so narrower bistability or an unstable periodic orbit was not excluded. We repeated onset refinement at the six inhibitory decay times used for spiking classifiers, holding other mean-field parameters fixed. Fig. 2I compares the mean-field Hopf eigenfrequency $f_{\text{Hopf}}$ with the median finite-drive spectral-peak frequency $f_{\text{peak}}$ across three spiking classifiers per decay time.

Classifier training

We defined six training groups, TR01–TR06, by architecture, data partition, readout and varied training settings (Table 1). Each label identifies a group of training regimes, not an individual network or execution. Frequency and cycle-participation measurements reused TR03; loop-insertion and spike-perturbation experiments reused unpenalised TR02 classifiers. Unless specified otherwise, trained-network evaluations used the common post-training evaluation protocol below.

For trained-network comparisons, the unit of replication was one separately initialized and trained network, with three replicates per condition (seeds 42–44). These seeds controlled network initialization and, for TR06 variable-rate training, the sampling of input rates. The training-data split and post-training evaluation test-image subset were fixed separately. Repeated interventions on a network were repeated measurements; individual weights, neurons, presentations and neuron–cycle pairs were not additional network replicates.

Each network’s outcome was calculated before aggregation across networks. Sample SD used the denominator $n - 1$, and SEM was calculated as $\frac{\text{SD}}{\sqrt{n}}$, where $n = 3$ is the number of training replicates.

MNIST images were flattened into 784-element vectors and their pixel intensities divided by 255. For the six TR01 reference classifiers in Fig. 3A–C, we used the full official 60,000-image training partition, divided into 54,000 optimization images and 6,000 validation images using a class-stratified split. For the remaining 84 classifiers in TR02–TR06, we sampled 7,000 images uniformly without replacement from the same official partition and split them into 6,300 optimization and 700 validation images. Sampling and class-stratified splitting used seed 42, giving fixed data partitions across training replicates within each group. TR02 included both penalised and unpenalised classifiers. Loop-insertion and spike-perturbation experiments reused unpenalised TR02 classifiers without retraining. The official 10,000-image test partition remained separate from optimization and checkpoint selection. Post-training evaluations in Figs. 3–9 used a common subset of 1,000 test images sampled uniformly without replacement using seed 42, without class stratification. This subset and its presentation order were fixed independently of training seed. Single-image probes and illustrative rasters used separately specified examples; continuous-stream evaluations sampled from the full official test partition, as described below.

Each normalized pixel drove an independent Bernoulli spike channel:

Here, $s_{j [ k ]} \in \{ 0 , 1 \}$ indicates an input spike from pixel $j$ at timestep $k$, $a_{j} \in [ 0 , 1 ]$ is its normalized intensity, and $p_{j}$ is its per-step spike probability. The maximum-pixel rate $r_{\text{input,max}}$ was ordinarily 25 Hz, and the integration timestep $\Delta t_{\text{sim}}$ is expressed in seconds in this equation. Draws were independent across pixels, timesteps and presentations, with at most one spike per channel per step.

Independent image presentations lasted nominally 200 ms and began with freshly initialized hidden and output states. Training generated fresh encoding draws on each presentation. Spike encodings used for post-training evaluation were matched across networks at common input settings, with independent perturbation randomness (Appendix B5). Variable input rates and state continuity between successive images in the streaming experiment are described below.

All classifiers in TR01–TR06 had ten dimensionless, current-based spiking leaky integrate-and-fire output units, one per digit class. Exponential decay and input integration followed standard leaky integrate-and-fire theory.[2] The input normalization and immediate subtractive reset specified our implemented readout convention. Their E→output synapses delivered weighted excitatory spikes directly, without a separate synaptic decay state or a conductance-dependent driving force. The pre-reset output state followed

Here, $c \in \{ 0 , \ldots , 9 \}$ denotes the digit class, $k$ the simulation-step index, $u_{c}$ the dimensionless state retained after reset and $\tilde{u}_{c}$ its pre-reset value. The dimensionless weight $W_{j c}^{\text{out}}$ connects hidden excitatory neuron $j$ to output class $c$, and $s_{j}^{E} [ k + 1 ] \in \{ 0 , 1 \}$ denotes its spike delivered in that step. The decay factor is $\beta_{\text{out}} = \exp ( - \Delta t_{\text{sim}} / \tau_{\text{out}} )$, with output decay time $\tau_{\text{out}} = 2$ ms and integration timestep $\Delta t_{\text{sim}}$; $h = \Delta t_{\text{sim}} / ( 1 \, \text{ms} )$ specifies the dimensionless timestep scaling of the input. A pre-reset state at or above the dimensionless threshold $\theta_{\text{out}} = 1$ produced one spike, after which the threshold was subtracted from the state:

Here, $s_{c}^{\text{out}} [ k + 1 ] \in \{ 0 , 1 \}$ is the output spike indicator, and $\mathbf{1} [ \cdot ]$ equals one when its condition holds and zero otherwise. Output units had no refractory period and began each independent presentation at zero.

For TR01–TR05 classification (Figs. 3–9), we scored classes using the mean pre-reset state of these spiking units, rather than their emitted spike counts. Class logits and predictions were

Here, $\tilde{u}_{c [ k ]}$ is pre-reset output state at step $k$, $N_{t}$ the presentation’s step count, $z_{c}$ the class logit and $\hat{y}$ the predicted digit. The mean-state readout retained the output units’ spike generation and subtractive reset. Exact output updates are in Appendix A2; TR06 classifiers instead used the spike-count readout described below.

All 90 included classifiers completed 50 training epochs using AdamW with learning rate $\eta = 0.0004$, zero weight decay and minibatches of 256 images. The objective was mean cross-entropy, supplemented by the activity penalty described below where applicable. We used full, untruncated surrogate-gradient backpropagation through the complete temporal unrolling of every presentation in each minibatch.[3] Before each optimizer update, we clipped the global Euclidean gradient norm to a maximum of 1. Trainable connection weights were subsequently clamped below at zero.

Spikes remained binary during forward simulation. During backpropagation, the threshold derivative was replaced by the fast-sigmoid surrogate

where $s$ is the spike indicator, $\xi$ is the numerical distance from threshold—expressed in millivolts for hidden neurons and normalized units for output neurons—and the surrogate slope $\beta = 1$. Gradients through each hidden excitatory and inhibitory membrane increment were additionally divided by the dimensionless voltage-gradient damping divisor $d_{\text{grad}} = 1 , 000$ in both loop-on and loop-off. Damping altered backward derivatives, not the forward update rule; output updates were undamped. Exact gradient paths are specified in Appendix A3.

After each training epoch, we evaluated the corresponding validation partition—6,000 images for TR01 and 700 for TR02–TR06—using three independently seeded spike-encoding draws. The same encoding seeds were reused across epochs, keeping stochastic inputs fixed for checkpoint comparisons. For TR06 classifiers, validation input-rate draws were also fixed across epochs. Cross-entropy and accuracy were averaged over validation images and encoding draws; the activity penalty was excluded from the validation loss.

Best-validation checkpoints minimized mean validation cross-entropy over epochs 1–50, with ties resolved by higher accuracy, then earlier epoch. Final checkpoints were from epoch 50. Table 1 specifies the checkpoint used for each analysis.

Classifier comparisons and recurrent coupling

For the accuracy–rate comparison (Fig. 3D), we added a one-sided quadratic penalty to the classification loss:

Here, $L_{\text{total}}$ is the objective, $L_{\text{CE}}$ mean classification cross-entropy, $B$ the minibatch size and $b$ the presentation index. The excitatory rate $r_{E , b}$ is that presentation’s excitatory spike count divided by excitatory population size and duration in seconds. Division by $1 \, \text{Hz}$ makes the excess dimensionless. We used dimensionless penalty coefficient $\lambda_{\text{rate}} = 0.041$ at ceilings $r_{E , \text{ceil}}$ of 1, 2.5, 5, 10 and 25 Hz; the unpenalised condition used $\lambda_{\text{rate}} = 0$. The penalty acts on each presentation before minibatch averaging and imposes no hard rate limit. Only input/readout weights were trained; recurrence stayed fixed, enabled in loop-on and disabled in loop-off.

We compared loop-off and loop-on under the six activity conditions defined above, giving 36 reduced-pool networks in TR02 (Fig. 3D). Final-checkpoint test accuracy and excitatory firing rates used the common post-training evaluation protocol and were summarized as means and SEM across three training replicates per condition. Separately, the six TR01 full-MNIST reference classifiers supplied the illustrative rasters and validation learning curves (Fig. 3A–C). Learning curves averaged three replicate accuracies at each epoch without uncertainty intervals; illustrative-example selection is specified in Appendix B4.

The accuracy–rate, loop-insertion and spike-perturbation evaluations used the 18 TR02 loop-off classifiers trained with $d_{\text{grad}} = 1 , 000$, together with the corresponding loop-on classifiers where applicable (Fig. 3, Fig. 4 and Fig. 7). The remaining classifier measurements used the networks and evaluation protocols specified in Table 1.

We inserted reciprocal inhibition into the three unpenalised TR02 loop-off classifiers, holding learned input and readout weights fixed without retraining (Fig. 4).

The dimensionless loop strength $s$ ranged from 0 to 1 in increments of 0.1. At each strength, we newly initialized the E→I and I→E matrices from Gaussian distributions with parent means $\mu_{E \to I} = s$ µS and $\mu_{I \to E} = 2 s$ µS, respectively, and parent standard deviations equal to 10% of each parent mean. We clamped negative draws to zero and divided weights by the presynaptic population size: 1,024 for E→I and 256 for I→E. Thus, $s = 0$ disabled the loop and $s = 1$ used the standard recurrent initialization.

At each strength, we measured post-training test accuracy and E/I firing rates, summarized as means and sample SD across the three reused classifiers.

We compared fixed recurrence at standard strength $s = 1$ with trainable recurrence initialized at $s = 1$, $0.1$ or $0$ (Fig. 5), giving 12 networks in TR05. All conditions trained input and readout weights without an activity penalty. Post-training test accuracy and E/I firing rates used the common protocol. These outcomes and the rhythmicity metric were summarized as means and SEM across three training replicates per condition. The rhythmicity metric used the final training diagnostic from one 200-ms encoding of the first official-test digit-0 image, generated with seed 0 and fixed across networks, without smoothing across epochs.

For TR05 recurrent-weight comparisons, initial matrices were reconstructed using the original initialization seeds and parameters. At initialization and epoch 50, we pooled the three complete matrices separately for each connection direction and condition, giving 786,432 entries per pooled distribution. We calculated the fraction of strictly positive weights and the arithmetic mean weight across all entries, including zeros. Because the three matrices had equal size, each network contributed equally to these pooled summaries. We also calculated the same two statistics separately for each network and calculated their SEM across the three networks, separately at initialization and epoch 50. These were descriptive summaries, without significance tests of training-induced weight changes. Jointly trained weights within a matrix were not independent training replicates.

Inhibitory timescale and cycle participation

TR03 comprised 18 separately trained loop-on classifiers at inhibitory decay times $\tau_{\text{GABA}} = 4.5 , 6 , 9 , 12 , 18$ and $27$ ms (Fig. 6A–B). Input and readout weights were trained without an activity penalty; recurrence remained fixed at standard strength. Evaluation retained each assigned decay time and measured post-training test accuracy, excitatory firing rate and spectral-peak frequency using the common protocols. Fig. 6A–B summarizes these outcomes as means and SEM across the three training replicates at each decay time.

We fitted the affine relationship

by ordinary least squares, giving equal weight to each of the six condition means. Here, $\bar{r}_{E}$ and $\bar{f}_{\text{peak}}$ are the across-replicate mean excitatory firing rate and spectral-peak frequency, both in hertz; $a$ is the fitted intercept in hertz and $\beta_{r f}$ is the dimensionless fitted slope. Fit quality was quantified by the coefficient of determination $R_{\text{fit}}^{2}$, using the centered total sum of squares. The mean-field comparison in Fig. 2I used medians of the same three per-network frequencies at each decay time.

We reused the TR03 classifiers’ recorded spikes to measure excitatory counts per cycle (Fig. 6C–H), without new inference or training. Inhibitory bursts were detected from smoothed population counts with peak separation scaled to each network’s spectral period (Appendix B3).

Cycle boundaries were midpoints between inhibitory-burst peaks, with edge intervals extending to presentation boundaries. Presentations without bursts were excluded. Each excitatory neuron–cycle pair was classified as containing 0, 1, 2 or ≥3 spikes. We normalized each network’s four counts separately across its presentations and averaged its fractions equally with the other two networks in its condition, or across all 18 networks for the overall summary. This gave each network equal weight irrespective of its detected cycle total and balanced the six decay conditions in the overall summary. These equal-network distributions were descriptive means without uncertainty intervals. Neither the distributions nor the rate–frequency fit treated neuron–cycle pairs as independent replicates. Appendix B specifies discretization and single-burst conventions.

Cycle-participation analysis excluded twelve of 18,000 presentations without detected inhibitory bursts, all from the TR03 seed-43 network at 27-ms inhibitory decay. Weak image inputs, comparatively weak input connections and low firing rates left excitatory activity too dispersed to recruit inhibition. The same images recruited inhibition in other replicates, suggesting an interaction between stimulus strength, network properties and stochastic encoding.

Spike perturbations

We separately applied deletion and insertion to the six unpenalised TR02 loop-off and loop-on classifiers without retraining (Fig. 7A–B), using the common post-training evaluation protocol and matched input encodings across perturbation conditions. Outcomes were summarized as means and sample SD across three training replicates per architecture.

Deletion independently suppressed each naturally generated excitatory or inhibitory spike with probability $p_{\text{del}} = 0 , 0.1 , \ldots , 1$, preserving the neuronal reset and refractory period associated with that spike.

For insertion, we first measured each network’s unperturbed mean per-neuron excitatory firing rate, $r_{E , 0}$, across the complete test subset. The nominal per-neuron insertion rate was $r_{\text{add}} = \alpha r_{E , 0}$, where the dimensionless dose $\alpha = 0 , 0.1 , \ldots , 2$. Independent Bernoulli events were sampled for every excitatory and inhibitory neuron at each timestep, with probability equal to the insertion rate in hertz multiplied by the timestep in seconds. The same nominal insertion rate applied to both populations. Transmission was capped at one spike per neuron per timestep, so insertion coinciding with an existing spike added no event. Inserted events triggered neither resets nor refractory periods.

Modified events supplied the readout and subsequent recurrent transmission. Insertion conditions were aggregated at matched baseline-relative doses.

We applied two count-preserving inhibitory replay perturbations to the three unpenalised TR02 loop-on classifiers (Fig. 8). Unperturbed inhibitory spikes recorded under the common post-training evaluation protocol were shifted and replayed, replacing naturally generated inhibitory outputs while excitatory activity and readout responses were recomputed with unchanged weights and identical inputs.

For the conditions reported in Fig. 8, independent-spike jitter assigned each inhibitory event a zero-mean Gaussian offset with standard deviation $\sigma = 0 , 0.5 , 1 , 2 , 5 , 9$ or $14$ ms. Group jitter assigned one shared Gaussian offset to all inhibitory events originating within each fixed 22.8-ms window of a presentation, using $\sigma = 0 , 1 , 3 , 7$ or $14$ ms. These windows comprised 228 simulation steps and were defined independently of detected bursts.

Each inhibitory neuron’s spike count was preserved within every presentation. The parameter $\sigma$ describes proposed offsets before boundary and collision adjustments; both arms shared the zero-jitter replay control. Exact transformations are specified in Appendix D.

We measured post-training test accuracy and excitatory and replayed inhibitory rates, reporting means across three reused classifiers without uncertainty intervals. Delivered inhibition could not respond to ongoing excitatory activity.

Simulation timestep comparison

We compared the 15 separately trained TR04 loop-on classifiers (Fig. 9): three training replicates (seeds 42–44) at each integration timestep of 0.05, 0.1, 0.2, 0.3 and 0.6 ms. All networks underwent 50 training epochs without an activity penalty, with recurrent weights fixed at standard strength. Excitatory and inhibitory refractory periods remained 1.2 and 0.6 ms, respectively, corresponding to integer step counts throughout.

Each epoch-50 checkpoint was evaluated at its training timestep on the common 1,000-image test subset. Nominally 200-ms presentations were truncated to complete simulation steps, giving realised durations of 199.8 ms at 0.3 and 0.6 ms and 200 ms otherwise. Mean per-neuron excitatory firing rates included all excitatory neurons and test presentations, using the realised duration. Accuracy and firing rate were summarised as means and SEM across the three training replicates at each timestep.

Continuous streams

To identify input rates that retained enough digit information for classification, we encoded MNIST pixels as stochastic spikes and filtered each pixel’s events through synaptic conductance and a non-spiking membrane. Each pixel’s mean voltage above rest formed a feature for nonlinear feed-forward decoders. Their held-out accuracy across input rates informed the rates selected for variable-rate training of the three TR06 classifiers later evaluated on continuous streams.

These classifiers used the common isolated-image training protocol without an activity penalty, training input and readout weights while keeping recurrence fixed at standard strength (Fig. 10). For each image presentation, the maximum-pixel input rate, $r_{\text{input,max}}$, was drawn with equal probability from the decoder-informed set: 0.5, 0.75, 1, 1.5, 2, 3, 5, 7.5, 10, 15 and 25 Hz. The rate remained constant throughout the presentation and was sampled independently for each image within a minibatch, then resampled on subsequent presentations.

TR06 classifiers used spiking leaky integrate-and-fire output units with the same decay, threshold and subtractive-reset dynamics described above. Their readout differed from ordinary classification: class logits were accumulated output-spike counts rather than mean pre-reset states:

Here, $s_{c}^{\text{out} [ k ]} \in \{ 0 , 1 \}$ indicates an output spike; class index $c$, logit $z_{c}$ and presentation length $N_{t}$ are defined above. Counts entered cross-entropy directly without decay or duration normalization, using the same surrogate derivative. Classification selected the largest final count, breaking ties by lowest class index.

Readout weights were initialized from a Gaussian distribution with dimensionless parent mean $\mu_{\text{out}} = 0.05$ and parent standard deviation $\sigma_{\text{out}} = 0.04$ (Appendix A4), clamped below at zero without fan-in normalization, and then trained alongside the input weights. Both sets of weights were held fixed during evaluation. Evaluations used each network’s best-validation checkpoint. The same three validation-selected TR06 classifiers were evaluated across all duration–rate conditions. During continuous streams, output state and accumulated counts were reset to zero at each supplied image boundary, while hidden neuronal and synaptic states continued between images.

For quantitative streaming evaluation (Fig. 10E–F), we crossed presentation durations of 25, 50, 100 and 200 ms with the eleven maximum-pixel input rates used during training. Each of the three frozen TR06 classifiers processed 40 five-image streams per duration–rate combination, in batches of five streams. This gave 200 classification decisions per network and condition, and 26,400 decisions overall. Duration and input rate remained constant within each stream, and images followed without gaps. Each stream began with freshly initialized hidden and output states; subsequent image boundaries used the state-handling procedure described above.

For each stream, five images were sampled uniformly without replacement from the full 10,000-image official MNIST test partition, without class stratification. Digit labels could repeat within a stream, and images could recur across streams. We prespecified one bank of 40 ordered streams and reused the same image indices for every network, duration and rate. Thus, conditions shared both the target images and the preceding-image order. Spike encoding used separate generators with distinct seeds for all 5,280 network–condition–stream combinations (Appendix B5). The 26,400 decisions are repeated evaluations of this image bank, not that many independently sampled images.

Accuracy included all 200 decisions per network and condition, including presentations with no output spikes. Fig. 10E shows means across the three training replicates without uncertainty intervals; Fig. 10F reuses the 200-ms evaluations and shows means and SEM across those same replicates. Presentation and count-accumulation windows had equal duration, so their effects were not separated. We did not compare continuing hidden state with a hidden-state-reset control or test autonomous image segmentation.

For continuous-stream evaluation, all networks shared one image bank but used different encoding draws. Its SEM therefore summarizes variation across these trained networks and their encodings, conditional on the sampled images. One bank and one encoding draw per network, condition and stream do not separately estimate image-sampling, encoding and training variability. Image pairing does not make the 200 sequential decisions per condition independent network replicates.

We reused the TR06 illustrative recording (Fig. 10A–D), selected independently of quantitative evaluation using a predefined five-correct criterion, met by the first candidate (Appendix B4). It was not a stream from the shared quantitative image bank.

To visualize the evolving readout within each presentation (Fig. 10D), we transformed cumulative output-spike counts into softmax shares:

Here, $n_{c [ k ]}$ counts class-$c$ output spikes through step $k$ since the supplied boundary; $q_{c [ k ]}$ is its dimensionless share, with $d$ indexing classes. These uncalibrated shares visualize evidence; the 0.5 line is not a decision threshold.

Software and hardware

The study software was developed iteratively with AI coding agents (OpenAI Codex and Claude Code) under the lead author’s direction. The lead author specified the scientific models and experimental protocols, directed implementation and revisions, and retained responsibility for software acceptance and its scientific use. Verification included automated unit, integration and numerical-regression testing. Simulator tests cover neuronal and synaptic dynamics, spike and refractory timing, surrogate gradients, checkpoint restoration and training continuation. Numerical checks include analytical solutions, preserved reference outputs, and comparisons with independent Brian2 implementations of conductance-based dynamics and snnTorch implementations of current-based dynamics under matched conditions. Experiment-level tests cover checkpoint selection, spike perturbations, analysis outputs and data provenance. The implementations and test suites are publicly available for inspection and reuse.

The experiment code is maintained in the Pinglab experiments repository, and the spiking neural network (SNN) simulator is maintained in the SNNLab simulator repository. Spiking simulation, training and evaluation used custom Python/PyTorch code; all 90 classifiers were trained with PyTorch 2.11.0+cu128 on NVIDIA A100-SXM4 GPUs with 80 GB of device memory, using one CUDA device per classifier. Mean-field calculations used Python 3.10.19, NumPy 2.2.6 and SciPy 1.15.3. For TR01, mean training-loop and validation times over 50 epochs were 402 and 158 min for loop-off and 429 and 161 min for loop-on, respectively, averaged across three seeds per architecture. Mean epoch times, including training and validation, were 11.2 min for loop-off and 11.8 min for loop-on over epochs 2–50, excluding first-epoch warm-up, reference-image diagnostics and subsequent logging or checkpoint-file writes.