033 — Mean field model

Abstract

nb025 characterised the PING rate floor empirically: a recruitment cliff where the I-loop engages, gamma at ≈25 Hz. Both are bifurcation behaviours, and this entry pins them down in a mean-field model — with three goals and no more: locate the Hopf, predict its frequency, and classify it super- or subcritical. The right object is the 4D conductance mean-field $(E, I, g_e^I, g_i^E)$: the Hopf is a linear fact about its $4\times4$ Jacobian (giving $f^\star \approx 24.5$ Hz, in the spiking network’s gamma band and tracking its dependence on $\tau_\text{GABA}$), and direct simulation classifies it supercritical. The tempting 2D rate model cannot host the rhythm at all — provably, by Bendixson–Dulac, and because the reduction that produces it is not even asymptotically valid — so the analysis stays in 4D throughout.

Methods

Why a 2D model is the wrong object

The textbook tool for a bifurcation in an E/I circuit is a 2D Wilson-Cowan rate model: two population rates $(E, I)$ with sigmoid gains $\Phi_{E,I}$ and instantaneous E↔I coupling. The PING architecture carries no recurrent E→E or I→I synapse ($W^{EE} = W^{II} = 0$), so the model is

It is the wrong object, for two reasons.

The reduction isn’t valid. Collapsing the full conductance-based dynamics to (1) means adiabatically eliminating the synaptic conductances — legitimate only if they relax fast relative to the rhythm. They do not: $\tau_\text{GABA} \approx 5$–$10$ ms is the same order as the gamma period itself ($\approx 25$ ms at 40 Hz). The inhibitory decay is not a slaved fast variable, it is the clock that sets the period. With no timescale separation to exploit, the elimination that produces (1) is unjustified — (1) is not a simplification of PING, it is a different system.

And even if you force it, it cannot oscillate. The divergence of the planar field (1) is pinned to a strictly negative constant,

because $\Phi_E$‘s argument carries no $E$ and $\Phi_I$‘s carries no $I$ — with no E→E or I→I synapse the gain $\Phi'$ multiplies only the off-diagonal cross-terms (it sets the rotation rate, never the divergence). By the Bendixson–Dulac theorem a planar flow whose divergence holds one sign on a simply connected region encloses no periodic orbit; so (1) sustains no limit cycle for any drive, coupling strength, gain shape, or amplitude — not merely no Hopf, but no oscillation anywhere in phase space. (The local face of the same fact: the Jacobian trace of (1) equals (2) everywhere, so it never reaches the $\mathrm{tr}(J) = 0$ a Hopf requires.) Consequently any 2D rate model that does oscillate has to import a destabiliser PING lacks — recurrent excitation in Wilson-Cowan, a cubic self-gain in FitzHugh-Nagumo — so it reproduces the rhythm through a fictional mechanism rather than approximating the real one.

What’s missing is the delay. What actually destabilises the silent state is the round-trip phase lag of the loop — E excites I after $\tau_\text{AMPA}$, I shunts E after $\tau_\text{GABA}$ — and (1), with instantaneous coupling, has discarded exactly it: inhibition that tracks $E$ with zero delay can only damp, never overshoot. Holding the lag is irreducibly more than two-dimensional — either keep the two synaptic conductances $(g_e^I, g_i^E)$ as first-order states (the 4D model below, whose synaptic poles rotate a complex pair across the imaginary axis at a Hopf), or write the delay explicitly (a delay-differential equation, whose transcendental characteristic quasi-polynomial carries infinitely many roots). So the analysis stays in 4D. “Local analysis near the Hopf” is done in the 4D coordinates — the Jacobian crossing below, and the amplitude law of §3.2 — never by swapping in a planar system.

Figure 1. Same drive, same kick — only the 4D model oscillates

The 4D mean-field model

We start with $N_E + N_I$ coupled conductance-based spiking neurons and reduce them, step by step, to four coupled ODEs for the population-mean rates $E(t), I(t)$ and the population-mean conductances $g_e^I(t), g_i^E(t)$. Each step replaces one degree of freedom with a smoother, lower-dimensional approximation; the cost of each is recorded at the end. The full chain is:

$N$-cell spiking COBA → homogeneous-coupling spike-train integrals → population-mean rates and conductances → driving-force linearisation → relaxational rate dynamics with f-I gain $\Phi$ → final 4D system in $(E, I, g_e^I, g_i^E)$.

Step 1. Single-cell COBA equations. Index E cells by $j \in \{1, \ldots, N_E\}$ and I cells by $k \in \{1, \ldots, N_I\}$. From the methods paper ar009 eq. (1) with the architecture constraints $W^{EE} = W^{II} = 0$ (no recurrent E→E, no I→I) and an external drive $I_{\text{ext},j}$ added on E:

The E cell has no $g_e^E$ term because the only excitation it receives in this analysis is $I_\text{ext}$ (the feedforward drive); the I cell has no $g_i^I$ because there is no I→I coupling. The conductances $g_{i,j}^E, g_{e,k}^I$ are themselves driven by presynaptic spikes through first-order linear filters:

with $s_j^E(t) = \sum_n \delta(t - t_{j,n}^E)$ the spike train of E-cell $j$ as a sum of Dirac pulses at its firing times (similarly for $s_k^I$). (D2a–b) are the continuous-time forms of eq. (8) and (7) in the methods backbone — exponential synapses with per-spike weight $W$ and decay $\tau$.

Step 2. Homogeneous coupling. The weight matrices in the simulator are random with finite variance. Replace each entry by its population mean: $W^{EI}_{kj} \to w^{EI}$ and $W^{IE}_{jk} \to w^{IE}$ for all $j, k$. The presynaptic-summation terms in (D2a–b) become

where we introduce the population-mean firing rates

Strictly, $E(t)$ and $I(t)$ as defined are sums of deltas. The mean-field interpretation is that under a smooth-rate ansatz — averaging over a short window or assuming a continuous spiking density — we treat $E, I$ as smooth functions, so $w^{EI} N_E E(t)$ becomes an ordinary continuous source term. (D3) drops two things: weight heterogeneity, and finite-size fluctuations. The latter is the stochastic noise that survives in any actual finite population: even with all weights at their mean values, the population rate $E(t) = N_E^{-1}\sum_j s_j^E(t)$ is a sum of $N_E$ independent point-process samples, so it has shot-noise variance $\mathrm{Var}[E(t)] \propto E(t)/N_E$ over a short window. The mean-field limit takes $N_E, N_I \to \infty$ at fixed $\tilde W^{EI}, \tilde W^{IE}$, in which the law of large numbers makes $E(t), I(t)$ converge to their expected values and the variance vanishes. At finite $N_E = 1024, N_I = 256$ the residual $O(N^{-1/2})$ fluctuations act as a noise term on (D4), which (a) smears the Hopf bifurcation — instead of a sharp $A^2 \propto (I_\text{ext} - I_\text{ext}^\star)$ onset, the spiking simulator shows a rounded shoulder with sub-threshold variance — and (b) maintains a low-amplitude noisy gamma oscillation just below criticality where the linearised mean-field predicts strict silence. Both effects reappear in the spiking simulations, with magnitude controlled by $1/\sqrt{N}$.

Because the source terms (D3a) and (D3b) no longer carry a cell index, every E cell sees the same $g_i^E$ and every I cell sees the same $g_e^I$. The per-cell conductance equations collapse to two population-mean equations. Defining lumped couplings

(carries-the-fan-in: each E spike contributes the per-edge weight, and there are $N_E$ presynaptic E cells), the conductance dynamics become

Two equations down from $N_E + N_I$. (Note: the absolute scale of $\tilde W^{EI}, \tilde W^{IE}$ folds into $W^{EI}, W^{IE}$ used in the final 4D system after Step 5.)

Step 3. Driving-force linearisation. The synaptic terms in (D1) still depend on $V$ through factors $(V - E_I)$ and $(V - E_E)$. Take $V$ to its resting value in the synaptic terms only (not in the leak or the threshold-crossing, which we’ll handle separately via the f-I curve). With $V_\text{rest} = -65$ mV and the reversal potentials from ar009,

(so $-(V - E_E) \approx -\Delta V_\text{exc} = +65$ mV — depolarising drive on I). The synaptic currents in (D1a–b) become

— each synaptic input is now linear in conductance, with the cell-population-independent prefactor folded into a constant. This is the crucial reduction: it kills the multiplicative coupling between membrane potential and synaptic conductance that makes COBA so different from CUBA, but only in the driving force. Shunting (the effect of $g$ on $\tau_\text{eff}$ and $g_\text{tot}$) is dropped at this step; this is what the table at the end records.

Step 4. Population rate from an f-I curve. Under (D5) the membrane equations (D1a–b) have the form $C_m \dot V = -g_L (V - E_L) + I_\text{syn}$ — exactly LIF with a synaptic current. A LIF cell receiving constant net input current $I$ fires at rate $\phi(I)$ where $\phi$ is its f-I curve. Replacing each cell’s spike output by its rate response gives

with effective input currents (from D1 with D5 substituted)

(I cells receive only excitation in this architecture; E cells receive external drive minus the GABA shunt.) The instantaneous-rate replacement (D6a) is only valid when input changes are slow compared to the membrane integration time. In reality, $E(t)$ relaxes toward the f-I-curve fixed point on its membrane time constant $\tau_E$. Encode this with an explicit relaxation:

where $\Phi_E, \Phi_I$ are the smooth steady-state gain functions (sigmoids in the numerics). Drives that change faster than $\tau_E, \tau_I$ are filtered; the rate doesn’t snap to the f-I curve instantaneously. Two more equations down — together with (D4a–b) we now have four equations in $(E, I, g_e^I, g_i^E)$.

Step 5. Absorb the driving-force constants. The fixed prefactors $\Delta V_\text{inh}, |\Delta V_\text{exc}|$ in (D7) and the fan-in scalings in (D4) are just multiplicative constants — they don’t carry dynamical information. Fold them into the couplings:

and likewise absorb the $|\Delta V_\text{exc}|$ inside the I-cell argument by redefining $g_e^I \mapsto g_e^I \cdot |\Delta V_\text{exc}|$ (similarly for $g_i^E$). The conductances now have units of input current rather than µS; the f-I curves take their argument directly. This change of variables is pure bookkeeping — no dynamics are lost.

The four equations after (D4) + (D7) + (D8) are exactly the 4D system shown next.

The 4D system

After Steps 1–5, the mean-field equations are

in state $(E, I, g_e^I, g_i^E)$. Whether this 4D state is minimal — or a 2D reduction would do — was settled in Why 2D reductions fail: the synaptic phase delays $\tau_\text{AMPA}, \tau_\text{GABA}$ are load-bearing, and a planar model that eliminates them cannot oscillate. The 4D state keeps them.

4D Jacobian

At a fixed point $(E^\star, I^\star, g_e^{I\star}, g_i^{E\star})$:

with $\Phi'_E, \Phi'_I$ evaluated at the fixed-point arguments.

Results

Of the three goals, two are linear — locating the Hopf and reading off its frequency are facts about the $4\times4$ Jacobian — and one is nonlinear — whether the bifurcation is super- or subcritical. We take them in that order.

Frequency

We locate the bifurcation numerically: sweep $I_\text{ext}$, track the fixed point (near-silent until the loop engages), and diagonalise $J$ at each step. The recruitment cliff is a Hopf bifurcation — the smallest $I_\text{ext}^\star$ at which a complex-conjugate eigenvalue pair reaches the imaginary axis ($\mathrm{Re}(\lambda) = 0$, $\mathrm{Im}(\lambda) \neq 0$), so the silent fixed point loses stability to a growing oscillation rather than to a static shift. The crossing frequency $\omega^\star = \mathrm{Im}(\lambda)$ is the gamma frequency, $f^\star = \omega^\star / 2\pi$.

For the canonical parameters ($\tau_E = 20$, $\tau_I = 5$, $\tau_\text{AMPA} = 2$, $\tau_\text{GABA} = 9$ ms; $W^{EI} = 80$, $W^{IE} = 60$):

That $f^\star \approx 24.5$ Hz sits squarely in the PING gamma band. It is not a parameter-free number — the mean-field’s abstract coupling units set the overall scale — but the dependence on the synaptic timescales is the predictive content, and across a $\tau_\text{GABA}$ sweep it tracks the spiking network’s gamma (Figure 3 below).

Tracking all four eigenvalues does double duty. At onset exactly one complex pair reaches the imaginary axis while the other modes stay well damped, so this is a simple, generic Hopf — not a double-Hopf where two pairs cross together (which would spawn tori the linear frequency cannot capture). That is precisely the condition under which a local near-onset analysis is valid. The complex plane shows the crossing and reads off the frequency at once:

Figure 2. 4D eigenvalues in the complex plane

Does the frequency track the synaptic clock? The mechanistic claim is that $\tau_\text{GABA}$ sets the period. Re-finding the Hopf at each $\tau_\text{GABA}$ in nb041’s retrained sweep and overlaying the spiking network’s measured $f_\gamma$:

Figure 3. Frequency vs inhibitory decay — mean-field vs spiking

Which coupling moves the cliff? Sweeping the recruitment threshold $I_\text{ext}^\star$ across the $(W^{EI}, W^{IE})$ plane recovers nb025’s empirical asymmetry from the equations:

Figure 4. Recruitment threshold across the coupling plane

Super or subcritical

The eigenvalue crossing says a limit cycle is born at $I_\text{ext}^\star$ but not whether it appears with zero amplitude (supercritical, soft onset) or with a finite jump and hysteresis (subcritical, hard onset). That is a nonlinear question, classically the sign of the first Lyapunov coefficient of the Hopf normal form. We read it off the full 4D field by direct simulation instead: integrate (3)–(4) from a small perturbation of the fixed point across a grid of $I_\text{ext}$ straddling $I_\text{ext}^\star$ and measure the steady-state peak-to-peak $E$ amplitude — the full top-to-bottom swing of the population rate once settled, $A = \max_t E(t) - \min_t E(t)$, which is exactly zero at a fixed point and positive the instant a limit cycle exists (no period-fitting or cycle-centre needed, the appeal of the direct measure). This needs no center-manifold reduction and no assumption that the local picture is 2D — and it is self-checking, since a hidden degeneracy would show up as beating or quasiperiodicity in the amplitude trace rather than the clean law below. Below $I_\text{ext}^\star$ the perturbation decays to zero; above it the amplitude grows continuously from zero, with $A^2 \propto (I_\text{ext} - I_\text{ext}^\star)$ (linear fit slope $3.2 \times 10^{-4}$, $R^2 = 0.997$). The square-root amplitude law and the absence of any jump or hysteresis mark a supercritical Hopf — the network slides smoothly from silence into a small-amplitude gamma cycle as drive crosses threshold, exactly the soft recruitment cliff of nb025.

Figure 5. Supercritical Hopf — amplitude onset by direct simulation

Why that shape is supercriticality. The square-root law is not a fit — it is forced by the Hopf normal form, which is what lets the measured curve stand in for the Lyapunov coefficient we declined to compute. Near the crossing the oscillation amplitude $r$ obeys a radial flow $\dot r = \lambda\,r + \ell_1\,r^3 + \dots$, with $\lambda \propto (I_\text{ext} - I_\text{ext}^\star)$ the growth rate of the now-unstable pair and $\ell_1$ the first Lyapunov coefficient. The limit cycle is the balance $\dot r = 0$, i.e. $r^2 = -\lambda/\ell_1$: a stabilising cubic ($\ell_1 < 0$) gives a real cycle of radius $r = \sqrt{-\lambda/\ell_1} \propto \sqrt{I_\text{ext} - I_\text{ext}^\star}$ — the linear instability pushes the amplitude out, the cubic reins it in, and they settle at a small finite radius. So ”$A^2$ linear in $(I_\text{ext} - I_\text{ext}^\star)$” is equivalent to $\ell_1 < 0$, which is the definition of supercritical; a subcritical Hopf ($\ell_1 > 0$) has a destabilising cubic, no small stable cycle, and an amplitude that appears with a jump and hysteresis — neither in the sweep. (It is the dynamical twin of a second-order phase transition: the gamma amplitude is the order parameter and the exponent $\tfrac12$ the mean-field critical exponent.)

The slope carries one more reading. Since $A \propto \sqrt{I_\text{ext} - I_\text{ext}^\star}$, the sensitivity $dA/dI_\text{ext} \propto (I_\text{ext} - I_\text{ext}^\star)^{-1/2}$ diverges at threshold: the amplitude rises with a vertical tangent, so the rhythm switches on smoothly (no jump) yet steeply — most of its amplitude appears in a thin sliver of drive just past $I_\text{ext}^\star$. Soft, but sharp: the dynamical content of nb025’s recruitment cliff.

The cycle the onset produces, just above threshold:

Figure 6. The 4D limit cycle near onset — E leads I

Discussion

The 4D Hopf (Figure 2) locates the recruitment threshold $I_\text{ext}^\star$ and a gamma frequency $f^\star \approx 24.5$ Hz from the synaptic timescales and coupling matrices — in the spiking band, reproducing its $\tau_\text{GABA}$-dependence (Figure 3) and nb025’s ”$W^{EI}$ moves the cliff” asymmetry (Figure 4), though flatter in absolute Hz than the spiking network. The 2D reduction shows why the full 4D state is needed: with $W^{EE} = 0$ it cannot destabilise the silent fixed point (Bendixson–Dulac, Figure 1), because the oscillation is driven by the synaptic phase delays the reduction throws away. The mean-field thus pins the mechanism of the $f_\gamma$ factor in nb041’s affine law $r_E = p \cdot f_\gamma$ (the per-cycle participation $p$ being the separate, spike-resolved result of nb046); the quantitative frequency calibration is approximate and left to the spiking network.

A note on the “use a 2D model” instinct, since it is the natural one: it is right about the analysis and wrong about the object. All three results here are local near the Hopf — the Jacobian crossing, the frequency, the amplitude scaling — and a Hopf is generically a 2D phenomenon. But the legitimate 2D-ness is the center-manifold normal form derived from the 4D field, not a Wilson-Cowan or FitzHugh-Nagumo system substituted for it. The appendix carries out exactly that reduction — Kuznetsov’s first Lyapunov coefficient computed from the 4D Jacobian and gain derivatives — and it agrees with the direct simulation on the sign ($\ell_1 < 0$), the frequency, and the amplitude slope to ≈ 4%. Working in 4D coordinates throughout keeps the local analysis honest and keeps the parameters mapped onto things the circuit actually has ($\tau_\text{GABA}$, $W^{EI}$) rather than onto a fictional recovery variable. “Reduce locally for the analysis, keep 4D as the object” is one position, not two.

Appendix — the normal form, derived from the 4D field

§3.2 settled super- vs subcriticality by simulation, with no theory required. This appendix does it the textbook way as well — the first Lyapunov coefficient — and checks the two agree. It also makes concrete what “go locally 2D near the Hopf” actually means, since that phrase is easy to misread.

The idea, in words. Near the Hopf the 4D system has four modes (eigen-directions of the Jacobian): one oscillating pair, barely stable, that is about to become the gamma rhythm, and two others that decay fast. “Fast” relative to the barely-stable gamma mode means that whatever state you start in, those two collapse almost immediately, and the long-run behaviour — the limit cycle and whether it switches on softly or with a jump — plays out on a 2D surface carrying the gamma mode. That surface is the center manifold; the equation for the dynamics on it is the normal form. So near onset the essential dynamics genuinely are 2D.

Two things to keep straight, because they are exactly where the “2D” word causes trouble:

• Those 2D coordinates are not the rates $(E, I)$. They are the amplitude and phase of the gamma mode — a direction that mixes all four variables.
• This is therefore not the Wilson-Cowan / FitzHugh-Nagumo rate model §2.1 rejected. That model eliminates the conductances (wrong — they carry the delay); the center manifold keeps the full 4D field and merely drops the directions that decay fast, so the delay — hence the frequency $\omega_0$ — stays in. The center manifold is the legitimate “go 2D,” derived from the 4D field rather than substituted for it.

What we need. Three things, all read off the 4D Jacobian $A$ at the Hopf:

• the frequency $\omega_0$ — the imaginary part of the critical (purely imaginary) eigenvalue, $\pm i\omega_0$;
• the critical eigenvectors $q$ and $p$: $q$ points along the gamma mode ($Aq = i\omega_0 q$), and $p$ is its dual partner ($A^{\!\top}p = -i\omega_0 p$, normalised so $\langle p,q\rangle = 1$), used to measure “how much gamma mode” a given state contains;
• how the field’s curvature — its quadratic and cubic terms — bends the gamma mode back on itself.

On the center manifold the dynamics become $\dot z = (i\omega_0 + \lambda)\,z + c_1\,z|z|^2$, where $z$ is the complex gamma amplitude, $\lambda \approx \lambda'(I_\text{ext}^\star)(I_\text{ext} - I_\text{ext}^\star)$ is the growth rate past threshold, and $c_1$ is the cubic term that saturates it. One number decides everything, $\ell_1 = \mathrm{Re}(c_1)/\omega_0$: negative, the cubic pushes back → supercritical (soft onset); positive, it pushes out → subcritical (hard onset). This is the same $r^2 = -\lambda/\mathrm{Re}(c_1)$ balance behind §3.2’s square-root law — $\ell_1$ is just its sign, computed rather than measured.

The curvature is almost all zero. The only nonlinear parts of the 4D field (3)–(4) are the two gain curves: the E-equation bends through $\Phi_E(I_\text{ext} - g_i^E)$, the I-equation through $\Phi_I(g_e^I)$, and the conductance equations are perfectly straight lines. So the quadratic form $B$ (the second derivatives) and cubic form $C$ (the third derivatives) each have a single non-zero slot — with state order $(E, I, g_e^I, g_i^E) = (x_1, x_2, x_3, x_4)$,

with $\Phi''$ (the gain’s curvature) and $\Phi'''$ (its rate of change of curvature) evaluated at each gain’s fixed-point drive — $I_\text{ext}^\star - g_i^{E\star}$ for E, $g_e^{I\star}$ for I.

The formula. Feeding these into Kuznetsov’s projection formula — the standard recipe for reading $\ell_1$ off a vector field of any dimension (Elements of Applied Bifurcation Theory, eq. 10.59):

It looks forbidding but it is pure bookkeeping: project the cubic curvature onto the gamma mode, then subtract the two ways the quadratic curvature feeds back through the off-mode directions (the matrix solves $A^{-1}$ and $(2i\omega_0 I - A)^{-1}$ route those quadratic “overtones” back onto the mode), and take the real part. Two linear solves and three dot products, evaluated once at the fixed point.

The result. Refining the Hopf to where the pair is exactly imaginary ($I_\text{ext}^\star = 1.288$, $\omega_0 = 0.152$ rad/ms, i.e. $f = 24.3$ Hz — matching the sweep’s 24.5 Hz) and evaluating:

Negative, so the analytic verdict is supercritical — exactly §3.2’s, now with the center-manifold-validity assumption made explicit rather than implicit (the “only one mode crosses” condition that licenses the whole reduction was certified back in Figure 2).

Does it match the simulation? Quantitatively. The normal form predicts more than a sign. The cycle radius is $r^2 = -\lambda/\mathrm{Re}(c_1)$, and at leading order the measured peak-to-peak E amplitude is $A = 4\,|q_E|\,r$ (with $q_E$ the E-component of the gamma mode $q$), so

With $\lambda'(I_\text{ext}^\star) = 0.085$ and the eigenvector, this predicts an onset slope of $3.36 \times 10^{-4}$ — against the $3.21 \times 10^{-4}$ measured by direct simulation (Figure 5): a ≈ 4% match. The small gap is leading-order truncation — the simulation runs at finite amplitude out to $I^\star + 0.5$, where the higher-order ($r^5$) terms the normal form drops start to matter.

The verdict. Sign, frequency, and amplitude slope all agree between the analytic normal form and the direct 4D simulation. So “go locally 2D near the Hopf” is real and legitimate here — and it is this center-manifold reduction of the 4D field (its eigenvectors $q, p$, its curvature forms $B, C$), not the rate model §2.1 ruled out.