Neural Mass Models: Dynamics & Applications

• Neural Mass Models are aggregate models that condense large neural population activity into a few state variables governed by ordinary differential equations.
• They encompass phenomenological formulations, second-order synaptic filtering, and exact mean-field reductions that dynamically represent firing rates.
• NMMs facilitate analysis of brain rhythms, seizure mechanisms, and cross-frequency coupling through systematic coarse graining and advanced numerical simulations.

Searching arXiv for recent and foundational papers on neural mass models to ground the article. arXiv search query: "neural mass models overview next generation neural mass models exact mean-field" Neural mass models (NMMs) are aggregate neural models designed to reproduce the collective dynamics of neuronal populations. Rather than simulating every neuron, they compress population activity into a small set of variables—typically population-averaged firing rates, postsynaptic potentials, synaptic currents, or mean membrane voltages—governed by ordinary differential equations. In contemporary usage, the term spans phenomenological lumped-parameter models, firing-rate equations, and exact mean-field reductions of spiking networks, and it serves as a bridge between neuron-level mechanisms and macroscopic signals such as EEG, MEG, and fMRI (Clusella et al., 2022, Castaldo et al., 4 Dec 2025).

1. Conceptual scope and coarse graining

NMMs have been used since the 1970s to model the coarse grained activity of large populations of neurons and synapses, especially for the study of brain rhythms (Coombes et al., 2016). A central premise is that a neuronal population can be represented by a low-dimensional dynamical system whose state variables summarize the collective response of many cells. In the most common formulations, those variables are population firing rates and mean voltages or postsynaptic potentials; in more detailed formulations they also include synaptic currents, adaptation variables, or synchrony variables (Castaldo et al., 4 Dec 2025).

A recent unifying view organizes NMMs as a ladder of increasing detail: harmonic oscillator, linear damped oscillator, Stuart-Landau oscillator, Wilson-Cowan, second-order synapse models such as Jansen-Rit-type systems, and next-generation models derived from quadratic integrate-and-fire populations (Castaldo et al., 4 Dec 2025). Across this ladder, oscillations are treated as emerging from a push-pull interaction between two or more neural populations, most often excitatory and inhibitory populations. The ladder is not merely pedagogical; it formalizes how additional biological detail enters through damping, nonlinear amplitude regulation, synaptic filtering, and dynamic rather than static input-output relations (Castaldo et al., 4 Dec 2025).

The problem of what constitutes a “population” is itself nontrivial. Computational work on self-similar integrate-and-fire functionality proposed a coarse graining from neurons to ensemble-nodes, ensemble-spikes, and ensemble-time steps, and identified unique rescaling parameters that preserve integrate-and-fire-like behavior: ensemble-node size $N_N \approx 10$, ensemble-spike threshold $N_s \approx 5$, and ensemble-step $N_T \approx 4$ (Amgalan et al., 2020). This suggests that the scale at which a node is treated as a neural mass can itself be a model parameter rather than a fixed anatomical primitive.

2. Classical and phenomenological formulations

Classical NMMs are phenomenological. They use coarse-grained variables such as mean membrane potential and firing rate, and they typically postulate a nonlinear transfer function, often sigmoidal, to convert aggregate input into output firing activity (Coombes et al., 2016). Representative examples include Wilson-Cowan, Zetterberg, Jansen-Rit, Liley, and the canonical microcircuit model (Coombes et al., 2016, Cooray et al., 2022).

In Jansen-Rit-type formulations, cortical dynamics are generated by interactions among pyramidal cells, excitatory interneurons, and inhibitory interneurons. One common form is

$$\begin{aligned} V_1(t) &= h_e(t) * \left[c_1 S(V_3(t)) + P(t)\right],\ V_2(t) &= h_e(t) * c_3 S(V_3(t)),\ V_3(t) &= h_e(t) * c_2 S(V_1(t)) - h_i(t) * c_4 S(V_2(t)), \end{aligned}$$

with sigmoid

$S(V)=\frac{2e_0}{1+\exp(r(V_0-V))}$

and postsynaptic response kernels

$$h_e(t)=Aat e^{-at}, \qquad h_i(t)=Bbt e^{-bt}, \qquad t\ge 0$$

(Pei, 7 Mar 2025). This architecture is mesoscopic: it represents cortical columns or local circuits, while maintaining explicit excitatory and inhibitory subpopulations.

A common phenomenological extension introduces second-order synaptic filtering. In one formulation,

$\tau_s^2 \ddot{u}=C\gamma r(t)-2\tau_s\dot{u}-u,$

with firing rate given by a static nonlinear transfer function,

$r(t)=\Phi[I(t)], \qquad I(t)=\kappa u(t)+p+I_E(t),$

and corresponding state equations

$$\begin{aligned} \tau_s \dot{s} &= z,\ \tau_s \dot{z} &= \Phi[Ks(t)+p+I_E(t)]-2z-s. \end{aligned}$$

This is the heuristic “NMM1” form analyzed for second-order synapses (Clusella et al., 2022).

Analytical tractability has motivated simplified nonlinearities. In piecewise linear Wilson-Cowan systems, periodic orbits can be constructed by patching matrix exponential solutions across switching manifolds, and stability can be analyzed with Floquet theory. In the discontinuous Heaviside limit, the system becomes a continuous-time switching network and stability must be treated with saltation matrices, placing the model in the theory of nonsmooth dynamical systems (Coombes et al., 2018).

3. Exact mean-field reductions and next-generation models

Next-generation NMMs replace the postulated static transfer function by macroscopic variables derived exactly from microscopic spiking dynamics. For quadratic integrate-and-fire (QIF) populations with second-order synapses, the exact mean-field state is written in terms of firing rate $r$, mean membrane potential $N_s \approx 5$0, and synaptic variables $N_s \approx 5$1: $N_s \approx 5$2 In this “NMM2” form, firing rate is not a static function of the instantaneous input but a state variable with its own dynamics (Clusella et al., 2022).

The relation between heuristic and exact models is asymptotic rather than general. In the infinitely slow synapse limit, with $N_s \approx 5$3, adiabatic elimination of the fast $N_s \approx 5$4 variables yields an exact QIF transfer function,

$N_s \approx 5$5

so that NMM2 reduces to NMM1. Under this mapping, the heuristic model is the adiabatic approximation of the exact model for constant inputs and infinitely slow synapses (Clusella et al., 2022). Outside that limit, the dynamic firing-rate description is essential.

An alternative exact reduction proceeds from the $N_s \approx 5$6-neuron. There the population state is described by a complex order parameter $N_s \approx 5$7 encoding synchrony, with macroscopic equations of the form

$N_s \approx 5$8

where the firing rate $N_s \approx 5$9 is derived rather than postulated (Coombes et al., 2016). This formulation explicitly tracks synchrony or phase coherence, and thereby supports event-related synchronization and desynchronization at the single-population level (Coombes et al., 2016).

The treatment of heterogeneity is a further defining issue. Lorentzian-distributed excitabilities are analytically favorable because they permit exact low-dimensional closure, and the voltage distribution remains Lorentzian with half-width at half-maximum $N_T \approx 4$0 (Pietras et al., 2024). However, the same work shows that Lorentzian heterogeneity can induce nonuniversal behavior when gap junction coupling is present, obscuring a diversity-induced transition to synchrony that appears for Gaussian and uniform heterogeneity (Pietras et al., 2024). This is one of the clearest demonstrations that analytical convenience and dynamical generality need not coincide.

4. Dynamical repertoire: oscillations, resonance, itinerancy, and chaos

NMMs support a broad phase-space structure including stationary points, limit cycles, and transitions among different modes of activity (Cooray et al., 2022). In the canonical microcircuit model, second-order perturbation analysis together with adiabatic approximation yields amplitude dynamics as a gradient flow,

$N_T \approx 4$1

and, with stochastic forcing,

$N_T \approx 4$2

The resulting stationary density

$N_T \approx 4$3

provides analytic expressions for occupancy probabilities and Kramers-type dwell times, so oscillatory itinerancy, semi-stability, and multi-stability can be treated beyond local linearization around a fixed point (Cooray et al., 2022).

The difference between heuristic and exact models becomes most visible in their oscillatory behavior. For inhibitory QIF populations, the exact mean-field model predicts self-sustained oscillations through a Hopf bifurcation, whereas the heuristic static-transfer model exhibits only damped resonance. For excitatory populations, the exact model shows resonant transient oscillations and stimulation-induced amplification, while the heuristic model relaxes in a node-like, non-oscillatory fashion (Clusella et al., 2022). In the exact model, stimulation of a pyramidal population induces resonant oscillatory activity whose peak frequency and amplitude increase with the self-coupling gain and the external excitatory input (Clusella et al., 2022).

Cross-frequency structure has been a major target of next-generation models. In coupled inhibitory populations with different synaptic time scales, the macroscopic dynamics include damped oscillations, periodic oscillations, quasi-periodic oscillations, and chaos; under bidirectional coupling, the same framework yields theta-gamma phase-phase and phase-amplitude coupling (Ceni et al., 2019). In QIF-based PING and ING set-ups driven by sinusoidal theta forcing, theta-nested gamma oscillations emerge near a Hopf bifurcation. Two forms were identified: perfectly phase-locked periodic states and imperfectly locked quasi-periodic or chaotic states, with locked states more frequent in the ING set-up (Segneri et al., 2020).

Large-scale network embedding expands the repertoire further. A 90-region whole-brain model with next-generation local nodes and empirically derived anatomical connectivity supports homogeneous resting and oscillatory states, traveling waves, cross-frequency coupling, period-doubling cascades, and deterministic chaos. Stability of homogeneous states can be analyzed by dispersion relations, Master Stability Function methods, and Floquet theory, and the resulting dynamics can have Kaplan-Yorke dimensions up to 80 out of 90 possible (Delicado et al., 3 Dec 2025). Sparse random networks of exact firing-rate models show a related contrast between undirected inhibitory systems, which produce heterogeneous stationary winner-takes-all patterns, and directed networks, which produce high-frequency rhythmic states corresponding to high-dimensional chaos with extensive properties (Clusella, 15 May 2026).

Thermodynamic-limit reductions do not eliminate finite-size effects. For next-generation NMMs, finite population size introduces shot noise as a stochastic term in the macroscopic equations, and the shot-noise power spectrum can show pronounced peaks at frequencies comparable to the mean firing rate. Even when weak in massively connected networks, these fluctuations can strongly affect dynamics through resonance (Klinshov et al., 2022).

5. Inference, inversion, and numerical simulation

Because many NMM variables are latent, estimation of hidden states and parameters is a central methodological problem. One data-driven approach trains a bidirectional LSTM on simulated EEG from the Jansen-Rit model, with a loss combining data fitting, model consistency, and parameter smoothness. On simulated test data it yields correlations with $N_T \approx 4$4 of around 0.99, is robust to noise, and can outperform a nonlinear Kalman filter when initial conditions are inaccurate; in real intracranial EEG it reveals changes in connectivity strength parameters at seizure onset (Liu et al., 2023).

For structural inference, an adapted sequential Monte Carlo approximate Bayesian computation scheme has been developed for a 6N-dimensional stochastic multi-population Jansen-Rit model. The model includes binary parameters $N_T \approx 4$5 that specify the presence or absence of directed couplings, and inference uses marginal density, spectral density, magnitude squared coherence, and cross-correlation as summary statistics. On simulated data the method recovers cascade, partial, and fully connected topologies; on seizure EEG it identifies higher activation, higher input noise, and a clearer, denser network structure during seizure periods (Ditlevsen et al., 2023).

Large-scale simulation raises a separate numerical problem: NMMs with distributed delays and stochastic forcing can become very high-dimensional. One solution formulates the dynamics as nonlinear Random Differential Equations coupled through a sparse three-dimensional Connectome Tensor that encodes both connection strengths and delay distributions. Semi-analytical Local Linearization is then used to integrate the dynamics while preserving dynamical fidelity, enabling efficient simulation of a single Zetterberg-Jansen-Rit cortical column and large interconnected populations of such columns (González-Mitjans et al., 2020).

These developments extend a longer tradition in which analytical simplifications were used to make network dynamics tractable. In piecewise linear networked Wilson-Cowan models with circulant connectivity, stability of the synchronous state reduces to a low-dimensional Floquet problem parameterized by eigenvalues of the interaction matrix; with discontinuous nonlinearities, the analysis must instead propagate perturbations through switching manifolds by saltation matrices (Coombes et al., 2018). The contrast illustrates a persistent tradeoff in NMM methodology between analytic accessibility and biological richness.

6. Applications, translational roles, and recurrent controversies

NMMs are widely used to study pathological transitions and stimulation protocols. In a slow-fast seizure model with pyramidal neurons, slow SOM$N_T \approx 4$6 interneurons, and fast PV$N_T \approx 4$7 interneurons, the pre-ictal regime is organized by multiple time-scale dynamics and by a bifurcation structure involving folds and Hopf points. Within that framework, intermediate stimulation frequencies $N_T \approx 4$8 Hz can abort seizures if the timescale difference is pronounced, and the efficacy of stimulation depends strongly on which population is targeted (Ersöz et al., 2020).

NMMs have also been used as models of neural coding. In a feed-forward architecture of Jansen-Rit columns optimized by a genetic algorithm, different input classes are represented as oscillations of differing phase, typically around $N_T \approx 4$9–$$\begin{aligned} V_1(t) &= h_e(t) * \left[c_1 S(V_3(t)) + P(t)\right],\ V_2(t) &= h_e(t) * c_3 S(V_3(t)),\ V_3(t) &= h_e(t) * c_2 S(V_1(t)) - h_i(t) * c_4 S(V_2(t)), \end{aligned}$$0 Hz, and that phase code can be decoded by delivering impulses of excitation synchronized with the relevant phase-shifted oscillation. The synchronized input produces higher oscillatory power, so class identity is transformed from phase to power (Pei, 7 Mar 2025).

Recent work has added a physics-inspired perspective. After a suitable time rescaling, next-generation NMMs become Hamiltonian in the limit of a homogeneous population or strong coupling, and for small but finite heterogeneity they are near-Hamiltonian. The resulting energy-like quantity explains why phase-plane orbits are near-ellipsoidal and predicts spike amplitude during bursting, suggesting a possible link between neural population dynamics and energy landscape theory used in the analysis of brain recordings (Andrean et al., 12 Sep 2025).

Several controversies recur across the literature. One is whether population firing should be represented by a static transfer function or by its own dynamical equation; another is whether analytically convenient heterogeneity distributions distort qualitative predictions; a third concerns the scale at which coarse graining preserves neuronal functionality. A unifying proposal addresses these tensions by presenting NMMs as a navigable ladder of models, each considered in isolation, under forcing, and in a coupled network, with the explicit aim of turning NMM choice from a subjective act into a principled design decision (Castaldo et al., 4 Dec 2025). This suggests that “neural mass model” is best understood not as a single formalism but as a structured family of reductions, approximations, and exact mean-field limits whose suitability depends on the target rhythm, data modality, and intervention regime.