Manuscript

Meta

1. Results: 1,526 words

2. Methods: 4,465 words

3. Captions: 1,638 words

4. Appendix A: 1,506 words

5. Appendix B: 1,235 words

6. Appendix C: 764 words

7. Appendix D: 379 words

TODO

1. Review P9 and Figure 5 with Tim.

Results

Reciprocal coupling creates a gamma-rhythmic low-rate regime

P1

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

P2

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, with 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).

P3

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 autocorrelation lobe–trough contrast (Appendix B2) near zero (Fig. 2A–C). With both pathways present, stronger coupling reduced excitatory firing to single-digit rates, sustained inhibitory firing and increased contrast towards one across a broad region.

P4

Selected rasters illustrate the same progression (Fig. 2D–F): autocorrelation contrast 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.

P5

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 ${𝐼}_{\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 ${𝑅}^2=0.999$ met the predefined numerical criteria for supercritical onset; criticality was not established analytically.

P6

We then compared how slowing inhibition affected oscillation frequency in the mean-field model and trained spiking classifiers. 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

P7

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. Both architectures reached approximately 90% validation accuracy, with dense loop-off firing and recurring loop-on volleys (Fig. 3A–C). 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). Tightening the activity-penalty ceiling progressively reduced excitatory firing and classification accuracy in both architectures, but their responses differed. 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%).

P8

We next asked whether adding reciprocal inhibition after training could produce sparse rhythmic activity while preserving classification accuracy. Adding reciprocal inhibition to trained 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.

P9

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 autocorrelation contrast relative to fixed loop-on recurrence. Final autocorrelation lobe–trough contrast 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. Weak contrast accompanied asymmetric reorganization in these three-replicate comparisons, without identifying a causal weight change.

Excitatory firing is organised by population-cycle participation

P10

We next varied the inhibitory decay time in separately trained loop-on networks 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 ${𝑅}^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]

P11

To examine how individual neurons participated in these population rhythms, we counted each excitatory neuron’s spikes 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

P12

We next tested how 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.

P13

We next examined whether classification performance depended on the simulation timestep by separately training and evaluating loop-on networks at five temporal resolutions. Classification accuracy remained near 90% across tested timesteps, although excitatory firing varied. This supports robustness after separate training at each timestep, without establishing fixed-weight convergence (Appendix A5).

P14

We next tested how inhibitory spike timing affected 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.

Loop-on networks classify continuously presented inputs

P15

We next tested whether loop-on classifiers could classify successive images while carrying their hidden state from one presentation to the next. A loop-on classifier trained across variable input rates correctly classified five successive digits presented at changing durations and input rates (Fig. 9A). Sparse excitatory firing (Fig. 9B) and inhibitory volleys (Fig. 9C) persisted throughout the stream. The readout accumulated spike-count evidence during each presentation (Fig. 9D). 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. 9E). At 200 ms, increasing input from 0.5 to 5 Hz raised accuracy from 26.3% to 82.7% (Fig. 9F). 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

P17

We combined untrained circuit simulations, trained MNIST classifiers, interventions on those classifiers and a separate mean-field model (Table 1). Inhibitory-timescale classifiers supplied frequency and cycle-participation measurements; unpenalised accuracy–rate classifiers supplied loop-insertion and spike-perturbation experiments. Unless specified otherwise, trained-network analyses used the common endpoint protocol below; Table 1 and Statistical reporting define checkpoints and replication.

Table 1: Experimental design and checkpoint policy. E/I: excitatory/inhibitory neurons. Each training condition comprised three replicates. Final: epoch 50. Best validation: lowest mean validation cross-entropy across three encoding draws, with ties resolved by higher accuracy, then earlier epoch. ¹ Mean-field results. ² Spiking-model results.

Spiking networks

P18

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

Here, $𝑉$ is membrane voltage and $𝑡$ is time. Membrane capacitance ${𝐶}_{𝑚}$ was 1 nF for excitatory neurons and 0.5 nF for inhibitory neurons; leak conductance ${𝑔}_{𝐿}$ was 0.05 and 0.1 µS, respectively. The leak reversal potential ${𝐸}_{𝐿}$ was −65 mV, while excitatory and inhibitory reversal potentials ${𝐸}_{𝐸}$ and ${𝐸}_{𝐼}$ were 0 and −80 mV. Synaptic conductances ${𝑔}_{𝐸}$ and ${𝑔}_{𝐼}$ varied with incoming spikes, as described below. At the spike threshold ${𝑉}_{\text{th}}=-50$ mV, neurons emitted a spike and returned to the reset potential ${𝑉}_{\text{reset}}=-65$ mV. Voltage remained at reset for an absolute refractory period ${𝜏}_{\text{ref}}$ of 1.2 ms in excitatory neurons and 0.6 ms in inhibitory neurons.

P19

Synaptic conductances decayed exponentially and increased upon incoming spikes:

Here, ${𝑔}_{𝑥\left[𝑘\right]}$ is conductance at step $𝑘$, ${𝑤}_{𝑥𝑗}$ the conductance increment from source $𝑗$, and ${𝑠}_{𝑗\left[𝑘\right]}\in \left\{0,1\right\}$ its delivered spike indicator. Decay times were ${𝜏}_{𝐸}={𝜏}_{\text{AMPA}}=2$ ms and ordinarily ${𝜏}_{𝐼}={𝜏}_{\text{GABA}}=6$ ms. Exponential-Euler voltage integration used $\mathrm{\Delta }{𝑡}_{\text{sim}}=0.1$ ms, except in the timestep study (Appendix A5). The inhibitory-timescale grid is specified below; exact updates and transmission ordering are in Appendix A1.

Spiking measurements used E/I refractory holds of 1.2/0.6 ms. The timestep comparison, separate mean-field calculation, quantitative continuous-stream evaluation and loop-off training and comparisons used the same refractory durations. The quantitative continuous-stream evaluation used a shared image bank.

P20

Classifier networks 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, ${𝑁}_{\text{pre}}$ is the presynaptic population size; input weights were initially zeroed with probability ${𝑞}_{\text{zero}}=0.95$ but remained trainable. Readout weights were dimensionless and initialized directly without fan-in normalization.

Classifier training

P21

MNIST images were flattened into 784-element vectors and their pixel intensities divided by 255. We sampled 7,000 images uniformly without replacement from the official training partition, then divided this subset into 6,300 optimization images and 700 validation images using a class-stratified split. Both sampling and splitting used seed 42, giving the same data partitions across training replicates. The official 10,000-image test partition remained separate from optimization and checkpoint selection. Endpoint evaluations in Figs. 3–8 and Appendix Fig. A1 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.

P22

Each normalized pixel drove an independent Bernoulli spike channel:

Here, ${𝑠}_{𝑗\left[𝑘\right]}\in \left\{0,1\right\}$ indicates an input spike from pixel $𝑗$ at timestep $𝑘$, ${𝑎}_{𝑗}\in \left[0,1\right]$ is its normalized intensity, and ${𝑝}_{𝑗}$ is its per-step spike probability. The maximum-pixel rate ${𝑟}_{\text{input,max}}$ was ordinarily 25 Hz, and the integration timestep $\mathrm{\Delta }{𝑡}_{\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. Endpoint encodings 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.

P23

The output layer contained ten leaky integrate-and-fire units, one per digit class. Each unit had a dimensionless state ${𝑢}_{𝑐}$, where $𝑐\in \left\{0,\dots ,9\right\}$ denotes the class, and an output decay time ${𝜏}_{\text{out}}=2$ ms. Weighted excitatory spikes drove this state. At each timestep, a pre-reset state at or above the output threshold ${𝜃}_{\text{out}}=1$ produced one spike, after which the threshold was subtracted from the state. Output units had no refractory period and began each independent presentation at zero.

For ordinary classification (Figs. 3–8 and Appendix Fig. A1), class logits and predictions were

Here, ${\widetilde{𝑢}}_{𝑐\left[𝑘\right]}$ is pre-reset output state at step $𝑘$, ${𝑁}_{𝑡}$ the presentation’s step count, ${𝑧}_{𝑐}$ the class logit and $\widehat{𝑦}$ the predicted digit. Exact output updates are in Appendix A streaming uses the spike-count readout below.

P24

We trained classifiers for 50 epochs using AdamW with learning rate $𝜂=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.[2] 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 $𝑠$ is the spike indicator, $𝜉$ is the numerical distance from threshold—expressed in millivolts for hidden neurons and normalized units for output neurons—and the surrogate slope $𝛽=1$. Gradients through each hidden excitatory and inhibitory membrane increment were additionally divided by the dimensionless voltage-gradient damping divisor ${𝑑}_{\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.

P25

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

Here, ${𝐿}_{\text{total}}$ is the objective, ${𝐿}_{\text{CE}}$ mean classification cross-entropy, $𝐵$ the minibatch size and $𝑏$ the presentation index. The excitatory rate ${𝑟}_{𝐸,𝑏}$ 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 ${𝜆}_{\text{rate}}=0.041$ at ceilings ${𝑟}_{𝐸,\text{ceil}}$ of 1, 2.5, 5, 10 and 25 Hz; the unpenalised condition used ${𝜆}_{\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.

P26

After each training epoch, we evaluated the 700 validation images using three independently seeded spike-encoding draws. The same encoding seeds were reused across epochs, keeping stochastic inputs fixed for checkpoint comparisons. For variable-rate 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.

Untrained circuit experiments

P27

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.

P28

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 ${𝜇}_{𝐸\to 𝐼}$ ranged from 0 to 3 µS in 0.3-µS increments, and the I→E parent mean ${𝜇}_{𝐼\to 𝐸}$ ranged from 0 to 6 µS in 0.6-µS increments, forming an $11\times 11$ grid. Recurrent weights used Gaussian draws with standard deviations equal to 10% of their respective 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 excitatory-population autocorrelation contrast over the remaining 900 ms. Illustrative rasters show the first 200 ms after burn-in at three selected grid conditions.

Activity measurements

P29

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

where ${𝑟}_{𝑃}$ is the firing rate of population $𝑃$, ${𝑛}_{\text{spike},𝑃}$ is its total spike count, ${𝑁}_{𝑃}$ is its neuron count and ${𝑇}_{\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 ${𝑓}_{\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.

P30

We quantified excitatory-population temporal structure using autocorrelation lobe–trough contrast (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. Contrast was defined as

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

Mean-field calculation

P31

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

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

The rate-relaxation constants ${𝜏}_{𝑟,𝐸}=20$ ms and ${𝜏}_{𝑟,𝐼}=5$ ms matched the respective membrane time constants. Synaptic decay constants were ${𝜏}_{\text{AMPA}}=2$ ms and, unless varied, ${𝜏}_{\text{GABA}}=6$ ms. Fixed driving-force magnitudes were $\mathrm{\Delta }{𝑉}_{\text{exc}}=65$ mV and $\mathrm{\Delta }{𝑉}_{\text{inh}}=15$ mV; summed conductance increments were ${𝐺}_{𝐸\to 𝐼}=1$ µS and ${𝐺}_{𝐼\to 𝐸}=2$ µS. The gains used the leak, reset and threshold parameters specified above, refractory periods ${𝜏}_{\text{ref},𝐸}=1.2$ ms and ${𝜏}_{\text{ref},𝐼}=0.6$ ms, and an effective voltage-noise scale ${𝜎}_{𝑉}=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.

P32

We continued fixed points across external currents ${𝐼}_{\text{ext}}=0$–4 nA and located onset at ${𝐼}_{\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 ${𝐽}_{\text{flow}}$ of the four-variable model; continuation and refinement details are given in Appendix C2. The onset frequency was

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

P33

We assessed criticality with 2-s upward/downward drive integrations, measuring peak-to-peak E-rate amplitude ${𝐴}_{\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 ${𝑓}_{\text{Hopf}}$ with the median finite-drive spectral-peak frequency ${𝑓}_{\text{peak}}$ across three spiking classifiers per decay time.

Classifier comparisons and recurrent coupling

P34

We compared loop-off and loop-on under the six activity conditions defined above, giving 36 networks (Fig. 3). Final-checkpoint test accuracy and excitatory firing rates used the common endpoint protocol. Unpenalised validation learning curves averaged replicate accuracies at each epoch; illustrative-example selection is specified in Appendix B4.

The accuracy–rate, loop-insertion and spike-perturbation evaluations used the 18 loop-off classifiers trained with ${𝑑}_{\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.

P35

We inserted reciprocal inhibition into the unpenalised loop-off classifiers from the accuracy–rate comparison, holding learned input and readout weights fixed without retraining (Fig. 4).

The dimensionless loop strength $𝑠$ 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 ${𝜇}_{𝐸\to 𝐼}=𝑠$ µS and ${𝜇}_{𝐼\to 𝐸}=2𝑠$ µS, respectively, and standard deviations equal to 10% of each 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, $𝑠=0$ disabled the loop and $𝑠=1$ used the standard recurrent initialization.

At each strength, we measured endpoint accuracy and E/I firing rates.

P36

We compared fixed recurrence at standard strength $𝑠=1$ with trainable recurrence initialized at $𝑠=1$, $0.1$ or $0$ (Fig. 5), giving 12 networks. All conditions trained input and readout weights without an activity penalty. Endpoint accuracy and E/I firing rates used the common protocol. Contrast 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 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

P37

We separately trained 18 loop-on classifiers at inhibitory decay times ${𝜏}_{\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 endpoint accuracy, excitatory firing rate and spectral-peak frequency using the common protocols.

We fitted the affine relationship

by ordinary least squares, giving equal weight to each of the six condition means. Here, ${\overline{𝑟}}_{𝐸}$ and ${\overline{𝑓}}_{\text{peak}}$ are the across-replicate mean excitatory firing rate and spectral-peak frequency, both in hertz; $𝑎$ is the fitted intercept in hertz and ${𝛽}_{𝑟𝑓}$ is the dimensionless fitted slope. Fit quality was quantified by the coefficient of determination ${𝑅}_{\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.

P38

We reused the same 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. Appendix B specifies discretization and single-burst conventions.

Spike perturbations and numerical resolution

P39

We separately applied deletion and insertion to unpenalised loop-off and loop-on classifiers without retraining (Fig. 7A–B), using the common endpoint protocol and matched input encodings across perturbation conditions.

Deletion independently suppressed each naturally generated excitatory or inhibitory spike with probability ${𝑝}_{\text{del}}=0,0.1,\dots ,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, ${𝑟}_{𝐸,0}$, across the complete test subset. The nominal per-neuron insertion rate was ${𝑟}_{\text{add}}=𝛼{𝑟}_{𝐸,0}$, where the dimensionless dose $𝛼=0,0.1,\dots ,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.

P40

We evaluated separately trained loop-on classifiers at matched training and evaluation timesteps while holding refractory durations fixed. The full protocol and results are given in Appendix A5.

P41

We applied two count-preserving inhibitory replay perturbations to unpenalised loop-on classifiers (Fig. 8). Unperturbed inhibitory spikes recorded under the common endpoint 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 $𝜎=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 $𝜎=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 $𝜎$ 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 endpoint accuracy and excitatory and replayed inhibitory rates. Delivered inhibition could not respond to ongoing excitatory activity.

Continuous streams

P42

Continuous-stream 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. 9).

For each presentation, the maximum-pixel input rate, ${𝑟}_{\text{input,max}}$, was sampled uniformly from 0.5, 0.75, 1, 1.5, 2, 3, 5, 7.5, 10, 15 and 25 Hz and held constant throughout that presentation. Rates were sampled independently for individual images within each minibatch and resampled on subsequent presentations.

Continuous-stream classifiers used the output units defined above, with class logits given by accumulated spike counts:

Here, ${𝑠}_{𝑐}^{\text{out}\left[𝑘\right]}\in \left\{0,1\right\}$ indicates an output spike; class index $𝑐$, logit ${𝑧}_{𝑐}$ and presentation length ${𝑁}_{𝑡}$ 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 drawn from a Gaussian distribution with dimensionless mean ${𝜇}_{\text{out}}=0.05$ and standard deviation ${𝜎}_{\text{out}}=0.04$, then clamped below at zero, without fan-in normalization. Evaluations used each network’s best-validation checkpoint with weights held fixed. The same three validation-selected 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.

P43

For quantitative streaming evaluation (Fig. 9E–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 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. 9F reuses the 200-ms evaluations from Fig. 9E aggregation and uncertainty follow the reporting protocol below. 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.

P44

We reused the illustrative recording 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, we transformed cumulative output-spike counts into softmax shares:

Here, ${𝑛}_{𝑐\left[𝑘\right]}$ counts class-$𝑐$ output spikes through step $𝑘$ since the supplied boundary; ${𝑞}_{𝑐\left[𝑘\right]}$ is its dimensionless share, with $𝑑$ indexing classes. These uncalibrated shares visualize evidence; the 0.5 line is not a decision threshold.

Statistical reporting and reproducibility

P45

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 variable-rate training, the sampling of input rates. The training-data split and endpoint 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. Untrained circuit probes used a single seed, and mean-field calculations were deterministic.

Each network’s outcome was calculated before aggregation across networks. Displayed error bars or bands represent SEM in Figs. 3, 5A–G, 6A–B and 9F and Appendix Fig. A1, and sample SD in Figs. 4 and 7A–B. Sample SD used the denominator $𝑛-1$, and SEM was calculated as $\frac{\text{SD}}{\sqrt{𝑛}}$, where $𝑛=3$ is the number of training replicates.

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.

Figures 8 and 9E show means without uncertainty intervals. Pooled weight summaries gave equal weight to each network because the matrices were equal-sized; their SEM used the three network-level summaries. Equal-network neuron–cycle distributions averaged separately normalized network fractions. These were descriptive means without uncertainty intervals. These summaries and fitted relationships do not treat neuron–cycle pairs as independent replicates.

Cycle-participation analysis excluded twelve of 18,000 presentations without a detected inhibitory burst, all from the seed-43 network at 27-ms inhibitory decay. These presentations contained excitatory activity, but it was too dispersed to recruit the inhibitory population. The affected images provided unusually weak input, while this network also had comparatively weak input connections and low firing rates. Because the same images recruited inhibition in the other network replicates, the failures appear to reflect a combination of weak stimuli, network-specific properties and stochastic encoding rather than a single cause.

All 84 included networks completed 50 training epochs.

P46

Spiking simulation, training and evaluation used custom Python/PyTorch code; all 84 classifiers were trained with PyTorch 2.11.0+cu128 on CUDA devices. Mean-field calculations used Python 3.10.19, NumPy 2.2.6 and SciPy 1.15.3. Experiment records retain stage-specific parameters, seeds, software/source revisions and original provenance for reused classifiers.