Abstract#

We asked how gamma oscillations begin in a population description of the PING circuit and which variables are essential. We swept drive, inhibitory timescale and noise, and compared the full rate model with simpler quasi-steady reductions.

Gamma emerged through oscillatory loss of stability; a reduced feedback system retained the oscillation, whereas the simplest reductions lost it. This supplies a mechanism for rhythmic recruitment, not proof that the trained spiking networks undergo the same bifurcation.

Results#

Hopf onset and frequency#

At the reference noise scale of 4 mV, one conjugate pair crossed at 0.594 nA with onset frequency 27.6 Hz. Up/down amplitudes nearly coincided, and the inhibitory-decay trend agreed qualitatively with separately measured spiking rhythms. We treat this as a candidate explanation of the exp025 — Accuracy and Firing Rate With and Without Inhibition, whose empirical marker is an inhibitory-rate crossing under input-weight scaling, not a fitted Hopf current. The model alone identifies neither that transition nor a minimum sustainable firing rate (Fig. 1).

Noise-scale onset sensitivity#

Across 3, 4, 5 and 6 mV, the sensitivity tests supported a reversible onset. The threshold spanned 0.59–0.75 nA while frequency remained near 27.6 Hz, and E amplitude at onset plus 0.4 nA fell from 8.5 to 3.7 Hz. At the 4 mV reference, equilibrium E/I rates were 4.1 and 0.7 Hz: the unstable equilibrium was low-rate, not silent. The noise parameter is free, so this sweep is not a calibration to spiking activity (Fig. 2).

Rates above Hopf onset#

The absolute cross-correlation peak lag was 5.2 ms. This magnitude is not a signed causal delay or a synaptic round-trip time (Fig. 3).

Four-state feedback-loop trajectories#

AMPA closely tracked E while the other variables showed larger phase offsets. These traces illustrate feedback timing; the ordering does not measure four independent transmission delays (Fig. 4).

Four-state trajectory projections#

The E–AMPA projection was narrow while other pairs enclosed larger areas. A periodic orbit is a curve and need not be planar; these projections do not demonstrate a centre manifold or prove that a particular pair closes the dynamics (Fig. 5).

Reduced-model responses#

The full model and AMPA-slaved three-variable reduction oscillated, whereas the rate-slaved two-variable probe decayed. Their onset frequencies were approximately 28 and 36 Hz, respectively. All six original-variable two-dimensional reductions had negative divergence. The three-variable minimum applies to this quasi-steady-state family, not to all possible models [1] (Fig. 6).

Methods#

We recomputed the deterministic population-rate model with 1.2/0.6-ms E/I refractory periods and a cancellation-resistant gain integral, then compared its onset frequencies with unchanged measurements from trained spiking networks.

Dataset#

Appendix: From spiking membranes to a population-rate closure#

Summary of COBA model.#

Here ${𝑉}_{𝑚}$ is membrane voltage (mV), ${𝐶}_{𝑚}$ capacitance (nF), ${𝑔}_{𝐿}$ leak conductance (µS), and ${𝐸}_{𝐿}$, ${𝐸}_{𝑒}$, ${𝐸}_{𝑖}$ are leak, excitatory and inhibitory reversal potentials (mV). ${𝑉}_{\text{th}}$ and ${𝑉}_{\text{reset}}$ are threshold and reset; $𝑠$ denotes a spike indicator in discrete time and a spike train in continuous time. $\mathrm{\Delta }{𝑡}_{\text{sim}}$ is the timestep, $𝑊$ a conductance increment per spike (µS), and superscripts identify the receiving E or I population.

The derivation starts from the COBANet model (exp100 — COBANet): conductance-based E and I membranes, a threshold-reset rule, and three exponential synapses (no E→E; I receives no inhibition):

Continuous-time form with tonic drive.#

Recast in continuous time. The synapses Equation 6 to Equation 8 are the exp-Euler form of first-order filters ${𝜏}_{\text{syn}}\dot{𝑔}=-𝑔+{𝜏}_{\text{syn}}\sum 𝑊𝑠$, used here as ODEs. At a constant input rate, ${𝑔}_{𝑒}^{𝐸}$ in Equation 6 settles to a steady mean, so its excitatory current into the E membrane Equation 3 is a near-constant depolarising drive; replacing it by a tonic current ${𝐼}_{\text{ext}}$ defines the swept control parameter. The E membrane then carries ${𝐼}_{\text{ext}}$ in place of ${𝑔}_{𝑒}^{𝐸}$:

Here ${𝑔}_{𝑖}^{𝐸}$ is the inhibition onto E and ${𝑔}_{𝑒}^{𝐼}$ the excitation onto I, each a continuous-time exponential filter of the presynaptic spikes:

with ${𝑠}^{𝐸}\left(𝑡\right),{𝑠}^{𝐼}\left(𝑡\right)$ the population spike trains and ${𝑊}^{𝐸𝐼},{𝑊}^{𝐼𝐸}$ the recurrent weight matrices; Equation 11 and Equation 12 are the continuous forms of Equation 7 and Equation 8.

Homogeneous coupling and population means.#

The motivating spiking network has ${𝑁}_{𝐸}=1024$ excitatory neurons and ${𝑁}_{𝐼}=256$ inhibitory neurons. These are population sizes, not extra state variables in the four-variable closure.

Now resolve the populations: index E neurons by $𝑗\in \left\{1,\dots ,{𝑁}_{𝐸}\right\}$ and I neurons by $𝑘\in \left\{1,\dots ,{𝑁}_{𝐼}\right\}$. The recurrent drive in Equation 11 and Equation 12 is the presynaptic sum ${𝑊}^{𝐸𝐼}{𝑠}^{𝐸}={\sum }_{𝑗}{𝑊}_{𝑘𝑗}^{𝐸𝐼}{𝑠}_{𝑗}^{𝐸}$ (and ${𝑊}^{𝐼𝐸}{𝑠}^{𝐼}={\sum }_{𝑘}{𝑊}_{𝑗𝑘}^{𝐼𝐸}{𝑠}_{𝑘}^{𝐼}$). Replace each random weight by its population mean, ${𝑊}_{𝑘𝑗}^{𝐸𝐼}\to {𝑤}^{𝐸𝐼}$ and ${𝑊}_{𝑗𝑘}^{𝐼𝐸}\to {𝑤}^{𝐼𝐸}$; the sums become

introducing the population-mean firing rates

A smooth-rate ansatz (short-window averaging) treats $𝐸\left(𝑡\right),𝐼\left(𝑡\right)$ as continuous, dropping weight heterogeneity and finite-size noise, the shot noise $\text{Var}\left[𝐸\left(𝑡\right)\right]\propto 𝐸\left(𝑡\right)/{𝑁}_{𝐸}$ for independent spike contributions at a fixed averaging window, with an analogous expression for I. Here $\text{Var}$ denotes variance across realizations. The corresponding typical fluctuation scale is $𝑂\left({𝑁}^{-1/2}\right)$, where $𝑁$ is population size. [(!) This scaling assumes independent or sufficiently weakly correlated contributions; recurrent synchrony can violate it. The earlier interpretation was that residual fluctuations at these finite population sizes smear onset and sustain weak noisy gamma below threshold. That remains a proposed explanation, not an effect isolated by this deterministic calculation.]

With no neuron index left, every E neuron sees the same ${𝑔}_{𝑖}^{𝐸}$ and every I neuron the same ${𝑔}_{𝑒}^{𝐼}$, collapsing the per-neuron conductances to population means. Defining lumped couplings

the conductance dynamics become

Two equations, down from ${𝑁}_{𝐸}+{𝑁}_{𝐼}$. (The fan-in scale $𝐺$ is converted to ${𝐽}_{𝐸\to 𝐼},{𝐽}_{𝐼\to 𝐸}$ in A.6.)

Running system, end of A.3: conductances are now two population means; the membrane is still per-neuron but sees those means:

Driving-force linearisation.#

The synaptic current is conductance times a driving force, $-𝑔\left({𝑉}_{𝑚}-{𝐸}_{\text{rev}}\right)$, a $𝑔$–${𝑉}_{𝑚}$ product, hence nonlinear. Freeze ${𝑉}_{𝑚}$ at rest, ${𝑉}_{\text{rest}}={𝐸}_{𝐿}=-65$ mV, in the driving force only (leak and threshold keep their full ${𝑉}_{𝑚}$-dependence, handled by the f-I curve in A.5). Each driving force becomes a fixed voltage gap:

The synaptic currents in Equation 9 and Equation 10 then lose their ${𝑉}_{𝑚}$-dependence and become proportional to conductance alone:

(inhibition pulls ${𝑉}_{𝑚}$ down, $\mathrm{\Delta }{𝑉}_{\text{inh,mag}}=+15$ mV; excitation pushes it up, $\mathrm{\Delta }{𝑉}_{\text{exc,mag}}=65$ mV). Removing the $𝑔$–${𝑉}_{𝑚}$ coupling reduces COBA to a current-based (CUBA) form; the cost is shunting: fixing ${𝑉}_{𝑚}$ neglects the fact that conductance also lowers the effective time constant (${𝜏}_{\text{eff}}={𝐶}_{𝑚}/{𝑔}_{\text{tot}}$, with ${𝑔}_{\text{tot}}$ the total membrane conductance).

Running system, end of A.4: the synaptic currents are now linear in conductance (no ${𝑉}_{𝑚}$ left in the driving force):

Population rate from an f-I curve.#

Under Equation 20 to Equation 23 the membrane equations Equation 9 and Equation 10 read ${𝐶}_{𝑚}\dot{{𝑉}_{𝑚}}=-{𝑔}_{𝐿}\left({𝑉}_{𝑚}-{𝐸}_{𝐿}\right)+{𝐼}_{\text{syn}}$, LIF with a synaptic current. A LIF neuron under constant net current $𝐼$ fires at its f-I rate $𝜑\left(𝐼\right)$; replacing each neuron’s spikes by that rate gives

with effective input currents (from Equation 9 and Equation 10 with Equation 20 to Equation 23 substituted)

(I receives only excitation; E receives the drive minus inhibitory current.) The instantaneous-rate replacement Equation 25 assumes slow inputs. The closure approximates the finite population response by relaxation on ${𝜏}_{𝐸}$ and ${𝜏}_{𝐼}$; this is [(!) a closure assumption, not an exact consequence of the single-neuron gain]:

where ${\mathrm{\Phi }}_{𝐸},{\mathrm{\Phi }}_{𝐼}$ are the smooth steady-state gain functions (the noisy LIF steady-state curve defined in Noisy LIF gain and parameter values). Two more equations down, together with Equation 17 and Equation 18, four equations in $\left(𝐸,𝐼,{𝑔}_{𝑒}^{𝐼},{𝑔}_{𝑖}^{𝐸}\right)$.

Running system, end of A.5: a closed 4D rate model in $\left(𝐸,𝐼,{𝑔}_{𝑒}^{𝐼},{𝑔}_{𝑖}^{𝐸}\right)$, constants not yet absorbed:

Absorb the driving-force constants.#

The prefactors $\mathrm{\Delta }{𝑉}_{\text{inh,mag}},\mathrm{\Delta }{𝑉}_{\text{exc,mag}}$ in Equation 28 and Equation 29 and the fan-in scalings in Equation 17 and Equation 18 are constants carrying no dynamics; fold them into the couplings:

[(!) Define current-valued coordinates ${ℎ}_{𝑒}^{𝐼}={𝑔}_{𝑒}^{𝐼}\mathrm{\Delta }{𝑉}_{\text{exc,mag}}$ and ${ℎ}_{𝑖}^{𝐸}={𝑔}_{𝑖}^{𝐸}\mathrm{\Delta }{𝑉}_{\text{inh,mag}}$ (nA). This invertible rescaling loses no dynamics. The figures show the original conductances $𝑔$ in µS; the equations below use $ℎ$ and current-valued couplings $𝐽$ (nA per spike).]

The 4D system#

After A.1–A.6, the mean-field equations are

in state $\left(𝐸,𝐼,{ℎ}_{𝑒}^{𝐼},{ℎ}_{𝑖}^{𝐸}\right)$. The tested quasi-steady reductions are examined in Appendix: Which variables can be eliminated? this is not a claim that four physical coordinates are the only possible description.

4D Jacobian#

At a fixed point $\left({𝐸}^{\ast },{𝐼}^{\ast },{ℎ}_{𝑒}^{𝐼\ast },{ℎ}_{𝑖}^{𝐸\ast }\right)$:

Here ${𝐽}_{\text{flow}}$ is the derivative of the vector field with respect to its state; ${\mathrm{\Phi }}_{𝐸}^{\prime },{\mathrm{\Phi }}_{𝐼}^{\prime }$ are gain derivatives with respect to input current, evaluated at the fixed-point arguments. Each linear mode evolves as ${𝑒}^{{𝜆}_{𝐽}𝑡}$, where ${𝜆}_{𝐽}$ is an eigenvalue: a negative real part decays and a positive real part grows. A simple Hopf requires one conjugate pair to cross with nonzero angular frequency while the other modes remain damped; criticality requires nonlinear information. Simultaneous crossings need a different analysis.

At the recorded crossing,

Here ${𝐼}_{\text{ext,Hopf}}$ is onset drive, ${𝜔}_{\text{Hopf}}$ angular frequency, and ${𝑓}_{\text{Hopf}}$ frequency in cycles per second. [(!) The factor 1000 converts milliseconds to seconds; it was missing from the earlier Hz equation.]

The eigenvalue plot can be read as a sequence of linear response tests. A point on the real axis is a non-oscillating mode; an off-axis conjugate pair oscillates while decaying to the left of the imaginary axis or growing to its right. The cyan crossing marks the change from damping to amplification. A double-Hopf has two simultaneously imaginary pairs; nearby nonlinear interactions can produce invariant tori with two angular phases [3]. [(!) This is general context, not an observed outcome here. The earlier assertion that a second pair crosses at higher drive is not supported by the recorded 0–4 nA sweep. Moreover, this model has

where $\text{tr}{𝐽}_{\text{flow}}$ is the sum of the four eigenvalues. Two simultaneously imaginary pairs would require zero trace, so a double-Hopf at an equilibrium is excluded for this particular four-filter model with finite positive time constants. This does not exclude every possible torus mechanism.]

Noisy LIF gain and parameter values#

Here $𝜇$ is mean input current (nA), ${𝜇}_{𝑉}$ its equivalent voltage (mV), $𝑢$ the dimensionless integration variable, and $\text{erf}$ the error function. $𝑄$ abbreviates the dimensionless integral; $𝑎$ and $𝑏$ are its scaled reset and threshold bounds. This is the same Siegert gain, split for readability. ${𝜏}_{𝑚}$ and ${𝜏}_{\text{ref}}$ are membrane and refractory times (ms); ${𝜎}_{𝑉}$ is the effective voltage-noise scale entering this formula, not a measured membrane standard deviation. The rate is in inverse milliseconds. The model uses ${𝐸}_{𝐿}={𝑉}_{\text{reset}}=-65$ mV, ${𝑉}_{\text{th}}=-50$ mV; E/I values are ${𝜏}_{𝑚}=\left(20,5\right)$ ms, ${𝑔}_{𝐿}=\left(0.05,0.10\right)$ µS and ${𝜏}_{\text{ref}}=\left(1.2,0.6\right)$ ms, respectively. The lumped conductance increments are ${𝐺}_{𝐸\to 𝐼}=1$ µS and ${𝐺}_{𝐼\to 𝐸}=2$ µS; fixed driving-force magnitudes are 65 and 15 mV. Fan-in-normalised mean weights give $𝐺=\overline{𝑤}𝑁$, where $𝑤$ is the mean presynaptic weight and $𝑁$ the number of presynaptic neurons. These values are inherited baseline settings, not fitted final-epoch weights.

In the spiking model’s strength notation,

Here ${𝐺}_{\text{ref}}=1$ µS is excitatory-to-inhibitory summed conductance and ${𝜌}_{\text{IE/EI}}=2$ is the dimensionless inhibitory/excitatory strength ratio. Fan-in normalization makes the individual mean weights $\frac{{𝐺}_{\text{ref}}}{{𝑁}_{𝐸}}$ and ${𝜌}_{\text{IE/EI}}\frac{{𝐺}_{\text{ref}}}{{𝑁}_{𝐼}}$. [(!) This restores the baseline parameter mapping; it does not identify these fixed couplings with the final trained weights.]

Quadrature used at most 200 subdivisions. For negative $𝑢$, we evaluated the integrand as $\text{erfcx}\left(-𝑢\right)=exp\left({𝑢}^2\right)\left(1+\text{erf}\left(𝑢\right)\right)$, where $\text{erfcx}$ is the scaled complementary error function; this avoids cancellation in $1+\text{erf}\left(𝑢\right)$. For nonnegative $𝑢$, the exponent remained capped at 700. The earlier calculation used 3/1.5-ms refractory periods and the cancelling expression, which produced integration warnings and percent-level gain errors at strong inputs. Both the refractory parameters and numerical evaluation changed in the new calculation; differences cannot be attributed solely to refractoriness. Independent evaluation of the gains at the recorded fixed points agreed within $9.1\times {10}^{-17}$ inverse milliseconds. The noise scale 3–6 mV was varied without calibration to spiking voltage statistics.

Appendix: Which variables can be eliminated?#

The four first-order filters preserve the sequence of excitation, recruitment and inhibition. Removing original variables by QSS substitution can destroy this feedback timing, but it is not the same operation as reducing the dynamics onto a centre manifold [1]. The following algebra retains the original reduction attempts, using the current-valued $ℎ$ coordinates from Absorb the driving-force constants.

1. Route A: the textbook Wilson-Cowan model (slave the conductances). The standard 2D tool is two rates with instantaneous coupling, the 4D model with instantaneous synaptic response at unchanged steady-state coupling. Slave each conductance to its filter’s steady value Equation 17 and Equation 18, ${ℎ}_{𝑒}^{𝐼}={𝜏}_{\text{AMPA}}{𝐽}_{𝐸\to 𝐼}𝐸$ and ${ℎ}_{𝑖}^{𝐸}={𝜏}_{\text{GABA}}{𝐽}_{𝐼\to 𝐸}𝐼$, and substitute into Equation 32:

   $${𝜏}_{𝐸}\dot{𝐸}=-𝐸+{\mathrm{\Phi }}_{𝐸}\left({𝐼}_{\text{ext}}-{𝜏}_{\text{GABA}}{𝐽}_{𝐼\to 𝐸}𝐼\right),$$

   (40)
   $${𝜏}_{𝐼}\dot{𝐼}=-𝐼+{\mathrm{\Phi }}_{𝐼}\left({𝜏}_{\text{AMPA}}{𝐽}_{𝐸\to 𝐼}𝐸\right).$$

   (41)

   Its divergence (the Jacobian trace),

   $$\frac{𝜕\dot{𝐸}}{𝜕𝐸}+\frac{𝜕\dot{𝐼}}{𝜕𝐼}=-\frac{1}{{𝜏}_{𝐸}}-\frac{1}{{𝜏}_{𝐼}}<0,$$

   (42)

   is a negative constant, so Bendixson–Dulac forbids a periodic orbit: no Hopf, for any drive or coupling on a simply connected region where the field is smooth. It has removed the two synaptic response lags; decay constants are filter response times, not fixed transmission delays. A controlled QSS approximation requires synapses fast relative to the full-model dynamics, not established here: ${𝜏}_{\text{GABA}}\approx 6$ ms is not negligible relative to the approximately $\frac{1000}{27.6}$ ms onset period.

2. Route B: quasi-steady-state the rates instead. The dual move: slave the rates, $𝐸={\mathrm{\Phi }}_{𝐸}\left({𝐼}_{\text{ext}}-{ℎ}_{𝑖}^{𝐸}\right)$ and $𝐼={\mathrm{\Phi }}_{𝐼}\left({ℎ}_{𝑒}^{𝐼}\right)$, into the conductance equations Equation 33, giving a 2D system in $\left({ℎ}_{𝑒}^{𝐼},{ℎ}_{𝑖}^{𝐸}\right)$:

   $${𝜏}_{\text{AMPA}}{\dot{ℎ}}_{𝑒}^{𝐼}=-{ℎ}_{𝑒}^{𝐼}+{𝜏}_{\text{AMPA}}{𝐽}_{𝐸\to 𝐼}{\mathrm{\Phi }}_{𝐸}\left({𝐼}_{\text{ext}}-{ℎ}_{𝑖}^{𝐸}\right),$$

   (43)
   $${𝜏}_{\text{GABA}}{\dot{ℎ}}_{𝑖}^{𝐸}=-{ℎ}_{𝑖}^{𝐸}+{𝜏}_{\text{GABA}}{𝐽}_{𝐼\to 𝐸}{\mathrm{\Phi }}_{𝐼}\left({ℎ}_{𝑒}^{𝐼}\right),$$

   (44)

   with the same negative-constant divergence,

   $$\frac{𝜕{\dot{ℎ}}_{𝑒}^{𝐼}}{𝜕{ℎ}_{𝑒}^{𝐼}}+\frac{𝜕{\dot{ℎ}}_{𝑖}^{𝐸}}{𝜕{ℎ}_{𝑖}^{𝐸}}=-\frac{1}{{𝜏}_{\text{AMPA}}}-\frac{1}{{𝜏}_{\text{GABA}}}<0,$$

   (45)

   so no cycle on a simply connected region with a smooth vector field: the displayed rate-slaved probe rings down (Figure 6). (These rates are the membrane variables, already reduced to an f-I rate.)

3. Route C: lump into fast and slow timescales. Slave the two fastest variables, the AMPA conductance (${𝜏}_{\text{AMPA}}=2$ ms) and the I rate (${𝜏}_{𝐼}=5$ ms), keeping the two slowest $\left\{𝐸,{ℎ}_{𝑖}^{𝐸}\right\}$ (${𝜏}_{\text{GABA}}=6$, ${𝜏}_{𝐸}=20$ ms):

   $${𝜏}_{𝐸}\dot{𝐸}=-𝐸+{\mathrm{\Phi }}_{𝐸}\left({𝐼}_{\text{ext}}-{ℎ}_{𝑖}^{𝐸}\right),$$

   (46)
   $${𝜏}_{\text{GABA}}{\dot{ℎ}}_{𝑖}^{𝐸}=-{ℎ}_{𝑖}^{𝐸}+{𝜏}_{\text{GABA}}{𝐽}_{𝐼\to 𝐸}{\mathrm{\Phi }}_{𝐼}\left({𝜏}_{\text{AMPA}}{𝐽}_{𝐸\to 𝐼}𝐸\right).$$

   (47)

   Trace $-1/{𝜏}_{𝐸}-1/{𝜏}_{\text{GABA}}<0$: no cycle. The split is forced anyway: the constants interleave, ${𝜏}_{\text{AMPA}}=2<{𝜏}_{𝐼}=5<{𝜏}_{\text{GABA}}=6<{𝜏}_{𝐸}=20$ ms, so “fast” and “slow” each mix a conductance with a rate.

4. All three fail for one structural reason. The network is a pure ring: the single loop $𝐸\to {ℎ}_{𝑒}^{𝐼}\to 𝐼\to {ℎ}_{𝑖}^{𝐸}\to 𝐸$, no recurrent E→E or I→I and no self-drive, so each variable’s only diagonal Jacobian term is its own decay and every gain ${\mathrm{\Phi }}^{\prime }$ sits off-diagonal. Eliminate any two variables and the 2D trace is $-1/{𝜏}_{𝑎}-1/{𝜏}_{𝑏}<0$, where ${𝜏}_{𝑎}$ and ${𝜏}_{𝑏}$ are the two remaining time constants; Bendixson–Dulac then rules out a cycle. Routes A–C are three of the $\left(\genfrac{}{}{0ex}{}{4}{2}\right)=6$ ways to pick the kept pair; the numerical study swept all six and none crossed. The negative-divergence argument applies to these original-variable QSS reductions. It does not exclude a nonlinear change of coordinates or a centre-manifold reduction of the same feedback mechanism.

   Adding recurrent E→E excitation or a cubic self-gain would change the diagonal dynamics and could escape this negative-divergence constraint. This motivated the earlier comparison with van der Pol/FitzHugh–Nagumo self-excitation oscillators. [(!) Such added terms would change the present physical-variable ring model, but they are not necessary for every two-dimensional representation of PING. A nonlinear centre-manifold reduction can describe the same local feedback dynamics without adding a physical E→E connection. The earlier universal claim about two-dimensional oscillators was too strong.]

5. Three dimensions survive. Slave only the fastest lag, the AMPA conductance ${ℎ}_{𝑒}^{𝐼}={𝜏}_{\text{AMPA}}{𝐽}_{𝐸\to 𝐼}𝐸$ (Figure 5), leaving a three-lag ring:

   $${𝜏}_{𝐸}\dot{𝐸}=-𝐸+{\mathrm{\Phi }}_{𝐸}\left({𝐼}_{\text{ext}}-{ℎ}_{𝑖}^{𝐸}\right),$$

   (48)
   $${𝜏}_{𝐼}\dot{𝐼}=-𝐼+{\mathrm{\Phi }}_{𝐼}\left({𝜏}_{\text{AMPA}}{𝐽}_{𝐸\to 𝐼}𝐸\right),$$

   (49)
   $${𝜏}_{\text{GABA}}{\dot{ℎ}}_{𝑖}^{𝐸}=-{ℎ}_{𝑖}^{𝐸}+{𝜏}_{\text{GABA}}{𝐽}_{𝐼\to 𝐸}𝐼,$$

   (50)

   This still Hopfs. Located like the 4D bifurcation (sweep ${𝐼}_{\text{ext}}$, diagonalise the $3\times 3$ Jacobian, find the complex-pair crossing), it gave ${𝐼}_{\text{ext,Hopf}}=0.69$ nA and ${𝑓}_{\text{Hopf}}=36$ Hz, both above the 4D values, in the numerical comparison. The six 2D reductions had no crossing in the sampled grid. The displayed probe compares 4D, 3D and the rate-slaved 2D model only (Figure 6); it is not a time-series panel of all six reductions.

6. Resolution: a centre manifold is a dynamical reduction, not a coordinate pair. [(!) Near a generic simple Hopf, the centre manifold is two-dimensional and tangent to the critical eigenspace; it need not be a plane. Its restricted vector field is a two-dimensional model of the local dynamics. A periodic orbit is a one-dimensional curve, and closed pairwise projections alone establish neither the manifold nor an autonomous two-variable closure. The three-variable minimum found here is restricted to the tested QSS ring family; it is compatible with the local two-dimensional centre-manifold description [1].]

The original dimensionality question remains useful: if the activity settles into a repeating rhythm, why keep four state variables? Amplitude and phase describe nearby oscillatory motion, while position projections show which coordinates nearly track each other. [(!) A closed loop need not lie in a plane, and a fixed periodic orbit itself needs only phase to locate a point. The nearly linear E–AMPA projection motivates testing AMPA slaving; it is not proof that the other coordinate pairs cannot parameterize a local manifold. The numerical QSS tests and the geometric question must be distinguished.]

Appendix: Numerical protocol and interpretation#

LSODA relative/absolute tolerances were $\left({10}^{-7},{10}^{-10}\right)$ for ramps, $\left({10}^{-9},{10}^{-12}\right)$ for waveforms, and $\left({10}^{-8},{10}^{-11}\right)$ for comparisons; maximum steps were 1, 0.25 and 0.5 ms, respectively. Brent refinement used absolute/relative tolerances $\left({10}^{-10},{10}^{-12}\right)$. The ramp began from the low-rate fixed point with a ${10}^{-3}$ ${\text{ms}}^{-1}$ E-rate kick; the waveform used the same kick. The 4D/2D comparison used a $2\times {10}^{-3}$ ${\text{ms}}^{-1}$ E kick. The common-drive ladder kicked E in 4D/3D by that amount and inhibitory conductance in 2D by $2\times {10}^{-3}$ µS. These probes are not matched perturbation-energy comparisons.

The spiking comparator comprised 18 independently trained networks: three seeds at each inhibitory decay of 4.5, 6, 9, 12, 18 and 27 ms. Final-epoch measurements used 1,000 fixed MNIST test images and 200 ms trials. Population traces were demeaned, full-trial Hann-window power spectra averaged over images, and a 5–150 Hz peak interpolated with its neighbouring bins before taking the three-network median. Onset eigenfrequencies and finite-drive spiking spectral peaks are different measurements. The comparison does not identify a causal contribution of gamma timing to classifier accuracy.

What the amplitude ramps test#

The measured amplitude is

where $𝐸\left(𝑡\right)$ is excitatory population rate and ${\mathrm{𝒯︀}}_{\text{obs}}$ is the final 500 ms observation window of each ramp step. Thus ${𝐴}_{\text{pp}}$ is peak-to-peak rate, in inverse milliseconds, not mean rate or oscillation power. Carrying each endpoint into the next drive tests whether the reached attractor depends on sweep direction.

In a supercritical Hopf, a stable small cycle emerges as the equilibrium loses stability. In a subcritical Hopf, the nearby cycle is unstable on the stable-equilibrium side and can bound its basin of attraction [3]. [(!) If a larger stable cycle also exists, a drive ramp can jump to it and remain there on reversal, producing a bistable hysteresis window. That larger cycle is an additional condition, not guaranteed by the local subcritical Hopf alone. Coincident sampled ramps support reversibility at their resolution; they do not establish absence of every unstable cycle.]

The supercritical normal form predicts the leading amplitude law

Here $\mathrm{\Delta }𝐼>0$ is excess drive (nA) and ${𝑐}_{\text{amp}}>0$ converts its square root into peak-to-peak rate; its units are ${\text{ms}}^{-1}$/$\sqrt{\text{nA}}$. The square-root law explains the formally unbounded onset slope in the ideal asymptotic description [3]. [(!) The recorded finite-range fit, with slope $1.1\times {10}^{-4}$ ${\text{ms}}^{-2}$/nA and ${𝑅}_{\text{fit}}^2=0.999$, is consistent with this law. It does not measure an infinite derivative or establish that most amplitude is acquired within a particular narrow drive band. Those earlier claims exceeded the sampled evidence.]

Mechanistic connections and their limits#

The motivating proposal was that the exp025 — Accuracy and Firing Rate With and Without Inhibition is a supercritical Hopf whose timescale comes from the E–I feedback loop. [(!) The present model supplies that candidate mechanism, but the empirical recruitment marker uses input-weight scaling and an inhibitory-rate crossing. No calibration equates it with this model’s current threshold. The equilibrium at onset had E/I rates 4.1/0.7 Hz and is not silent.]

The original frequency interpretation was that slower inhibition acts as the clock and that synchronous excitatory volleys sharpen the spiking rhythm, particularly at short inhibitory decay. [(!) The descending frequency trends support a timescale connection, but the spiking networks were separately retrained. Synchrony was not isolated as the cause of the mismatch, and the mean-field curve crosses above the spiking curve at 27 ms rather than remaining below throughout.]

The waveforms illustrate E recruitment of I followed by inhibition of E, with near-sinusoidal rates. The earlier account identified the measured 5.2 ms lag with E leading I and with a synaptic round trip. [(!) The recorded scalar is the absolute cross-correlation peak lag. It loses the sign and cannot by itself establish a causal or round-trip delay. AMPA and GABA decay times are filter response times, not fixed transmission delays. The loop interpretation remains a mechanism to test, not a delay measurement recovered from that scalar.]

References#

1. W. Zhang, V. Kirk, J. Sneyd, and M. Wechselberger. “Changes in the criticality of Hopf bifurcations due to certain model reduction techniques in systems with multiple timescales.” The Journal of Mathematical Neuroscience 1, 9 (2011). doi:10.1186/2190-8567-1-9
2. K. Kreutz-Delgado. “Mean Time-to-Fire for the Noisy LIF Neuron: A Detailed Derivation of the Siegert Formula.” arXiv (2015). doi:10.48550/arXiv.1501.04032
3. Y. A. Kuznetsov. Elements of Applied Bifurcation Theory, second edition. Springer (1998), sections 3.4, 5.2 and 8.6.