DMFT model

DMFT model#

exp033 · 28 May 2026 · pdf

Abstract#

A mean-field account of exp025′s recruitment cliff. In the 4D conductance DMFT the cliff is a Hopf bifurcation of the silent fixed point. With COBANet’s own LIF f-I curve as the population gain and couplings read off the biophysics (no fitted scale), the Jacobian eigenvalues place the threshold at $𝐼_{ext}^{*} = 0.6$ nA, the crossing pair’s imaginary part gives a gamma rhythm at $𝑓^{*} \approx 27.8$ Hz set by the synaptic timescales, and a hysteresis sweep shows the onset is supercritical: continuous and reversible.

Methods#

1. The 4D DMFT model#

1.1 Summary of COBA model.#

We start from the COBANet model (ar003 §2): conductance-based E and I membranes, a threshold-reset rule, and three exponential synapses (no E→E; I receives no inhibition):

𝐶_{𝑚}^{𝐸} {\dot{𝑉}}^{𝐸} = - 𝑔_{𝐿}^{𝐸} (𝑉^{𝐸} - 𝐸_{𝐿}) - 𝑔_{𝑒}^{𝐸} (𝑉^{𝐸} - 𝐸_{𝑒}) - 𝑔_{𝑖}^{𝐸} (𝑉^{𝐸} - 𝐸_{𝑖}) (1)

𝐶_{𝑚}^{𝐼} {\dot{𝑉}}^{𝐼} = - 𝑔_{𝐿}^{𝐼} (𝑉^{𝐼} - 𝐸_{𝐿}) - 𝑔_{𝑒}^{𝐼} (𝑉^{𝐼} - 𝐸_{𝑒}) (2)

𝑠_{𝑡 + 1} = 𝟙 [𝑉 \geq 𝑉_{th}], 𝑉 \leftarrow 𝑉_{reset} if 𝑠_{𝑡 + 1} = 1 or refractory (3)

𝑔_{𝑒, 𝑡 + 1}^{𝐸} = 𝑒^{- Δ 𝑡 / 𝜏_{AMPA}} 𝑔_{𝑒, 𝑡}^{𝐸} + 𝑊_{in} 𝑠_{𝑡}^{inp} (4)

𝑔_{𝑖, 𝑡 + 1}^{𝐸} = 𝑒^{- Δ 𝑡 / 𝜏_{GABA}} 𝑔_{𝑖, 𝑡}^{𝐸} + 𝑊_{ie} 𝑠_{𝑡}^{𝑖} (5)

𝑔_{𝑒, 𝑡 + 1}^{𝐼} = 𝑒^{- Δ 𝑡 / 𝜏_{AMPA}} 𝑔_{𝑒, 𝑡}^{𝐼} + 𝑊_{ei} 𝑠_{𝑡}^{𝑒} (6)

1.2. Continuous-time form with tonic drive.#

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

𝐶_{𝑚}^{𝐸} {\dot{𝑉}}^{𝐸} = - 𝑔_{𝐿}^{𝐸} (𝑉^{𝐸} - 𝐸_{𝐿}) - 𝑔_{𝑖}^{𝐸} (𝑉^{𝐸} - 𝐸_{𝑖}) + 𝐼_{ext} (7)

𝐶_{𝑚}^{𝐼} {\dot{𝑉}}^{𝐼} = - 𝑔_{𝐿}^{𝐼} (𝑉^{𝐼} - 𝐸_{𝐿}) - 𝑔_{𝑒}^{𝐼} (𝑉^{𝐼} - 𝐸_{𝑒}) (8)

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

𝜏_{AMPA} {\dot{𝑔}}_{𝑒}^{𝐼} = - 𝑔_{𝑒}^{𝐼} + 𝜏_{AMPA} 𝑊^{𝐸 𝐼} 𝑠^{𝐸} (𝑡) (9)

𝜏_{GABA} {\dot{𝑔}}_{𝑖}^{𝐸} = - 𝑔_{𝑖}^{𝐸} + 𝜏_{GABA} 𝑊^{𝐼 𝐸} 𝑠^{𝐼} (𝑡) (10)

with $𝑠^{𝐸} (𝑡), 𝑠^{𝐼} (𝑡)$ the population spike trains and $𝑊^{𝐸 𝐼}, 𝑊^{𝐼 𝐸}$ the recurrent weight matrices; (9)–(10) are the continuous forms of (5)–(6).

1.3. Homogeneous coupling and population means.#

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

\sum_{𝑗} 𝑊_{𝑘 𝑗}^{𝐸 𝐼} 𝑠_{𝑗}^{𝐸} ⟶ 𝑤^{𝐸 𝐼} \sum_{𝑗} 𝑠_{𝑗}^{𝐸} = 𝑤^{𝐸 𝐼} 𝑁_{𝐸} 𝐸 (𝑡) (11)

\sum_{𝑘} 𝑊_{𝑗 𝑘}^{𝐼 𝐸} 𝑠_{𝑘}^{𝐼} ⟶ 𝑤^{𝐼 𝐸} \sum_{𝑘} 𝑠_{𝑘}^{𝐼} = 𝑤^{𝐼 𝐸} 𝑁_{𝐼} 𝐼 (𝑡) (12)

introducing the population-mean firing rates

𝐸 (𝑡) \equiv \frac{1}{𝑁_{𝐸}} \sum_{𝑗 = 1}^{𝑁_{𝐸}} 𝑠_{𝑗}^{𝐸} (𝑡), 𝐼 (𝑡) \equiv \frac{1}{𝑁_{𝐼}} \sum_{𝑘 = 1}^{𝑁_{𝐼}} 𝑠_{𝑘}^{𝐼} (𝑡) . (13)

A smooth-rate ansatz (short-window averaging) treats $𝐸 (𝑡), 𝐼 (𝑡)$ as continuous, dropping weight heterogeneity and finite-size noise, the shot noise $Var [𝐸 (𝑡)] \propto 𝐸 (𝑡) / 𝑁_{𝐸}$ that vanishes only as $𝑁_{𝐸}, 𝑁_{𝐼} \to \infty$ . At finite $𝑁_{𝐸} = 1024, 𝑁_{𝐼} = 256$ the residual $𝑂 (𝑁^{- 1 / 2})$ fluctuations smear the Hopf onset and sustain a weak noisy gamma below threshold, both seen in the spiking simulations.

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

{\tilde{𝑊}}^{𝐸 𝐼} \equiv 𝑤^{𝐸 𝐼} 𝑁_{𝐸}, {\tilde{𝑊}}^{𝐼 𝐸} \equiv 𝑤^{𝐼 𝐸} 𝑁_{𝐼} (14)

the conductance dynamics become

𝜏_{AMPA} {\dot{𝑔}}_{𝑒}^{𝐼} = - 𝑔_{𝑒}^{𝐼} + 𝜏_{AMPA} {\tilde{𝑊}}^{𝐸 𝐼} 𝐸 (15)

𝜏_{GABA} {\dot{𝑔}}_{𝑖}^{𝐸} = - 𝑔_{𝑖}^{𝐸} + 𝜏_{GABA} {\tilde{𝑊}}^{𝐼 𝐸} 𝐼 (16)

Two equations, down from $𝑁_{𝐸} + 𝑁_{𝐼}$ . (The fan-in scale $\tilde{𝑊}$ folds into $𝑊^{𝐸 𝐼}, 𝑊^{𝐼 𝐸}$ in 1.6.)

Running system, end of 1.3: conductances are now two population means; the membrane is still per-cell but sees those means:

\begin{aligned} 𝐶_{𝑚}^{𝐸} {\dot{𝑉}}_{𝑗}^{𝐸} & = - 𝑔_{𝐿}^{𝐸} (𝑉_{𝑗}^{𝐸} - 𝐸_{𝐿}) - 𝑔_{𝑖}^{𝐸} (𝑉_{𝑗}^{𝐸} - 𝐸_{𝑖}) + 𝐼_{ext} \\ 𝐶_{𝑚}^{𝐼} {\dot{𝑉}}_{𝑘}^{𝐼} & = - 𝑔_{𝐿}^{𝐼} (𝑉_{𝑘}^{𝐼} - 𝐸_{𝐿}) - 𝑔_{𝑒}^{𝐼} (𝑉_{𝑘}^{𝐼} - 𝐸_{𝑒}) \\ 𝜏_{AMPA} {\dot{𝑔}}_{𝑒}^{𝐼} & = - 𝑔_{𝑒}^{𝐼} + 𝜏_{AMPA} {\tilde{𝑊}}^{𝐸 𝐼} 𝐸 \\ 𝜏_{GABA} {\dot{𝑔}}_{𝑖}^{𝐸} & = - 𝑔_{𝑖}^{𝐸} + 𝜏_{GABA} {\tilde{𝑊}}^{𝐼 𝐸} 𝐼 \end{aligned}

1.4. Driving-force linearisation.#

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

Δ 𝑉_{inh} \equiv 𝑉_{rest} - 𝐸_{𝑖} = - 65 - (- 80) = 15 mV (17)

Δ 𝑉_{exc} \equiv 𝑉_{rest} - 𝐸_{𝑒} = - 65 - 0 = - 65 mV

The synaptic currents in (7)–(8) then lose their $𝑉$ -dependence and become proportional to conductance alone:

- 𝑔_{𝑖}^{𝐸} (𝑉_{𝑗}^{𝐸} - 𝐸_{𝑖}) \approx - 𝑔_{𝑖}^{𝐸} Δ 𝑉_{inh} (18)

- 𝑔_{𝑒}^{𝐼} (𝑉_{𝑘}^{𝐼} - 𝐸_{𝑒}) \approx - 𝑔_{𝑒}^{𝐼} Δ 𝑉_{exc} = + 𝑔_{𝑒}^{𝐼} \cdot | 𝐸_{𝑒} - 𝑉_{rest} | (19)

(inhibition pulls $𝑉$ down, $Δ 𝑉_{inh} = + 15$ mV; excitation pushes it up, $| Δ 𝑉_{exc} | = 65$ mV). Removing the $𝑔$ – $𝑉$ coupling reduces COBA to a current-based (CUBA) form; the cost is shunting: with $𝑉$ fixed we ignore that conductance also lowers the effective time constant ( $𝜏_{eff} = 𝐶_{𝑚} / 𝑔_{tot}$ ).

Running system, end of 1.4: the synaptic currents are now linear in conductance (no $𝑉$ left in the driving force):

\begin{aligned} 𝐶_{𝑚}^{𝐸} {\dot{𝑉}}_{𝑗}^{𝐸} & = - 𝑔_{𝐿}^{𝐸} (𝑉_{𝑗}^{𝐸} - 𝐸_{𝐿}) - 𝑔_{𝑖}^{𝐸} Δ 𝑉_{inh} + 𝐼_{ext} \\ 𝐶_{𝑚}^{𝐼} {\dot{𝑉}}_{𝑘}^{𝐼} & = - 𝑔_{𝐿}^{𝐼} (𝑉_{𝑘}^{𝐼} - 𝐸_{𝐿}) + 𝑔_{𝑒}^{𝐼} | Δ 𝑉_{exc} | \\ 𝜏_{AMPA} {\dot{𝑔}}_{𝑒}^{𝐼} & = - 𝑔_{𝑒}^{𝐼} + 𝜏_{AMPA} {\tilde{𝑊}}^{𝐸 𝐼} 𝐸 \\ 𝜏_{GABA} {\dot{𝑔}}_{𝑖}^{𝐸} & = - 𝑔_{𝑖}^{𝐸} + 𝜏_{GABA} {\tilde{𝑊}}^{𝐼 𝐸} 𝐼 \end{aligned}

1.5. Population rate from an f-I curve.#

Under (17)–(19) the membrane equations (7)–(8) read $𝐶_{𝑚} \dot{𝑉} = - 𝑔_{𝐿} (𝑉 - 𝐸_{𝐿}) + 𝐼_{syn}$ , LIF with a synaptic current. A LIF cell under constant net current $𝐼$ fires at its f-I rate $𝜑 (𝐼)$ ; replacing each cell’s spikes by that rate gives

𝐸 (𝑡) \approx 𝜑_{𝐸} (𝐼_{eff}^{𝐸} (𝑡)), 𝐼 (𝑡) \approx 𝜑_{𝐼} (𝐼_{eff}^{𝐼} (𝑡)) (20)

with effective input currents (from (7)–(8) with (17)–(19) substituted)

𝐼_{eff}^{𝐸} (𝑡) = 𝐼_{ext} (𝑡) - 𝑔_{𝑖}^{𝐸} (𝑡) Δ 𝑉_{inh} (21)

𝐼_{eff}^{𝐼} (𝑡) = 𝑔_{𝑒}^{𝐼} (𝑡) | Δ 𝑉_{exc} | (22)

(I receives only excitation; E receives the drive minus the GABA shunt.) The instantaneous-rate replacement (20) holds only for slow inputs; in reality $𝐸 (𝑡)$ relaxes toward the f-I fixed point on $𝜏_{𝐸}$ , which we encode explicitly:

𝜏_{𝐸} \dot{𝐸} = - 𝐸 + Φ_{𝐸} (𝐼_{ext} - 𝑔_{𝑖}^{𝐸} Δ 𝑉_{inh}) (23)

𝜏_{𝐼} \dot{𝐼} = - 𝐼 + Φ_{𝐼} (𝑔_{𝑒}^{𝐼} | Δ 𝑉_{exc} |) (24)

where $Φ_{𝐸}, Φ_{𝐼}$ are the smooth steady-state gain functions (COBANet’s LIF f-I curve in the numerics; see Locating the Hopf). Two more equations down, together with (15)–(16), four equations in $(𝐸, 𝐼, 𝑔_{𝑒}^{𝐼}, 𝑔_{𝑖}^{𝐸})$ .

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

\begin{aligned} 𝜏_{𝐸} \dot{𝐸} & = - 𝐸 + Φ_{𝐸} (𝐼_{ext} - 𝑔_{𝑖}^{𝐸} Δ 𝑉_{inh}) \\ 𝜏_{𝐼} \dot{𝐼} & = - 𝐼 + Φ_{𝐼} (𝑔_{𝑒}^{𝐼} | Δ 𝑉_{exc} |) \\ 𝜏_{AMPA} {\dot{𝑔}}_{𝑒}^{𝐼} & = - 𝑔_{𝑒}^{𝐼} + 𝜏_{AMPA} {\tilde{𝑊}}^{𝐸 𝐼} 𝐸 \\ 𝜏_{GABA} {\dot{𝑔}}_{𝑖}^{𝐸} & = - 𝑔_{𝑖}^{𝐸} + 𝜏_{GABA} {\tilde{𝑊}}^{𝐼 𝐸} 𝐼 \end{aligned}

1.6. Absorb the driving-force constants.#

The prefactors $Δ 𝑉_{inh}, | Δ 𝑉_{exc} |$ in (23)–(24) and the fan-in scalings in (15)–(16) are constants carrying no dynamics; fold them into the couplings:

𝑊^{𝐸 𝐼} \equiv {\tilde{𝑊}}^{𝐸 𝐼} \cdot | Δ 𝑉_{exc} |, 𝑊^{𝐼 𝐸} \equiv {\tilde{𝑊}}^{𝐼 𝐸} \cdot Δ 𝑉_{inh} (25)

and absorb $| Δ 𝑉_{exc} |$ into the I-cell argument by redefining $𝑔_{𝑒}^{𝐼} \mapsto 𝑔_{𝑒}^{𝐼} \cdot | Δ 𝑉_{exc} |$ (similarly $𝑔_{𝑖}^{𝐸}$ ). The conductances now carry current units and the f-I curves take their argument directly: a change of variables, no dynamics lost.

1.7 The 4D system#

After 1.1–1.6, the mean-field equations are

𝜏_{𝐸} \dot{𝐸} = - 𝐸 + Φ_{𝐸} (𝐼_{ext} - 𝑔_{𝑖}^{𝐸}), 𝜏_{𝐼} \dot{𝐼} = - 𝐼 + Φ_{𝐼} (𝑔_{𝑒}^{𝐼}) (26)

𝜏_{AMPA} {\dot{𝑔}}_{𝑒}^{𝐼} = - 𝑔_{𝑒}^{𝐼} + 𝜏_{AMPA} 𝑊^{𝐸 𝐼} 𝐸, 𝜏_{GABA} {\dot{𝑔}}_{𝑖}^{𝐸} = - 𝑔_{𝑖}^{𝐸} + 𝜏_{GABA} 𝑊^{𝐼 𝐸} 𝐼 (27)

in state $(𝐸, 𝐼, 𝑔_{𝑒}^{𝐼}, 𝑔_{𝑖}^{𝐸})$ . Whether 4D is minimal, or a 2D reduction would do, is settled in the appendix, Is it really 4D?

1.8 4D Jacobian#

At a fixed point $(𝐸^{*}, 𝐼^{*}, 𝑔_{𝑒}^{𝐼 *}, 𝑔_{𝑖}^{𝐸 *})$ :

𝐽 = (\begin{matrix} - 1 / 𝜏_{𝐸} & 0 & 0 & - Φ_{𝐸}^{'} / 𝜏_{𝐸} \\ 0 & - 1 / 𝜏_{𝐼} & Φ_{𝐼}^{'} / 𝜏_{𝐼} & 0 \\ 𝑊^{𝐸 𝐼} & 0 & - 1 / 𝜏_{AMPA} & 0 \\ 0 & 𝑊^{𝐼 𝐸} & 0 & - 1 / 𝜏_{GABA} \end{matrix}) (28)

with $Φ_{𝐸}^{'}, Φ_{𝐼}^{'}$ evaluated at the fixed-point arguments.

2. Locating the Hopf#

2.1 The Hopf condition and sweep#

Below the cliff the network is silent: nudge it and the perturbation dies away. Above the cliff the same nudge grows into a sustained rhythm. The switch is a Hopf bifurcation of the silent fixed point: the smallest drive $𝐼_{ext}^{*}$ at which the fixed point stops damping oscillations and starts amplifying them.

Linear stability makes this precise. Each mode of the linearised system evolves as $𝑒^{𝜆 𝑡}$ , where $𝜆$ is an eigenvalue of the Jacobian (28): a negative real part decays, a positive one grows. The Hopf is the instant an oscillating mode (a complex-conjugate pair) crosses from the left half-plane to the right ( $Re 𝜆 = 0$ , $Im 𝜆 \neq 0$ ), so the silent state gives way to a growing oscillation rather than a static shift. To find it we sweep $𝐼_{ext}$ , track the fixed point, and diagonalise $𝐽$ at each step; $𝐼_{ext}^{*}$ is the first drive where the leading pair reaches the axis (Figure 2). The local analysis holds only if a single pair crosses while the rest stay damped (a simple Hopf); two pairs crossing at once (a double-Hopf) would seed tori the linearisation cannot see. Timescales: $𝜏_{𝐸} = 20$ , $𝜏_{𝐼} = 5$ , $𝜏_{AMPA} = 2$ , $𝜏_{GABA} = 6$ ms.

2.2 Gain and couplings from COBANet#

The gain $Φ$ and the couplings come from COBANet, not from a fit. For $Φ$ we use the LIF f-I curve: a cell under white-noise input of mean $𝜇$ and standard deviation $𝜎_{𝑉}$ fires at the Siegert rate

𝜑 (𝜇) = {[𝜏_{ref} + 𝜏_{𝑚} \sqrt{𝜋} \int_{(𝑉_{reset} - 𝜇_{𝑉}) / 𝜎_{𝑉}}^{(𝑉_{th} - 𝜇_{𝑉}) / 𝜎_{𝑉}} 𝑒^{𝑢^{2}} (1 + erf 𝑢) d 𝑢]}^{- 1}, 𝜇_{𝑉} = 𝐸_{𝐿} + 𝜇 / 𝑔_{𝐿} (29)

with COBANet’s $𝜏_{𝑚}, 𝑔_{𝐿}, 𝑉_{th}, 𝑉_{reset}, 𝜏_{ref}$ per population. The couplings are $𝑊^{𝐸 𝐼} = 𝑤^{𝐸 𝐼} 𝑁_{𝐸} | Δ 𝑉_{exc} |$ and $𝑊^{𝐼 𝐸} = 𝑤^{𝐼 𝐸} 𝑁_{𝐼} Δ 𝑉_{inh}$ (eqs 14, 25): the fan-in normalisation ( $𝑤 \mapsto 𝑤 / 𝑁$ ) makes $𝑤^{𝐸 𝐼} 𝑁_{𝐸}, 𝑤^{𝐼 𝐸} 𝑁_{𝐼}$ the ei-strength values $𝑠$ and $𝑟 𝑠$ (≈ 1 and 2 µS), times the driving forces (17). With $Φ$ in physical units $𝐼_{ext}$ is a current in nanoamps: the absolute scale is fixed by the biophysics, not chosen. The membrane-noise std $𝜎_{𝑉}$ is the one free parameter, set to 4 mV; the located Hopf is insensitive to it ( $𝑓^{*}$ moves under 1 Hz over $𝜎_{𝑉} \in [3, 6]$ mV).

3. Calculating its frequency#

At the crossing $𝜆 = \pm 𝑖 𝜔^{*}$ , so $𝑓^{*} = 𝜔^{*} / 2 𝜋$ is read from the Jacobian at threshold. To test whether inhibitory decay sets the period, we re-find the Hopf at each $𝜏_{GABA}$ of exp041′s retrained sweep and compare $𝑓^{*}$ to the spiking $𝑓_{𝛾}$ (Figure 3).

4. Classifying sub- or supercritical#

4.1 The hysteresis sweep#

Above the Hopf a limit cycle exists; whether it appears gently or abruptly decides whether exp025′s recruitment cliff is the graded, reversible onset that picture assumes. We test it with a hysteresis sweep: step $𝐼_{ext}$ quasi-statically up through $𝐼_{ext}^{*}$ and back down, at each step integrating (26)–(27) to steady state from the previous step’s end state (so any coexisting cycle is carried along), and record the peak-to-peak amplitude $𝐴 = {max}_{𝑡} 𝐸 - {min}_{𝑡} 𝐸$ .

A stable cycle born at threshold makes the rising and falling branches coincide, with $𝐴$ returning to zero at $𝐼_{ext}^{*}$ : a reversible, supercritical onset. An unstable cycle sitting below threshold instead acts as a basin boundary: the silent state jumps to a distant large-amplitude cycle that survives as the drive is lowered back, so the branches split into a hysteresis loop with a bistable window: subcritical (Figure 4).

4.2 The cycle waveform#

The same integration gives the cycle’s shape above onset: $𝐸$ and $𝐼$ are near-sinusoidal, with the E burst leading the I burst by the round-trip synaptic delay: E recruits I, I shunts E a few milliseconds later (Figure 5).

Results#

Three panels. A: the four Jacobian eigenvalues in the complex plane, coloured by drive, with one conjugate pair crossing the imaginary axis at the Hopf. B: the hysteresis sweep, rising and falling branches coinciding. C: gamma frequency falling with tau_GABA for the mean-field and the exp041 spiking network. — Figure 1: Summary of the analysis; the detailed panels follow. A, the Hopf: one eigenvalue pair crosses into the right half-plane at $𝐼^{*} = 0.6$ nA, fixing a gamma rhythm at $𝑓^{*} \approx 27.8$ Hz (full plot, Figure 2). B, the onset is supercritical and reversible, the up/down branches coinciding (Figure 4). C, the predicted frequency falls with $𝜏_{GABA}$ , tracking the spiking measurement (Figure 3). exp025′s recruitment cliff is a supercritical Hopf, set by the synaptic time constants, derived below from COBANet’s f-I curve and biophysical couplings with no fitted scale.

Frequency#

Sweeping $𝐼_{ext}$ and diagonalising $𝐽$ at each step locates the Hopf at

𝐼_{ext}^{*} = 0.6 nA, 𝜔^{*} = 0.175 rad/ms, 𝑓^{*} = \frac{𝜔^{*}}{2 𝜋} \approx 27.8 Hz,

in the PING gamma band, from COBANet’s f-I curve and biophysical couplings with no fitted scale. One complex pair crosses the imaginary axis while the others stay damped, a simple Hopf (Figure 2). The predictive content is the dependence on the synaptic timescales: across a $𝜏_{GABA}$ sweep $𝑓^{*}$ tracks the spiking network’s gamma qualitatively, not quantitatively (Figure 3).

The four 4D Jacobian eigenvalues in the complex plane across the drive sweep, coloured dark to bright with I_ext. One complex-conjugate pair lifts off the real axis and crosses the imaginary axis at plus/minus i omega-star; the other two eigenvalues stay in the left half-plane. — Figure 2: **How to read it.** Each dot is one eigenvalue $𝜆$ of the Jacobian (28). Its horizontal position is the growth rate $Re 𝜆$ (negative means that mode decays, positive means it grows), and the dotted line at $Re 𝜆 = 0$ is the stability boundary: everything to its left is damped. Its vertical position is the oscillation frequency $Im 𝜆$ (zero = a non-oscillating mode, on the real axis). The Jacobian is $4 \times 4$ , so each drive contributes four dots; sweeping $𝐼_{ext}$ (colour, dark→bright) drags them along the curves shown. Follow the upper and lower arms: one complex-conjugate pair lifts off the real axis and marches rightward, touching the boundary at $\pm 𝑖 𝜔^{*}$ (cyan circles) when $𝐼_{ext}^{*} = 0.6$ nA. That crossing is the Hopf, where the network tips from damping to amplifying, and the height of the crossing sets the rhythm, $𝑓^{*} = 𝜔^{*} / 2 𝜋 \approx 27.8$ Hz. The other two eigenvalues stay left of the line throughout, so a single oscillation goes unstable: a simple Hopf. (A second pair reaches the axis only at far higher drive, above the recruitment cliff.)

Gamma frequency versus tau_GABA. Both the calibrated mean-field f-star and the exp041 spiking f-gamma fall monotonically as tau_GABA increases; the mean-field curve runs below the spiking curve across the sweep, the gap largest at short tau_GABA and narrowing at long tau_GABA. — Figure 3: Both fall monotonically with $𝜏_{GABA}$ : the inhibitory decay is the clock in both. The match is qualitative, not quantitative. The calibrated mean-field is *flatter* and runs below the spiking network across the sweep, the gap largest at short $𝜏_{GABA}$ ( $27.8$ vs $43.9$ Hz at the canonical $6$ ms) and narrowing as $𝜏_{GABA}$ grows. The reduction captures the mechanism but not the full sensitivity, expected since it drops the spike-synchrony of the E-volley that sharpens the cycle most at short $𝜏_{GABA}$ .

Super or subcritical#

The rising and falling sweeps coincide (Figure 4): the amplitude grows continuously from zero at $𝐼_{ext}^{*}$ and retraces exactly on the way down, with no hysteresis (loop width $0$ nA, branch gap $2 \times 10^{- 6}$ ) and no cycle coexisting with the silent state. This is a supercritical Hopf, the recruitment cliff of exp025. The rising branch obeys the predicted $𝐴^{2} \propto (𝐼_{ext} - 𝐼_{ext}^{*})$ (slope $1.1 \times 10^{- 4}$ , $𝑅^{2} = 0.999$ ), so $𝐴 \propto \sqrt{𝐼_{ext} - 𝐼_{ext}^{*}}$ and $d 𝐴 / d 𝐼_{ext}$ diverges at threshold: most of the amplitude appears in a narrow band of drive above $𝐼_{ext}^{*}$ .

Peak-to-peak E amplitude versus I_ext for a quasi-static up-and-down ramp of the drive across I-star. The rising branch (drive increasing) and falling branch (drive decreasing) coincide, amplitude growing continuously from zero at the threshold with no hysteresis loop. — Figure 4: Quasi-static up/down ramp of the drive across $𝐼_{ext}^{*}$ . The rising branch (black, drive increasing) and falling branch (red, drive decreasing) coincide: the cycle turns on and off at the same drive, with amplitude growing continuously from zero. No loop, no bistable window: the reversible onset of a supercritical Hopf. (A subcritical onset would show the two branches separating into a hysteresis loop.)

Above onset, E leads I by ≈ 5.3 ms, the loop delay the 2D reduction omits.

The E rate (black) and I rate (red) over one limit cycle just above onset, on twin y-axes against time. Both are near-sinusoidal and the E burst leads the I burst by a few milliseconds. — Figure 5: $𝐸$ (black) and $𝐼$ (red) over the limit cycle at $𝐼_{ext} = 𝐼^{*} + 0.4$ nA. The E burst leads the I burst by ≈ 5.3 ms: E recruits I, I shunts E a few milliseconds later, and the cycle repeats, the loop delay the 2D reduction omits.

Appendix#

Is it really 4D?#

A Hopf gives a closed-loop trajectory, so the long-run motion is planar, so there should be a 2D description. This appendix collects the attempts. Verdict: a 2D description exists (a centre manifold), but the vector field does not reduce below three.

A Hopf gives planar (two-dimensional) dynamics, so 4D looks suspect. A limit cycle is a closed loop, traversed with two coordinates (amplitude and phase), so the motion lies in a plane. If it is 2D, a 2D model should exist; steps 3–5 hunt for it.
The geometry confirms the motion is 2D (Figures 6–7). The four variables cycle in loop order ( $𝐸$ , then $𝑔_{𝑒}^{𝐼}$ , then $𝐼$ , then $𝑔_{𝑖}^{𝐸}$ , each lagging the last, Figure 6). Projected onto every variable pair (Figure 7), the cycle is one closed loop on a thin 2D sheet, a centre manifold (the loop is a 1D curve on it). Only $(𝐸, 𝑔_{𝑒}^{𝐼})$ nearly collapses to a line (AMPA trails the E rate); the rest enclose real area, so no variable pair can stand in for the sheet.
Route A: the textbook Wilson-Cowan model (slave the conductances). The standard 2D tool is two rates with instantaneous coupling, the 4D model with infinitely fast synapses. Slave each conductance to its filter’s steady value (15)–(16), $𝑔_{𝑒}^{𝐼} = 𝜏_{AMPA} 𝑊^{𝐸 𝐼} 𝐸$ and $𝑔_{𝑖}^{𝐸} = 𝜏_{GABA} 𝑊^{𝐼 𝐸} 𝐼$ , and substitute into (26):
$𝜏_{𝐸} \dot{𝐸} = - 𝐸 + Φ_{𝐸} (𝐼_{ext} - 𝜏_{GABA} 𝑊^{𝐼 𝐸} 𝐼), (30)$ $𝜏_{𝐼} \dot{𝐼} = - 𝐼 + Φ_{𝐼} (𝜏_{AMPA} 𝑊^{𝐸 𝐼} 𝐸) . (31)$
Its divergence (the Jacobian trace),
$\frac{𝜕 \dot{𝐸}}{𝜕 𝐸} + \frac{𝜕 \dot{𝐼}}{𝜕 𝐼} = - \frac{1}{𝜏_{𝐸}} - \frac{1}{𝜏_{𝐼}} < 0, (32)$
is a negative constant, so Bendixson–Dulac forbids a periodic orbit: no Hopf, for any drive or coupling. It has dropped the round-trip delay (E excites I after $𝜏_{AMPA}$ , I shunts E after $𝜏_{GABA}$ ): zero-lag inhibition damps but cannot overshoot. The reduction is valid only if synapses are fast relative to the rhythm, which they are not ( $𝜏_{GABA} \approx 6$ ms is the gamma period’s order).
Route B: quasi-steady-state the rates instead. The dual move: slave the rates, $𝐸 = Φ_{𝐸} (𝐼_{ext} - 𝑔_{𝑖}^{𝐸})$ and $𝐼 = Φ_{𝐼} (𝑔_{𝑒}^{𝐼})$ , into the conductance equations (27), giving a 2D system in $(𝑔_{𝑒}^{𝐼}, 𝑔_{𝑖}^{𝐸})$ :
$𝜏_{AMPA} {\dot{𝑔}}_{𝑒}^{𝐼} = - 𝑔_{𝑒}^{𝐼} + 𝜏_{AMPA} 𝑊^{𝐸 𝐼} Φ_{𝐸} (𝐼_{ext} - 𝑔_{𝑖}^{𝐸}), (33)$ $𝜏_{GABA} {\dot{𝑔}}_{𝑖}^{𝐸} = - 𝑔_{𝑖}^{𝐸} + 𝜏_{GABA} 𝑊^{𝐼 𝐸} Φ_{𝐼} (𝑔_{𝑒}^{𝐼}), (34)$
with the same negative-constant divergence,
$\frac{𝜕 {\dot{𝑔}}_{𝑒}^{𝐼}}{𝜕 𝑔_{𝑒}^{𝐼}} + \frac{𝜕 {\dot{𝑔}}_{𝑖}^{𝐸}}{𝜕 𝑔_{𝑖}^{𝐸}} = - \frac{1}{𝜏_{AMPA}} - \frac{1}{𝜏_{GABA}} < 0, (35)$
so no cycle: it rings down (Figure 8). (These rates are the membrane variables, already reduced to an f-I rate.)
Route C: lump into fast and slow timescales. Slave the two fastest variables, the AMPA conductance ( $𝜏_{AMPA} = 2$ ms) and the I rate ( $𝜏_{𝐼} = 5$ ms), keeping the two slowest ${𝐸, 𝑔_{𝑖}^{𝐸}}$ ( $𝜏_{GABA} = 6$ , $𝜏_{𝐸} = 20$ ms):
$𝜏_{𝐸} \dot{𝐸} = - 𝐸 + Φ_{𝐸} (𝐼_{ext} - 𝑔_{𝑖}^{𝐸}), (36)$ $𝜏_{GABA} {\dot{𝑔}}_{𝑖}^{𝐸} = - 𝑔_{𝑖}^{𝐸} + 𝜏_{GABA} 𝑊^{𝐼 𝐸} Φ_{𝐼} (𝜏_{AMPA} 𝑊^{𝐸 𝐼} 𝐸) . (37)$
Trace $- 1 / 𝜏_{𝐸} - 1 / 𝜏_{GABA} < 0$ : no cycle. The split is forced anyway: the constants interleave, $𝜏_{AMPA} = 2 < 𝜏_{𝐼} = 5 < 𝜏_{GABA} = 6 < 𝜏_{𝐸} = 20$ ms, so “fast” and “slow” each mix a conductance with a rate.
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 $Φ^{'}$ sits off-diagonal. Eliminate any two variables and the 2D trace is $- 1 / 𝜏_{𝑎} - 1 / 𝜏_{𝑏} < 0$ ; Bendixson–Dulac then rules out a cycle. Routes A–C are three of the $(\binom{4}{2}) = 6$ ways to pick the kept pair; the runner sweeps all six and none crosses. A 2D model that does oscillate has added a destabiliser PING lacks: recurrent excitation or a cubic self-gain, i.e. a positive diagonal term. We do not add one: PING has no $𝑊^{𝐸 𝐸}$ to supply it, and bolting one on would make a self-excitation oscillator (van der Pol / FitzHugh–Nagumo), not PING: a different mechanism for the gamma, not this network’s.
Three dimensions survive. Slave only the fastest lag, the AMPA conductance $𝑔_{𝑒}^{𝐼} = 𝜏_{AMPA} 𝑊^{𝐸 𝐼} 𝐸$ (step 2′s near-degenerate pair), leaving a three-lag ring:
$𝜏_{𝐸} \dot{𝐸} = - 𝐸 + Φ_{𝐸} (𝐼_{ext} - 𝑔_{𝑖}^{𝐸}), (38)$ $𝜏_{𝐼} \dot{𝐼} = - 𝐼 + Φ_{𝐼} (𝜏_{AMPA} 𝑊^{𝐸 𝐼} 𝐸), (39)$ $𝜏_{GABA} {\dot{𝑔}}_{𝑖}^{𝐸} = - 𝑔_{𝑖}^{𝐸} + 𝜏_{GABA} 𝑊^{𝐼 𝐸} 𝐼, (40)$
This still Hopfs. Located like the 4D bifurcation (sweep $𝐼_{ext}$ , diagonalise the $3 \times 3$ Jacobian, find the complex-pair crossing), it gives $𝐼_{ext}^{*} = 0.7$ nA and $𝑓^{*} = 37$ Hz, both above the 4D values, since dropping the AMPA lag stiffens the loop. The same sweep finds no crossing for either 2D reduction. Figure 8: 4D and 3D sustain the rhythm; all three 2D reductions ring down.
Resolution: the 2D description is the centre manifold, not a coordinate pair. The 2D description that exists is the centre manifold, the plane of the two critical, oscillatory eigenvectors, tilted across all of $(𝐸, 𝐼, 𝑔_{𝑒}^{𝐼}, 𝑔_{𝑖}^{𝐸})$ : a plane in the eigenbasis, not any coordinate pair. So it is genuinely 2D yet unreachable by dropping two physical variables: a coordinate projection lands on the wrong plane, where the gains go off-diagonal and Bendixson–Dulac kills the cycle (steps 3–6). The attractor is a 1D loop on a 2D manifold; the vector field does not reduce below three, since a Hopf needs three first-order lags in the $𝐸 \to 𝐼 \to 𝐸$ ring to build the destabilising phase. 4D is the natural model, 3D the floor, and the 2D centre manifold where the cycle lives, not a model in the original variables.

Four stacked time series over one limit cycle sharing a time axis, in loop order E rate, g_e^I, I rate, g_i^E. Each variable peaks after the one above it, one round trip of the ring per gamma cycle; g_e^I tracks the E rate almost rigidly. — Figure 6: The four state variables over the limit cycle ( $𝐼_{ext} = 𝐼^{*} + 0.4$ nA), in loop order $𝐸 \to 𝑔_{𝑒}^{𝐼} \to 𝐼 \to 𝑔_{𝑖}^{𝐸}$ . Each peaks after the one above it: $𝑔_{𝑒}^{𝐼}$ tracks the $𝐸$ rate almost rigidly (the near-degenerate pair of step 2), then drives $𝐼$ , which fills $𝑔_{𝑖}^{𝐸}$ , which shunts $𝐸$ : one round trip of the ring per gamma cycle.

The 4D limit cycle projected onto all six variable pairs, one closed loop per panel. The E versus g_e^I panel collapses almost to a line; the other five pairs enclose genuine area. — Figure 7: The 4D limit cycle ( $𝐼_{ext} = 𝐼^{*} + 0.4$ nA) projected onto all six variable pairs. Every projection is one closed loop: the trajectory lives on a 2D centre manifold. The $(𝐸, 𝑔_{𝑒}^{𝐼})$ panel collapses almost to a line: the fast AMPA conductance tracks the E rate, which is why slaving it (the 3D reduction) preserves the dynamics, while the other pairs enclose genuine area and cannot be collapsed.

The g_i^E deviation from its fixed point after a small kick at a common drive, for three models. The 4D full model (black solid) and the 3D AMPA-slaved reduction (black dashed) both sustain a limit cycle; the 2D rate-slaved reduction (red) rings down to the fixed point. — Figure 8: The reductions tested by simulation: $𝑔_{𝑖}^{𝐸}$ after a small kick at a common drive ( $𝐼_{ext} = 1$ nA, above every threshold). The 4D model (black, $𝑓^{*} = 28$ Hz) and the 3D AMPA-slaved reduction (black dashed, $𝑓^{*} = 37$ Hz) both sustain the limit cycle; the 2D rate-slaved reduction (red) rings down to its fixed point: no Hopf, as eq (35) requires.