Next-Generation Neural Mass Models

• NGNMMs are low-dimensional mean-field descriptions that exactly reduce spiking neuron networks to macroscopic equations representing firing rates, synchrony, and membrane potentials.
• They refine classical models by replacing static sigmoidal transfer functions with dynamic, coupled differential equations that capture resonant oscillations, bifurcations, and cross-frequency coupling.
• NGNMM frameworks extend to include adaptation, synaptic plasticity, and finite-size effects, enabling whole-brain embeddings and analysis of complex spatiotemporal dynamics.

Next-generation neural mass models (NGNMMs) are low-dimensional mean-field descriptions of large populations of spiking neurons that are derived exactly, rather than heuristically, from microscopic dynamics. In the formulations developed for $\theta$-neurons and quadratic integrate-and-fire (QIF) neurons, the macroscopic state is represented by dynamical variables such as the population firing rate, mean membrane potential, synaptic activity, and, in some formulations, synchrony. Because these models are exact macroscopic descriptions of underlying microscopic spiking neurodynamics in the thermodynamic limit, they occupy a distinct position relative to classical neural mass models based on sigmoidal transfer functions, and they have become a central framework for the analysis of oscillations, cross-frequency coupling, working memory, adaptation, finite-size fluctuations, and whole-brain dynamics (Coombes et al., 2016, Ceni et al., 2019, Clusella et al., 2022).

1. Exact mean-field reductions and canonical state variables

The defining feature of NGNMMs is exact reduction from heterogeneous spiking networks to a closed set of macroscopic equations. In the QIF setting, the reduction proceeds from a continuity equation for the probability density of membrane potentials, typically under Lorentzian-distributed excitabilities, and uses the Ott–Antonsen ansatz to collapse the infinite-dimensional system to a few collective variables. For inhibitory populations with exponential synapses, the macroscopic description contains three variables per population: instantaneous firing rate $r^{(l)}$, mean membrane potential $v^{(l)}$, and mean synaptic activity $s^{(l)}$ (Ceni et al., 2019). In the second-order synapse formulation, the exact model is written as

$$\tau_m \dot r = \frac{\Delta}{\pi \tau_m} + 2rv,\qquad \tau_m \dot v = \eta - (\pi r \tau_m)^2 + v^2 + \tau_m J s + I_E(t),$$

with synaptic filtering

$$\tau_s \dot s = z,\qquad \tau_s \dot z = r - 2z - s,$$

so that firing rate and mean membrane potential co-evolve with filtered synaptic input rather than being linked by a static nonlinearity (Clusella et al., 2022).

An alternative exact route uses networks of $\theta$-neurons. There the Ott–Antonsen reduction yields coupled equations for the Kuramoto order parameter $Z(t)$, which measures synchrony, and the synaptic conductance $g$. In that formulation the synaptic response is written as

$$Qg = \kappa f(Z),\qquad Q = \left(1+\frac{1}{\alpha}\frac{d}{dt}\right)^2,$$

and the firing rate is not postulated as a sigmoid but derived exactly as

$r^{(l)}$0

This introduces population synchrony as an explicit dynamical variable and allows event-related synchronisation and desynchronisation to be represented at the single-population level (Coombes et al., 2016).

The exactness of these reductions is conditional. The standard QIF reductions are exact in the thermodynamic limit for globally coupled populations with Lorentzian heterogeneity, and the resulting state variables are directly interpretable in terms of mesoscopic observables. In particular, the mean membrane potential is retained explicitly, which later work exploited as a proxy for local field potential and electroencephalographic signals during working-memory tasks (Taher et al., 2020).

2. Relation to classical neural mass models

Classical neural mass models are typically built around a static nonlinear transfer function, often sigmoid, that maps synaptic input to output firing rate. NGNMMs reject that approximation as a general principle: the firing rate is a dynamical variable obeying coupled nonlinear differential equations. The comparison between a heuristic model with second-order synapses (NMM1) and the exact QIF mean-field model (NMM2) shows that the two become mathematically equivalent only in the infinitely slow synapse limit, so NMM1 is an approximation to NMM2 only in that regime (Clusella et al., 2022).

For physiologically realistic parameter values, this equivalence breaks down. The exact model reproduces self-sustained oscillations in inhibitory interneuron networks via a supercritical Hopf bifurcation, whereas the heuristic model does not. The exact model also exhibits resonant amplification of firing rate in excitatory populations under oscillatory drive, with peak frequency and amplitude increasing with self-coupling gain and external excitatory input; the heuristic model instead remains low-pass and does not produce the same band-pass resonance (Clusella et al., 2022).

Aspect	Heuristic NMM1	Next-generation NMM2
Firing rate	Static function of input	Dynamical variable
Formal status	Approximation in the infinitely slow synapse limit	Exact mean-field reduction for QIF populations
Inhibitory oscillations	Does not reproduce self-sustained oscillations in the compared setting	Exhibits self-sustained oscillations via Hopf bifurcation
Response to periodic input	Low-pass filtered	Resonant and parameter-dependent

A related distinction concerns synchrony. In classical models, synchrony is typically not a state variable. In the $r^{(l)}$1-neuron NGNMM, by contrast, synchrony is represented by $r^{(l)}$2, and the firing rate becomes a consequence of synchrony rather than a postulated wave-to-pulse transform. This is one reason why NGNMMs can reproduce event-related synchronisation and desynchronisation at the population level (Coombes et al., 2016).

3. Oscillations, bifurcations, and cross-frequency coupling

A major application of NGNMMs is the mechanistic description of collective oscillations. In a single self-inhibitory population with exponential synapses, collective oscillations emerge through a super-critical Hopf bifurcation. Their existence depends crucially on three parameters: synaptic time scale, synaptic coupling, and excitability. Finite synaptic time is essential, strong self-inhibition is required, and reduced heterogeneity favors sustained oscillations; the frequency is controlled primarily by the synaptic time scale, with weaker dependence on coupling and strong sensitivity to disorder (Ceni et al., 2019).

For two inhibitory populations with different synaptic time scales, the same framework yields a much broader repertoire. In a master–slave configuration, the model supports damped oscillations toward a stable focus, periodic oscillations, quasi-periodic oscillations, and chaos. When the disparity in synaptic time scales is sufficiently pronounced, quasi-periodic and chaotic collective oscillations emerge; the reported macroscopic chaos is low-dimensional, with one positive Lyapunov exponent and a Kaplan–Yorke dimension slightly greater than $r^{(l)}$3 (Ceni et al., 2019).

Bidirectionally coupled inhibitory populations can autonomously generate theta–gamma cross-frequency coupling. Two forms were identified: phase–phase coupling, expressed as $r^{(l)}$4 phase locking such as $r^{(l)}$5, and phase–amplitude coupling, in which the amplitude of gamma oscillations is modulated by theta phase. External sinusoidal forcing at theta frequency enhances and widens the parameter region exhibiting theta–gamma coupling, increases resilience to excitability disorder, and can induce a transition from phase–phase to phase–amplitude coupling as disorder increases (Ceni et al., 2019).

A closely related set of results was obtained in forced PING and ING configurations. Near a Hopf bifurcation, a sinusoidal theta drive produces theta-nested gamma oscillations that always display phase–amplitude coupling. Two regimes were distinguished: a perfectly phase-locked periodic regime and an imperfectly locked regime with quasi-periodic or chaotic dynamics. The locked states were found to be more frequent in the ING set-up. In agreement with in vitro experiments, nested oscillations occur for forcing frequencies in the range $r^{(l)}$6 Hz, their amplitudes grow proportionally to forcing amplitude, and gamma activity is clearly modulated by theta phase. A documented limitation is that the gamma-power peak does not shift to higher frequencies as theta frequency increases unless noise amplitude or forcing amplitude are also increased (Segneri et al., 2020).

4. Adaptation, short-term plasticity, finite-size effects, and cognitive dynamics

NGNMMs have been extended to include spike-frequency adaptation (SFA), short-term synaptic plasticity (STP), and finite-size fluctuations without abandoning low dimensionality. In the mesoscopic SFA formulation for a single population,

$r^{(l)}$7

adaptation acts as a slow inhibitory feedback. In excitatory populations it favors population bursts; in inhibitory populations it hinders tonic population spiking. For two symmetrically coupled populations, SFA generates new collective regimes with cross-frequency coupling between the fast synaptic time scale and the slow adaptation time scale, including anti-phase slow-fast nested oscillations and symmetric or asymmetric bursting. A reduction of SFA was reported to increase theta frequency while decreasing gamma frequency, analogous to cholinergic modulation effects described for hippocampus and olfactory cortex (Ferrara et al., 2022).

STP extends the exact QIF reduction to synaptic facilitation and depression. In the synaptic working-memory model, the macroscopic variables remain the firing rate and mean membrane potential for each population, while the effective recurrent coupling becomes $r^{(l)}$8, with mesoscopic facilitation and depression variables obeying Tsodyks–Markram-type dynamics. Within this framework, memory access and loading are associated with stimulus-locked transient oscillations followed by steady-state activity in the $r^{(l)}$9–$v^{(l)}$0 band; facilitation and depression together support either synaptic reactivation or persistent activity; memory juggling and competition emerge already when only two items are loaded; and gamma power increases with the number of loaded items, whereas theta and beta power show non-monotonic behaviors. The model also yields an analytic expression for maximal memory capacity and identifies the mean membrane potential as a suitable proxy for memory load (Taher et al., 2020).

The STP setting also reveals a rich slow–fast geometry. In a next-generation neural mass model with plastic synapses, bursting can be organized by classical canards, torus canards, mixed-type canards, and spike-adding transitions. Near the singular limit, jump-on canards can block a continuous transition to bursting; in more biologically plausible regimes, the transition becomes continuous and bursts emerge through consecutive spike-adding transitions. Numerical evidence indicates that the same mechanisms organize bursting in the underlying QIF network, owing to the exactness of the mean-field limit (Taher et al., 2021).

A separate extension incorporates finite-size effects through shot noise. For large but finite populations, the discrete nature of spikes adds a stochastic term of order $v^{(l)}$1 to the macroscopic equations. The resulting shot noise is colored rather than white: its power spectrum can exhibit pronounced peaks at frequencies comparable to the mean firing rate, and resonance with collective damped oscillations can make its dynamical impact crucial even when its amplitude is weak in large massively connected networks (Klinshov et al., 2022).

5. Heterogeneity, Lorentzian tractability, and limits of universality

Lorentzian heterogeneity is central to most exact NGNMM reductions, but its special role is not merely technical convenience. For Lorentzian-distributed inputs, the dynamic mean-field reduction is two-dimensional, and there is an exclusive relationship between Lorentzian distributed inputs and Lorentzian distributed voltages. In stationary states, the voltage distribution is again Lorentzian, and its width coincides with the population firing rate. These properties explain why Lorentzian heterogeneity is analytically favorable and why the resulting firing-rate equations are especially tractable (Pietras et al., 2024).

The same study also establishes that exact mean-field reductions are not restricted to Lorentzian heterogeneity. For arbitrary unimodal heterogeneities, stationary firing rates and voltages can be written explicitly, and dynamic reductions can be constructed for distributions with finite complex poles, such as rational and $v^{(l)}$2-Gaussian approximations. The difference is dimensional: only the Lorentzian case reduces to a two-dimensional system; Gaussian or uniform heterogeneity requires higher-dimensional reductions that are more difficult to analyze (Pietras et al., 2024).

A common misconception is that Lorentzian-based NGNMMs are automatically representative of all heterogeneity classes. The comparison across heterogeneities shows that with chemical synapses the qualitative dynamics are usually similar and differences are mostly quantitative, such as shifts in bifurcation or synchronization boundaries. With gap junction coupling, however, Lorentzian heterogeneity induces nonuniversal behavior: it obscures a diversity-induced transition to synchrony that is present for Gaussian or uniform heterogeneity when the mean input is negative (Pietras et al., 2024). This marks an important boundary on the generality of the standard low-dimensional reduction.

6. Whole-brain embeddings, computation, and analytical outlook

NGNMMs are increasingly embedded in large-scale network models. One strand of work addresses the numerical problem directly through sparse Connectome Tensors that encode both connection strength and distributed delays, together with Local Linearization for semi-analytical integration of high-dimensional random differential equations. This framework was demonstrated on a single Zetterberg–Jansen–Rit cortical column and on networks of 1000 such columns, showing how distributed delays reshape spectra and phase portraits while maintaining computational tractability (González-Mitjans et al., 2020). Although not specific to one local model class, this computational infrastructure is directly relevant to large NGNMM networks with realistic connectivity and delay structure.

A more direct whole-brain application uses a 90-region anatomically coupled model in which each region is a next-generation PING module with excitatory and inhibitory QIF populations and exponentially decaying synapses. Stability analysis of homogeneous states and their transverse modes links unstable directions to the emergence of rich spatiotemporal patterns. Compared with classical neural mass models, the NGNMM formulation yields a broader dynamical repertoire in both homogeneous and heterogeneous regimes, including periodic traveling waves, high-dimensional spatiotemporal chaos with Kaplan–Yorke dimension up to $v^{(l)}$3 in a 90-node network, and gamma oscillations whose amplitude is modulated by slower rhythms through anatomical coupling alone (Delicado et al., 3 Dec 2025).

Analytical work has also begun to reformulate NGNMMs in the language of physics. After a state-dependent time rescaling, NGNMMs with instantaneous coupling become Hamiltonian in the homogeneous or strong-coupling limit, with an explicit energy-like quantity

$v^{(l)}$4

for the standard case. This energy-like function explains why phase-plane orbits are near-ellipsoidal, predicts spike amplitude during bursting, and offers a mechanistic bridge to energy landscape theory used in analyses of brain recordings (Andrean et al., 12 Sep 2025).

In the broader taxonomy of neural population models, NGNMMs now occupy the rung often labeled “NMM2”: models in which the transfer relation is dynamic rather than static. The proposed “ladder of detail,” extending from harmonic oscillators and Wilson–Cowan systems to second-order synaptic models and exact QIF mean-fields, situates NGNMMs as the point where micro-to-meso derivation, dynamic transfer functionals, and biophysical interpretability coincide (Castaldo et al., 4 Dec 2025). This suggests that the principal significance of NGNMMs lies not only in any single application, but in their role as a mathematically explicit bridge between microscopic spiking mechanisms, mesoscopic collective variables, and large-scale brain dynamics.