Rosetta Stone of Neural Mass Models

Continue to Appendix (Terminology, Coordinate Transformations, Bifurcations, Kuramoto, and more) →

List of Abbreviations

AIT

Algorithmic information theory

AMPA

-amino-3-hydroxy-5-methyl-4-isoxazolepropionic acid (receptor)

BOLD

Blood-oxygen-level-dependent (signal)

CSD

Current source density

DBS

Deep-brain stimulation

DCM

Dynamic causal modelling

DDE

Delay differential equation

DHO

Damped harmonic oscillator

dMRI

Diffusion MRI

DMF

Dynamic mean-field (model, Wong–Wang/Deco)

E–I

Excitatory–inhibitory (population motif)

EC

Effective connectivity

EEG

Electroencephalography

EIF

Exponential integrate-and-fire (neuron)

ERP

Event-related potential

FC

Functional connectivity

FCD

Functional connectivity dynamics

fMRI

Functional magnetic resonance imaging

GABA

-aminobutyric acid (receptor)

HFO

High-frequency oscillation

HH

Hodgkin–Huxley (model)

HWHM

Half-width at half-maximum

ING

Interneuron gamma (mechanism)

JR

Jansen–Rit (neural mass model)

LaNMM

Laminar neural mass model

LFP

Local field potential

LIF

Leaky integrate-and-fire (neuron)

MEG

Magnetoencephalography

M/EEG

Magneto- and electroencephalography

MOM

Metastable oscillatory mode

MOU

Multivariate Ornstein–Uhlenbeck (process)

MPR

Montbri’o–Paz’o–Roxin (mean-field reduction)

MRI

Magnetic resonance imaging

MSF

Master stability function

NMDA

N-methyl-D-aspartate (receptor)

NMM

Neural mass model

NMM1

First-generation neural mass model (second-order synapses, static sigmoid)

NMM2

Next-generation neural mass model (exact mean-field, dynamic relation)

ODE

Ordinary differential equation

OU

Ornstein–Uhlenbeck (process)

PDE

Partial differential equation

PING

Pyramidal–interneuron gamma (mechanism)

PSD

Power spectral density

PSP

Postsynaptic potential

PV

Parvalbumin (interneuron type)

QIF

Quadratic integrate-and-fire (neuron)

RG

Renormalization group

RNG

Random number generator

SC

Structural connectivity

SDDE

Stochastic delay differential equation

SDE

Stochastic differential equation

SEEG

Stereo-electroencephalography

SL

Stuart–Landau (oscillator / Hopf normal form)

SOM

Somatostatin (interneuron type)

tACS

Transcranial alternating-current stimulation

tDCS

Transcranial direct-current stimulation

tES

Transcranial electrical stimulation

TMS

Transcranial magnetic stimulation

TVB

The Virtual Brain (simulator)

VIP

Vasoactive intestinal peptide (interneuron type)

WILCO

Wilson–Cowan (model)

tocdepth2

1 Introduction

“Understanding is the ability to see one thing in many ways.”
—R. P. Feynman

Neural oscillations span every measurable scale, from subthreshold membrane resonances to the rhythms recorded by MEG/EEG and fMRI. While these dynamics may reflect critical computational roles or pathological signatures, bridging such empirical observations with a heterogeneous theoretical landscape remains a practical challenge.

Mathematical neural mass models (NMMs) provide a coarse-grained, biophysically informed description of population dynamics that bridges microcircuit mechanisms with macroscopic signals and enables principled inference and prediction. In this context, a “neural mass" refers to a “lumped", zero-dimensional representation in space: each region or circuit is summarized by a small set of population-averaged state variables, often (but not always) expressed in a firing-rate format.

A mass is associated with a single point in space — or, more commonly in whole-brain modeling, with a parcel of a brain parcellation that stands in for an extended cortical or subcortical region. The whole brain is then represented as a discrete network of such point-like masses, with anatomical connectivity encoded in their pairwise couplings. When continuous spatial variation is required — for example, traveling waves, cortical patterning, or laminar gradients — this discrete description generalizes to neural field models [1]. The discrete parcel indices are replaced by a continuous spatial coordinate, and the finite system of ordinary differential equations (ODEs) becomes a partial differential equation (PDE; or integro-differential equation) in space and time.

Earlier syntheses have shaped the field from complementary angles: Ermentrout (1998) [2] emphasized spatiotemporal pattern formation; Ashwin et al. (2016) [3] surveyed oscillatory network dynamics and phase-reduction perspectives. In parallel, however, connectome-based whole-brain modelling has become a standard way to link anatomy, delays, and neuroimaging observables [4,5,6,7,8,9], and recent exact mean-field reductions have begun to connect spiking-network models directly to low-dimensional mass equations. As a result, many scientists now move between oscillator normal forms, rate-based neural masses with explicit synaptic filtering/transfer functions, and exact spiking-to-mass reductions — without a shared notation in which the moves compose.

Links between these formalisms exist across the literature, but they are scattered and often use different conventions for coupling, forcing, and biological interpretation. The goal of this review is to provide a practical Rosetta Stone: a scoped, step-by-step mapping that starts from canonical oscillator theory and shows how successive modelling choices—damping, forcing, nonlinear saturation, synaptic filtering, and transfer functions—lead to the neural mass models most commonly used in large-scale brain modelling. Our aim is to provide a coherent “ladder” that (i) makes assumptions explicit, (ii) standardizes how exogenous inputs and inter-areal coupling enter at the node level, and (iii) offers concrete recipes for assembling and interpreting connectome-coupled network models and their perturbational responses.

We emphasize the limited scope of this exercise. The mappings we provide are local bridges: explicit reductions valid near specific bifurcations (typically Hopf), under specific limits (weak coupling, strong attraction, exact mean-field), or under specific structural assumptions (current-based coupling, identical populations). This review is not an exhaustive survey, and we do not propose a single global theory unifying all neural mass formalisms across all parameter regimes. Local bridges between models already exist in the literature but are obscured by mismatched conventions; a step-by-step walk across model space—in a single shared notation, with the assumptions of each bridge made explicit—turns these scattered equivalences into practical translation rules.

To this end, the paper’s structure mirrors a spiral curriculum: each incremental step in model complexity is introduced first in its isolated, single-mass form and subsequently generalized to the networked or coupled scenario essential for whole-brain modeling.

We begin with the undamped harmonic oscillator—the archetype of pure phase dynamics—and sequentially introduce dissipation, external forcing, and nonlinearities. We emphasize how a basic push-pull motif underlies its oscillatory dynamics. This systematically leads us to the Stuart–Landau oscillator (SL), whose characteristic cubic nonlinearity delivers amplitude regulation with robust limit-cycle behavior. From this pivotal point, we pause to discuss some basic elements needed to establish firm connections with biology: synapses, which transform and delay signals arriving at populations, and transfer functions, which shape the response of accumulated synaptic perturbations into output firing rates. With this at hand, we jump to the Wilson-Cowan (WILCO) model, which provides a limit cycle linking back to the Hopf bifurcation in the SL model. Here we (typically) interpret abstract amplitude and phase coordinates as firing rates with biologically meaningful excitatory–inhibitory (E–I) population interactions, with the transfer function being the dynamical element. Second-order synaptic filters with static transfer functions naturally yield more nuanced models focusing on the dynamics of post-synaptic potentials, such as the Jansen–Rit and laminar neural-mass families (NMM1). This class of models can reproduce empirically observed alpha–gamma oscillatory interactions and respond realistically to physiological and pharmacological interventions. At each step a shared push–pull dynamical structure recurs, recognizing it across formalisms is what makes the cross-walk possible. Finally, we discuss next-generation neural mass models that can be derived exactly from spiking-network limits under mean-field assumptions—building on earlier exact macroscopic descriptions for theta-neuron networks and exemplified by the Montbrió–Pazó–Roxin reduction for quadratic integrate-and-fire (QIF) neuron populations—and we use this class (NMM2) to show how “dynamic” transfer functions arise from microscopic spiking mechanisms [10,11,12]. Beyond this lineage, the conductance-based, propagation-speed, decision-circuit, and dynamic causal modelling families listed in block (B) of Table Table 1.1 each sit at a one-step extension of the rate/PSP lineage; we credit them but do not re-derive them in detail.

1.1 Key Concepts

We collect here the definitions used throughout the paper.

Neurons, synapses, mesoscale. We focus throughout on the mesoscale, or “lumped,” description. Motivated by experimental measures like Local Field Potentials (LFPs) which aggregate activity over space, we replace the state of individual cells with variables representing the statistical mean of the population [13,14]. In this framework: a) input spikes become a continuous average firing rate input (spikes per second). b) synaptic dynamics act as a linear low-pass filter, smoothing average inputs into an average membrane potential perturbation. c) the sharp firing threshold of single neurons is smoothed by population heterogeneity and timing (variance in thresholds and states, event timing), resulting in a continuous, non-linear sigmoid transfer function that maps average membrane potential to an output average firing rate.

Computational modeling concepts

1. Neural mass vs. neural field. A neural mass model is a lumped (point) description of a neuronal population (e.g. a cortical column, parcel, or nucleus), in which the collective activity of many neurons is summarized by a finite set of population-averaged state variables. These variables may represent firing rates, mean membrane potentials, synaptic currents, or other mesoscopic quantities, and are typically governed by low-dimensional ordinary differential equations (ODEs), stochastic differential equation (SDEs) or, when propagation delays are retained, delay differential equations (DDEs). A neural field model is the spatially continuous counterpart, in which the same variables depend on position and interact through spatial kernels or PDE/integral operators. Both masses and fields can be rate-based, voltage-based, conductance-based, or derived as mean-field limits of spiking networks; thus, mass should be understood as lumped vs. spatially extended, not as a synonym for rate. We use a two-variable representation as the minimal autonomous system that supports oscillation. Real neural mass models are typically higher-dimensional; the pair is the recurring rotational sub-block, not a model in itself.

   In this review, the concept of “neural mass model” is used in two related but distinguishable ways in the literature. In the Wilson–Cowan / Freeman / Lopes da Silva / Jansen–Rit / Wendling / laminar neural mass model (LaNMM) lineage, and in its exact next-generation extension by Montbri’o, Paz’o & Roxin (MPR / NMM2), “neural mass” is the family name of a specific modelling tradition: lumped population descriptions written in terms of firing rates and/or post-synaptic potentials, with current-based synaptic interactions through linear filters, and a population-level transfer relation linking input current to output firing rate. Within this tradition the transfer relation may be a static sigmoid (Jansen–Rit, Wendling, LaNMM), a dynamically filtered sigmoid (Wilson–Cowan, where the firing rate evolves through ), or an exact dynamical relation derived from the spiking microscale (MPR / NMM2). In the dynamical-systems and mathematical-neuroscience traditions, by contrast, “neural mass” simply means lumped (as opposed to neural field, spatially extended), and the term covers any zero-dimensional population description regardless of whether its dynamical variables are firing rates, mean voltages, conductances, or order parameters of an exact spiking-network reduction, and regardless of whether synaptic interactions are current- or conductance-based. Under this broader definition, the Liley conductance-based mean-field [15], the Robinson corticothalamic propagation model [16,17], the Wong–Wang / Dynamic-Mean-Field reduction [18,19], and the David–Friston canonical microcircuit [6] are all neural mass models. We adopt this broader (lumped-vs.-field) definition throughout the review. The derivations and cross-walks that follow, however, concentrate on the current-based rate/PSP lineage and its exact next-generation reduction, because that is the sub-family within which the chain of dynamical equivalences — phase oscillator damped oscillator Stuart–Landau Wilson–Cowan NMM1 NMM2 — is tightest, and where a unified notation yields the most leverage. The conductance-based, propagation-speed, decision-circuit, and dynamic causal modelling (DCM) lineages are credited in block (B) of the roadmap Table and revisited in the synthesis of §§8, but they are not re-derived in detail here.

2. Oscillators and oscillations. We use “oscillator” and “oscillation” in two related but non-identical senses. In the dynamical-systems sense, an oscillator is an autonomous system that supports persistent rhythmic motion, most commonly through an attracting periodic orbit (a stable limit cycle), but also through neutrally stable centers in idealized conservative models. In the data-analysis sense, an oscillation refers to a signal component with a dominant timescale (often visible as a narrow-band spectral peak). These notions overlap but need not coincide: for example, a damped linear focus driven by noise can produce noise-sustained quasi-cycles (spectral peaks without a deterministic limit cycle), whereas chaotic systems are aperiodic yet may still show prominent spectral structure. In algorithmic information theory terms, we would say that a signal is an oscillation if it can be efficiently compressed by exploiting its approximate periodicity. This reflects the scientific observer’s perspective on modeling the phenomenon (see the Appendix U1 and the topology of the Stuart-Landau equation for a more in-depth discussion).

   Harmonic oscillators are prototypical oscillatory systems, linear and conservative models of a mass-spring system with an angular frequency and a sinusoid whose amplitude is fixed by initial conditions [20].

   Limit-cycle oscillators are nonlinear and dissipative; after transients, they settle onto a stable closed orbit.

3. A node is one neural mass; a network is a set of nodes connected by synaptic links, which can encode axonal delays and gains. Coupling just two nodes is enough for symmetry-breaking, phase locking, and collective bifurcations.

4. Push–pull motif as a local dynamical decomposition. Across the models reviewed here, oscillations typically arise from an interplay between (i) a variable that tends to drive activity away from a reference state (“push”) and (ii) a variable (or process) that provides a restoring or damping effect (“pull”). In a harmonic oscillator these are displacement and momentum (or two quadrature variables); in E–I population models they correspond to excitatory and inhibitory components of the loop; in complex normal forms they appear as the real and imaginary parts of . Importantly, in biophysically detailed conductance-based synapses, synaptic currents depend on the state through reversal potentials (“shunting”): the effective sign and strength of an interaction can therefore change with membrane potential. Throughout, the push–pull language should be read as an operating-point-dependent effective interaction (often after linearization and/or under current-based approximations), rather than as a universal, state-independent statement that “excitation always pushes” and “inhibition always pulls”.

5. Forcing is any external drive that breaks the autonomy of a node: input from another node in the network (which we call coupling below), electrical stimulation, pharmacological modulation, sensory pulses, or broadband synaptic noise from sources not explicitly modelled. When the unforced system has a strongly attracting limit cycle and forcing is weak, its leading-order effect is typically a modulation of phase captured by a phase response curve [21]. When attraction is weak and/or forcing is strong, amplitude deviations and qualitative changes of the attractor can no longer be neglected, and a phase-only description may fail.[22]

6. Coupling is forcing generated by other nodes: one neural mass serves as a time-dependent drive for another. Coupling may be instantaneous or delayed, linear or nonlinear, and weak or strong. In the weak-coupling/strong-attraction regime, coupling often reduces to an effective interaction between phases (yielding generalized Kuramoto-type descriptions); outside that regime, amplitude effects (and possibly multistability) become central. External stimulation can be viewed as a degenerate form of coupling in which the “partner” signal is prescribed by the experimenter or device.[23,24].

Mathematical tools

1. Fixed points and linear stability. Given a dynamical system with and at least continuously differentiable, a fixed point satisfies . Unless stated otherwise, “stability” in this review refers to linear stability: the local behavior is determined by the eigenvalues of the Jacobian (and, in delay systems, by the corresponding characteristic spectrum). A stable fixed point attracts nearby states, while an unstable one repels nearby states [25].

2. Hopf bifurcation (also Hopf-Andronov, HA). A two-dimensional system can begin or cease to oscillate in four ways ; we focus on HA bifurcation throughout, and refer the reader to [25,26,27] for the others. A Hopf bifurcation occurs when a complex-conjugate pair of eigenvalues of crosses the imaginary axis, exchanging the stability of a fixed point.

3. Differential Linear operator (): Every synapse in the model behaves like a linear filter: it receives an incoming firing-rate trace and converts it into a post-synaptic potential (PSP) . Written explicitly, this is a linear operation (convolution, see Appendix App. I) equation x(t)= K[r(t)] , equation or, equivalently, equation r(t)= L[x(t)] , equation where is the synapse’s impulse-response kernel and the inverse operator of , . Framing the dynamics in terms of makes the filter’s gain, decay time and delay visible in a handful of coefficients and shows how the population amplifies, attenuates or phase-shifts small perturbations [25]. In the simplest push–pull oscillator, there are two such operators, one excitatory and one inhibitory; extending the network adds one row per additional connections between populations.

4. Nonlinearity: Models beyond the harmonic oscillator include nonlinear terms. They reflect the transformation of synaptic inputs into firing rates by the neuronal population. This is necessarily a nonlinear function because firing rates are bounded above and below.

5. Phase reduction and phase–amplitude reductions. When a system possesses a stable limit cycle with strong attraction and is subject to weak forcing/coupling, its dynamics can often be reduced to a phase equation , leading to phase-oscillator networks. If amplitude excursions matter (e.g. near weak attraction, near bifurcations, or under stronger perturbations), one can extend this to phase–amplitude descriptions (e.g. phase–isostable reductions) that retain the leading amplitude modes in addition to phase. We will point out, as models become more biophysical, when phase-only reductions are appropriate and when full-state (amplitude-including) descriptions are required.

2 Linear Response and Phase-Only Limits

Linear and phase-only descriptions are a natural starting point in any physics-oriented review of neural mass models, for three concrete reasons. First, networks of linear (or linearised) damped oscillators admit closed-form expressions for stationary covariances, lagged covariances, and power spectra in terms of the Jacobian and the connectome, yielding analytic predictions for empirical observables—functional connectivity (FC), MEG/EEG power spectral densities—that are routinely fitted to data [30,31]. Second, every nonlinear neural mass model introduced later in this review (Stuart–Landau, Wilson–Cowan, NMM1, NMM2) reduces near a stable focus to a linear damped oscillator, so the present section supplies the universal local description that subsequent models inherit. Third, the linear network on a connectome is the natural setting for graph-spectral and connectome-harmonic decompositions, in which the dynamics diagonalise in the eigenbasis of the connectivity Laplacian and link directly to large-scale brain modes observed in MEG/EEG [32]. At its core, an oscillation requires (i) at least two dynamical variables—or one complex variable whose real and imaginary parts exchange “energy” or activity, pushing or pulling on the another—and (ii) a mechanism to rotate or cycle continuously through phase space. The simplest realisation of these requirements is the undamped, phase-only oscillator, in which a single phase variable advances uniformly while the amplitude remains fixed. This phase reduction is itself an explicit modelling choice—valid under weak forcing and strong attraction to a limit cycle—and serves throughout this review as the pedagogical starting point onto which damping, external forcing, and nonlinear amplitude dynamics are progressively grafted [33,34].

2.1 Undamped (Phase‐Only) Oscillator

The basal mechanism leading to oscillations in the activity of neuronal populations relies on the interplay between two subpopulations of neurons, one of them inhibitory and the other excitatory. We can represent the activity of these two populations by two variables and , respectively. In the simplest description, we can assume that decreases linearly, whereas increases also linearly. Assuming that the proportionality coefficient is the same in the two cases, we obtain a simple set of two coupled linear differential equations:

Defining

one obtains directly a rescaled model in polar form:

This is simply a phase (undamped) oscillator with natural angular frequency , which in our case determines the synaptic/membrane time constants of the populations. From (2.4)-(2.5) we immediately obtain , while . Thus this system is characterized by a constant (conserved) amplitude and a linearly advancing phase .

Equations (2.1)-(2.2) provide the fundamental push-pull motif underlying oscillations in all our models: whenever is positive, it drives negative, acting like an inhibitory force on ; when becomes negative, the sign flips and is driven upward, as if released from inhibition into excitation. In turn, a positive pushes upward—exciting —while a negative pulls downward, inhibiting . This continuous alternation of “push” and “pull” ensures that neither variable drifts off balance: instead, they chase each other around in a perfect circle of constant amplitude.

Introducing the complex variable

one finds

so that

Its solution is , from which and can be determined.

2.1.1 Effect of forcing

To illustrate the effect of a tonic drive (deterministic or stochastic), we augment the undamped oscillator (2.8) with a complex forcing term . First, consider a constant :

A purely real produces a DC offset , around which the trajectory circulates (see Appendix section §H.5 for more details). When varies in time, Eqs. (2.9)–(2.10) describe how the instantaneous bias modulates both amplitude and phase. To see this, we note that in the complex representation,

the term gently shifts the oscillator’s center and can transiently entrain or phase‐shift the trajectory without altering its underlying circular geometry (see Appendix §H.5 for further details on the effects of forcing).

2.1.2 Internal and external contributions to forcing; electric fields

Forcing appears in all models, internal (from other nodes, i.e., coupling) or external (unmodelled nodes, neuromodulation, electric fields). We allow for the possibility of forcing to be of a stochastic nature, indicating this with a hat notation (). Next, we divide forcing contributions between internal contributions from the network model in which the node lies (, which we include in models typically as coupling), and external to it (). We also separate the latter into external driving of physiological origin () and perturbations from an external electric field ,

where is some operator or function of the electric field. Thus, we will normally model as a sum of a deterministic component — typically from inputs from external populations and an electric field contribution from transcranial electrical stimulation (tES), transcranial magnetic stimulation (TMS) or deep brain stimulation (DBS), for example —, and an additive zero-mean stochastic contribution to forcing .

For example, in the case of weak electric fields at low frequencies (transcranial electrical stimulation or tES), the contribution from the electric field is of the form,

Here, the electric field effect is linear through a vectorial coupling constant () and captures ensuing membrane perturbations [35,36,37,38,39].

We will use the notation in Eq. (2.12) in the following sections, adding a node index if needed (i.e., for the forcing on node ).

2.1.3 Network of undamped harmonic oscillators

Here, we provide the general coupled network model and show how it connects to the Kuramoto model presented below. Following the additive forcing model for coupling, the equation for a network of oscillators is

To connect the coupled linear oscillator network in (2.15) with a phase‐only description, we write each state in polar form

Substituting into (2.15) and separating real and imaginary parts yields coupled equations for amplitudes and phases. With , the amplitude dynamics read

and the phase dynamics

which already has Kuramoto–Sakaguchi structure, but with time-dependent amplitudes . A genuine phase-only model is recovered as an explicit reduction choice: when the amplitudes relax to a fixed spatial profile (a slow manifold approached when the linear spectrum has a single marginal mode and all others decay), the amplitude equation (2.17) sets and (2.18) reduces to the constant-coupling form of Eq. (2.19) below.

A convenient way to make this connection more concrete is to add a weak uniform decay term, , and to interpret the resulting dynamics in terms of the spectrum of the linear operator [40,41]. For appropriate choices of and , all but one eigenmode decay, while a single collective mode remains marginal and rotates at a common angular frequency. In this collective oscillation regime, the amplitudes relax to a fixed spatial profile , so that (2.18) reduces asymptotically to a genuine phase model with constant effective couplings,

This construction shows how a linear complex network with global phase symmetry naturally generates a low‑dimensional manifold of phase dynamics that is well captured by Kuramoto–type equations [42,43,44]. At the same time, the Kuramoto model is usually taken in the opposite, more general direction: as a phenomenological phase description for weakly coupled limit cycles with symmetry, where the coupling matrix and coupling function are not constrained to arise from any particular underlying linear operator.

2.1.4 The Kuramoto model

The canonical model for studying collective phase dynamics is the Kuramoto model — “the hydrogen atom of synchronisation” due to its simplicity, analytical tractability, and extraordinary reach across disciplines [45,34,46,47].

Starting from the undamped limit cycle in Eq. (2.4), Kuramoto assumed that each oscillator interacts with all others only through the difference of their phases. For a population of units this yields

where is the natural frequency of oscillator and is a global coupling constant.

As additional motivation of this model, the sine term is the first (and usually dominant) Fourier component of any ‑periodic interaction function[48], and retaining only this term might capture most of the essential physics. In the classic all‑to‑all case with a unimodal, symmetric frequency distribution and pure sine coupling (no phase lag), increasing past a critical value produces a continuous (second‑order) transition to synchrony that is solvable exactly in the thermodynamic limit () [49]. More generally—depending on the frequency distribution, the presence of phase‑lag or higher harmonics in the coupling, and the network structure—the transition can be abrupt and hysteretic (first‑order) or continuous. No other phase‑only coupling law has proved as analytically transparent or universally useful.

For whole-brain simulations, we couple phases over a weighted network and include drive and noise:

Here is the connectivity matrix, a global coupling gain, the intrinsic angular frequencies, the delays, fixed phase lags, and the external forcing. In whole-brain applications, is most commonly estimated from diffusion magnetic resonance imaging (dMRI; tractography from white-matter streamlines) and from tract lengths and an effective conduction velocity .

The box below summarises the parameters of the phase-oscillator description and gives guidelines for using it to model neural populations.

2.1.5 Applications

Historically, the phase‑oscillator framework grew out of Winfree’s program in theoretical biology, which showed how large populations of weakly coupled limit‑cycle oscillators with heterogeneous natural frequencies could synchronize and be reduced to phase dynamics [50]. Kuramoto then analysed an idealized continuum limit with sinusoidal coupling, deriving a closed mean‑field description via a complex order parameter and predicting a continuous transition from desynchrony to partial synchrony [51]. His 1984 monograph consolidated the theory and notation used today, and later expositions placed the model as a paradigm for collective synchronization across physics, chemistry, and neuroscience [42,52,53].

Kuramoto‑type networks have become ubiquitous in computational neuroscience, appearing in everything from large‑scale cortex‑as‑a‑graph models that reproduce zero‑lag synchrony, realistic blood-oxygen-level-dependent (BOLD) fluctuations, and structured MEG amplitude envelopes[54,5,55,56] to theories of how rhythmic sensory inputs entrain cortical circuits at the level of individual columns[57]. This popularity stems from the model’s tight mapping to neural ingredients: each oscillator’s intrinsic frequency can represent the diverse frequency content of cortical rhythms, and the coupling constant collects the net gain of synaptic (as well as ephaptic or subcortical) interactions among neurons or brain regions. At the macroscopic scale, the emergence of phase coherence provides a principled analogue of the global signals captured by EEG/MEG[54,58].

As coupling increases, heterogeneous oscillators undergo a transition from desynchrony to partial synchrony[51,53], offering a mechanistic lens on how large‑scale neural oscillations (e.g., alpha) can arise gradually as effective interactions strengthen. Within connectome‑constrained implementations, this framework has clarified how zero‑lag synchrony can occur across distant cortical areas, how hubs facilitate inter‑modular coordination, and how network activity explores metastable, clustering regimes that wax and wane over time[59,60,61,62].

The same phase‑only formalism is well suited to external drive: because inputs enter as phase biases, periodic stimulation can lock phases over predictable bandwidths, accounting for stimulus‑locked entrainment and phase‑reset phenomena in EEG/MEG[57]. This logic extends beyond the cortex—for example, circadian networks in the suprachiasmatic nucleus are naturally captured as forced phase‑oscillator ensembles[63].

To bridge toward biological realism, neuroscience variants relax idealizations by introducing sparse/weighted structural connectomes, heterogeneous propagation delays, stochasticity, and higher‑harmonic couplings. These “realism knobs” generate chimeras (coexistence of synchronized and desynchronized populations), metastable clustering, and frequency‑dependent phase‑lag structure—while retaining a tractable phase core[60,61,64,65,66]. Importantly, when these extras are set to zero, the models reduce to the canonical Kuramoto equation, preserving interpretability and analytical leverage.

Finally, the phase‑only reduction is biologically insightful because many cortical E/I loops operate near a limit cycle, so the theory cleanly separates when a population fires (phase) from how strongly it fires (amplitude). This turns questions about perception, attention, communication‑through‑coherence[67], and intervention into questions about phase alignment and effective coupling—precisely the levers neuromodulators and stimulation can adjust. In practice, Kuramoto‑type reductions have informed control strategies in deep‑brain stimulation and desynchronization therapies[68,69,70,71], and physiologically motivated variants link coupling and phase shifts to transmitters and hemodynamics[72]; related approaches have begun to quantify coupling abnormalities in clinical cohorts[73].

2.2 Damped Harmonic Oscillator (DHO)

Introducing a real damping coefficient ( for decay, for growth) into Eq. (2.5) (unforced case) the equations of motion become

which, in Cartesian coordinates, reads

The fundamental “push–pull” interaction between and persists in the damped case: the term continues to inhibit (or disinhibit) exactly as in the undamped case, while drives . On top of this reciprocal coupling, each variable now experiences its own leakage (if ) or intrinsic amplification (if ) through the terms and . As a result, the two‐node loop still chases its 90° phase offset, but gradually spirals inward under leak or outward under growth. This interplay—orthogonal E/I feedback combined with uniform decay or gain—captures how real neural circuits blend balanced excitation and inhibition with membrane‐ or synaptic‐level dissipation (or recruitment), producing damped (or growing) oscillations rather than perfect, constant‐amplitude rotations.

Letting yields the compact complex form

2.2.1 Effect of Forcing

To introduce a constant (or very slowly varying) drive, one again writes

This moves the only fixed point of the system from to

From the linear stability analysis, trajectories spiral into (for ) or out of (for ) the off‐center equilibrium (see details in Appendix App. C). In polar coordinates (see Eqs. (2.22)–(2.23)), can exactly balance the leakage term . A constant real ensures that does not decay to zero, mirroring how a steady current injection holds a neuron at a depolarized potential. More generally, time‐dependent encodes transient pushes and pulls that can transiently boost amplitude, shift phase, or perturb the equilibrium away from its natural damped focus—preparing the oscillator for subsequent coupling effects in network models.

2.2.2 Coupling: Network of damped harmonic oscillators

Consider oscillators that, when uncoupled, each satisfy

where (real) is the common damping () or growth () rate, is the natural frequency, and is any external (complex‐valued) input. Introducing all‐to‐all diffusive coupling among these units yields the network equation

Here, each denotes the intrinsic frequency of node , is the global coupling strength, and the term represents a diffusive interaction that pulls each oscillator toward the network centroid . Consequently, this all‐to‐all linear network is precisely the linearized limit of the Stuart–Landau network (introduced in the next section).

For whole‑brain simulations in the damped regime, we place a linear resonator at each node and couple them through a weighted connectome with optional delays, additive drive, and noise:

Here is the complex state of node ; sets linear damping (decay time ) and is the intrinsic angular frequency; is an optional global coupling gain sometimes used to fit models; are (possibly directed) connectome weights; are propagation delays (set if ignored) and is the external forcing; Setting and taking recovers the all‑to‑all form in Eq. (2.29).

With linear coupling and real , the network coupling term in Eq. (2.30) unpacks as

i.e., the connection links push to push and pull to pull at equal strength, with no cross-talk between the push of one node and the pull of another. Allowing to be complex rotates the coupling by and mixes the two channels across nodes — the Sakaguchi–Kuramoto phase-lag structure that emerges in the polar reduction (2.18), with at a carrier frequency identifying the offset with a conduction delay [74]. The interpretation is mode-agnostic: may be a pair of firing rates, a post-synaptic potential and its time derivative, or any pair of variables that play the push–pull role of the local rhythm. Linear coupling in thus corresponds to current-based, frequency-independent inter-areal drive linearised around a stable operating point.

2.2.3 Applications

Equation (2.26) is the small-amplitude (linearised) limit of the canonical Stuart–Landau (SL) oscillator, whose full version is introduced below. In this regime, each node behaves as a damped complex oscillator whose real and imaginary parts jointly encode a local E/I loop: the real part plays the role of in-phase activity, the imaginary part the quadrature component, and the linear coefficient sets the decay time back to equilibrium while fixes the resonance frequency.

Linearity means that second-order statistics — instantaneous and lagged covariances, cross-spectra, and power spectral densities — admit closed-form expressions involving the Jacobian and the connectome (with optional delays), so functional connectivity and spectra follow without simulation. This is worked out explicitly for the SL whole-brain model by Ponce-Alvarez and Deco, who derive analytic formulas and show their accuracy near the stable focus [75].

Because of that tractability, the linear damped-oscillator network has become a standard baseline in whole-brain modeling pipelines: (i) as a linearized SL model with analytically tractable functional connectivity (FC) and power spectral density (PSD) for rapid parameter sweeps and state comparisons; (ii) as a multivariate Ornstein–Uhlenbeck (MOU) model on the connectome to fit time-shifted covariances and estimate effective connectivity; and (iii) as graph-spectral neural-field approximations in which connectome eigenmodes diagonalize the dynamics. Representative examples include the SL-linearisation for closed-form statistics and fast grid searches [75], MOU-based estimation of directed effective connectivity from fMRI [31], and analytic graph-neural-field or spectral-graph models that link Laplacian eigenmodes to MEG/EEG spectra and spatial patterns [76,77,78]. Together, these works show that much of the large-scale resting-state phenomenology of the brain can be captured by linear (or linearised) dynamics, a point underscored by systematic model-selection studies arguing that macroscopic resting-state activity is often best described by linear models [79].

Biologically, this phase-linear description retains the features most relevant at mesoscopic scales. The damping reflects local E/I balance and neuromodulatory tone; captures dominant local timescales; coupling rescales long-range synaptic, ephaptic or subcortical gain; and inputs represent afferents or stimulation. In this light, stimulation acts as designed forcing that reveals the network’s linear frequency response (and hence entrainment bandwidths). In contrast, pharmacology primarily shifts effective gain or damping and can therefore alter susceptibility and resonance.

Stimulation and pharmacology as perturbations of the linearised focus.

Linearising any of the network models in this review around a stable focus and writing the deviation recovers the complex damped network of Eq. (2.30), , with Jacobian , , , and Gaussian white noise of covariance . For the stationary covariance satisfies the Lyapunov equation , and the cross-spectral density factorises through the resolvent :  [30,31]. A pharmacological intervention enters as a parameter shift inherited from changes in local gain/damping () and effective coupling (). To first order the covariance and PSD respond as

A drug that nudges a region toward Hopf shrinks the real part of the corresponding eigenvalue of and sharpens the spectral peak as the inverse of the eigenvalue’s distance to the imaginary axis; an intervention that rescales effective coupling redistributes power across modes through the eigenstructure of . Stimulation reads the same equation with active. A periodic drive (e.g., transcranial alternating current stimulation (tACS) at carrier frequency ) admits the steady-state response : the same resolvent that sets the spectrum also sets the entrainment bandwidth and its spatial selectivity, with regional response amplified along eigenmodes of whose imaginary part lies near . Pharmacology and stimulation are two readouts of the same linearised susceptibility, identifiable from the same fits [80].

Finally, the linear network sits naturally at the base of a modeling “ladder”, where it has already provided evidence for the functional relevance of oscillations in the brain [81]. When nonlinear terms are reintroduced (full SL), one recovers amplitude dynamics, multistability, and turbulence-like regimes exploited to explain state-dependent changes (e.g., wake vs. sedation, psychedelic modulation) and nonequilibrium signatures. Yet even there, linear-response or linear-noise approximations remain invaluable for deriving analytic perturbation metrics (e.g. fluctuation–dissipation measures of nonequilibrium) and for mapping empirical changes onto interpretable parameter shifts [82,83,84,85,86].

3 Stuart–Landau Oscillator (SL)

In polar form writing yields

Here, the linear terms generate growth and rotation, while the amplitude‐dependent nonlinearities self‐regulate both amplitude and frequency, stabilizing the limit cycle at . More specifically, is the linear growth () or decay () rate, is the intrinsic oscillation frequency, governs nonlinear amplitude saturation, introduces nonlinear frequency modulation, coupling amplitude and phase. Equations ((3.1), (3.2)) represent the normal form of a Hopf bifurcation[90], which gives rise to the sustained oscillations characteristic of the SL (or sometimes called Hopf) model (see Figure Fig. 3.3 and Appendix App. C for details).

The amplitude equation drives toward

while the phase evolves at an amplitude‐dependent rate .

In Cartesian form with and , one obtains

which makes explicit once again the dominating, push-pull motif in the first-order terms. In complex form, in turn, we have:

Effect of Forcing To model constant or time-varying synaptic inputs, we add a complex forcing term ,

A constant displaces the limit cycle to a new fixed point , implicitly given by

generally yielding . Thus, provides a biologically plausible mechanism for input‐dependent modulation of both amplitude and frequency, although it may also alter the stability of the attractor (see Appendix §H.5).

Coupling Extending to a network of diffusively coupled SL oscillators, each node obeys

where is the intrinsic frequency of node , the coupling strength, and any external input. The diffusive term synchronizes the network by pulling each oscillator toward the ensemble mean, while nonlinear saturation ensures bounded amplitudes, giving rise to collective phenomena such as synchronization, amplitude death, and clustering.

For whole‑brain simulations with nonlinear nodes, we couple Stuart–Landau oscillators over a weighted connectome with optional delays, additive drive, and noise:

where is the complex state; and are the local SL parameters; is a global coupling gain; is the (possibly directed) connectivity matrix (set for all‑to‑all if no connectome is used); are propagation delays (set if ignored) and is the external forcing. We can rewrite these equations in polar coordinates to reconnect with the Kuramoto model [91]. For phase-oscillator networks with distance-dependent conduction delays—directly relevant to whole-brain SL models with heterogeneous tract lengths—Budzinski et al. [92] derive analytical predictions for the resulting spatiotemporal patterns by analysing the spectrum of a delay operator combined with the adjacency matrix, providing closed-form access to the locked, twisted, and wave-like states observed in delayed Kuramoto-type whole-brain simulations.

3.1 Applications

The Stuart–Landau equation arose as the weakly nonlinear amplitude equation for the laminar–turbulent transition in shear flows [93,94,95] and is now understood as the normal form of a Hopf bifurcation. Any sustained oscillation that emerges from an instability with weak nonlinearity admits a local SL description, which is what makes SL the natural meeting point for the rate, oscillator, and exact mean-field formalisms compared in this review.

In computational neuroscience, these same properties make SL a natural generative model for whole‑brain dynamics. At the node level, the bifurcation parameter controls the local working point—negative gives noise‑driven, damped fluctuations; yields large, susceptible excursions; and positive produces self‑sustained rhythms—while sets the intrinsic timescale. Network coupling on the structural connectome, together with finite conduction delays and stochastic drive, then determines how local oscillators form transient coalitions, travel as waves, or lock into metastable patterns. In resting‑state fMRI, SL networks tuned close to the edge of the Hopf bifurcation reproduce both static functional connectivity (FC) and its time variability (FCD), with the best fits occupying a narrow corridor of high metastability that effectively defines a dynamical cortical core [96]. Incorporating realistic delays in the same framework explains how fast local generators can coalesce into slow, spatially organized metastable oscillatory modes (MOMs) that appear and dissolve at reduced collective frequencies, linking anatomy to itinerant large‑scale patterns observed across modalities [97].

Electrophysiologically, SL models provide a bridge between structure and spectra. Allowing one or multiple resonant channels per region improves the correspondence to resting MEG: SL networks account for band‑limited envelope correlations, the location of spectral peaks, and the transient alignment of modes, with multi‑frequency instantiations offering the best cross‑band fits [98]. These same phase–amplitude dynamics, when embedded on the connectome with delays, rationalize why modest shifts in global coupling or delay reorganize band‑limited power and envelope FC without changing anatomy, providing a mechanistic map from structure to oscillatory phenomenology [97,99]. Multimodal comparisons (fMRI+MEG) using common structural priors further demonstrate that SL’s small, interpretable parameter set can jointly capture FC, FCD, and transient mode structure across measurements [7].

Casting SL network behavior in the language of turbulence sharpens the operational regime. When the model is fit to empirical amplitude turbulence and FC, the best‑performing working point is typically subcritical fluctuating—just below the Hopf bifurcation—where susceptibility and information‑encoding capacity are maximal; strength‑dependent perturbations in this regime reveal richer responsiveness than in supercritical limit cycling [85,100]. A particularly clean neuroscience instantiation of these ideas is provided by Freyer et al. [101], who showed that an SL-type model with multiplicative noise operating near a subcritical Hopf bifurcation reproduces the scale-invariant bistability of the alpha rhythm observed in resting-state MEG. In their formulation, noise drives stochastic switching between a low-amplitude focus and a coexisting limit cycle, generating the heavy-tailed amplitude distributions and -like spectral signatures characteristic of empirical alpha activity. This canonical example illustrates that noise near criticality is an organizing principle that supports physiologically realistic multistability in neural mass dynamics, and motivates the subcritical-fluctuating working point favored by SL whole-brain fits. Using closely related SL formulations, local‑coherence (“turbulence”) readouts distinguish wake, sleep, anesthesia, and pharmacologically altered states, positioning SL as a compact generative lens on state‑dependent information flow [102]. The same perturbative logic extends to psychedelics: after fitting LSD and placebo states, SL-based in silico stimulations predict enhanced sensitivity to strong inputs and characteristic turbulence signatures under LSD, consistent with empirical observations of expanded dynamical repertoires [103,104]. A complementary perspective frames such state changes through the geometry and complexity of the underlying state space, with plasticity-induced reshaping of attractors as the mechanistic substrate of psychedelic-driven dynamics [105].

Biologically, SL’s parameters expose the very levers experiments can turn. Interpreting the bifurcation parameter as net local E/I gain and neuromodulatory tone frames drugs as parameter shifts that move regions toward or away from oscillatory onset; interpreting external input as forcing turns stimulation into a probe of the network’s susceptibility, mode by mode. Personalizing SL fits thus yields subject-specific maps of working points and delays that generate testable predictions about which regions and frequencies will respond most — and how those responses reorganize whole-brain dynamics. Building on this, dynamic sensitivity analysis quantifies how localized stimulation perturbs subject-specific SL/WC fits and identifies the regions whose targeted modulation most efficiently drives transitions between brain states [106]. Finally, because SL admits a linear approximation near the stable focus, one can derive closed‑form spectra and covariances for rapid exploration and uncertainty quantification, then reintroduce nonlinearity as needed; this provides a transparent bridge between analytic tractability and the rich nonlinear phenomena that motivate the model in the first place [30].

4 Interlude: Synapses and Transfer Functionals

4.1 Synaptic Dynamics and operator formalism

Neurons have two main forms of communication: chemical and electrical. The first happens through neurotransmitters emitted by the presynaptic neuron after a spike [107]. The second is bound to the existence of gap-junctions between nearby cells [108]. It is estimated that chemical synapses drive the vast majority of neural communication in the mammalian brain, and for this reason, electrical coupling has often been considered as a secondary character in neural dynamics (see, however,  [109,110,111]). Thus, earlier work on NMMs focused on chemical coupling only.

For a single synaptic connection, the dynamics of neurotransmitter binding can be modelled as kinetic reactions [112,113,114]. Ultimately, both experimental and modeling results show that, typically, upon receiving a spike, the effect of neurotransmitter release to a postsynaptic neuron will first undergo an exponential increase of post-synaptic potential, followed by an exponential decay. Therefore, the postsynaptic-potential (mV) of a single, isolated neuron is usually modeled by the second order ordinary differential equation [114,115]

where and are the rise and decay times (ms), is the amplitude of the PSP (mV/kHz), and is the input term, modeling the arrival of presynaptic spikes.

While Eq. (4.1) corresponds to that of a single neuron , its linearity allows us to quantify the mean synaptic activity of neurons receiving and reacting identically to the same inputs as . Then, we obtain the same equation for the mean synaptic activity,

Interestingly, Eq. (4.2) shares the same structure as the equation of a damped, driven harmonic oscillator, which, as we will show in the following sections, allows it to display oscillations and exhibit diverse dynamical responses.

For a single pulse received at time zero (), the solution to (4.2) is

Usually, different types of neurotransmitters will result in different timescales. Table Table 4.2 contains some putative quantities for these parameters. Considering this variety of time scales, there are two limiting cases of Eq. (4.2) that are widely used in the literature. In some cases, one can simplify this equation by assuming that the rise and decay times are identical, . Then Eq. 4.2 reads

The solution of this equation for a single pulse at time zero () reads

which is sometimes referred to as alpha synapse.

Time constants for the principal receptor types — -amino-3-hydroxy-5-methyl-4-isoxazolepropionic acid (AMPA), N-methyl-D-aspartate (NMDA), and -aminobutyric acid (GABA, GABA) — are summarised in Table Table 4.2.

On the other hand, if one considers that the rise time is almost instantaneous, , then Eq. 4.2 reads

which is just an equation for exponential decay, thus a single pulse at gives the solution

Other types of neurotransmitter kinetics might result in more complex synaptic dynamics [112,113]. Of particular importance is the glutamate NMDA receptors, whose dynamics also depend on the voltage of the postsynaptic neuron [114,115].

The synapse as filter: operator form

How should the synapse equations be read? To shed some light on this, we recast Eq. (4.2) into operator form by rewriting it first, _r_d s(t) + (_r+_d) s +s = r(t), and defining the synapse linear operator

Then Eq. (4.2) can be written succinctly as

The synapse operator can be derived using the properties of linear time-invariant systems, and the impulse response of the neural mass has the form of Eq. (4.7), which is a model supported by experimental and theoretical studies [116,117,118,119] (see Appendix App. I for details).

The synapse operation can be cast as a causal filter. Since is a linear differential operator, using proper boundary conditions, we can define its inverse (another linear operator),

This shows that the synapse response (membrane voltage perturbation) is a filtered version of its firing rate input.

The same logic applies to Eq. 4.6, but this time the operator is first order,

Synapses are linear operators: filters that transduce incoming firing rate inputs into currents or PSPs. The effect of such a filter is to distort and delay the incoming signals. Their action can be characterized by the impulse response, i.e., how it responds to a sharp (delta function) input, . In the case of Wilson-Cowan, where the operator is first order, the response is simply a decaying exponential (see Fig. Fig. 4.4).

In more realistic models, such as Jansen-Rit (JR) [120], Wendling [121], or LaNMM [122,123], the linear operator is second order. The second-order operator is essential for realistically capturing synaptic conductance dynamics, as biological PSPs do not rise instantaneously. While a simple first-order exponential decay adequately describes passive membrane voltage relaxation, it neglects the finite time required for neurotransmitter binding, ion channel opening, and resulting conductance increase. Introducing separate rise () and decay () time constants (which may be the same numerically), the second-order operator accurately represents this two-phase process: a rapid initial conductance increase followed by slower conductance closure. This explicitly accounts for physiological delays and ensures PSP dynamics include a biologically realistic temporal scale, crucial for correctly modeling neuronal integration and timing-dependent neural computations.

Synaptic dynamics depend on the firing rate of the presynaptic population, which we now turn to.

4.2 Current-based versus conductance-based synapses

Throughout this review we adopt a current-based description of synaptic interaction, in which the synaptic input is treated as an additive perturbation to a postsynaptic state variable through the linear filter . This is a deliberate idealisation. Real chemical synapses are conductance-based: the synaptic current depends explicitly on the postsynaptic membrane potential through the reversal potential of the relevant ion channel,

where is the time-dependent synaptic conductance, generated by the same kinetic mechanisms that produce the second-order PSP shape derived above [112,113]. The dependence on is the source of shunting: as the membrane potential approaches , the driving force vanishes and the synapse delivers no current despite the elevated conductance, and if overshoots the effective sign of the interaction reverses. An “inhibitory push” can therefore become a “pull” or vanish entirely depending on the operating point, a property absent from the current-based filter. For the cross-walk developed in this review, the current-based form is retained because the rate-based and exact-mean-field neural mass models compared in §§§5–§7 are themselves formulated in current-based terms, and because near a stable operating point the conductance-based interaction linearises to a current-based effective coupling with state-independent effective gain. The conductance-based / shunting Liley-type mass listed in block (B) of the roadmap Table retains the full -dependence and is the natural starting point for readers requiring explicit shunting effects, including the voltage-dependence of NMDA-receptor dynamics noted at the end of §§4.1.

4.3 Transfer functional

Experimental results characterize the firing rate of a population of neurons as a function of its total input current or membrane potential. This approach extends the concept of -I curve, used to characterize the dynamics of single neurons [124,14], to a neural population. Accordingly, a population receiving an input current produces a firing rate given by

where is some operator (the transfer functional) that describes the process, including delays in response and non-linearity. For simplicity, this is sometimes simplified to a static and instantaneous nonlinear relation between the input and the output of a neural mass. The functional becomes then a function with

which receives the name of transfer function.

There are several possible choices for the shape of the transfer function (see Table Table 4.3). A common choice is the sigmoid function considered by Wilson and Cowan [125],

where is half of the maximum firing rate of each neuronal population, is the value of the current when the firing rate is , and determines the slope of the sigmoid at the central symmetry point . Wilson and Cowan argued that such a sigmoidal shape accounts for heterogeneity in either neural connectivity or excitability threshold. Freeman proposed another expression that he fitted to experimental data, obtaining a similar sigmoid shape [134,135].

Further theoretical and experimental research has validated that, indeed, neurons subject to an increasing input current produce a firing rate that, in many cases, can be captured by a nonlinear function. However, these functions do not always follow a sigmoid function. For instance, Amit and Brunel derived the transfer function that follows from considering a network of leaky integrate-and-fire neurons subject to Gaussian noise [127] Similarly, transfer functions for exponential integrate-and-fire, quadratic integrate-and-fire and others have been derived analytically [129,131,130,132]. We provide the transfer functions of some of these cases in Table Table 4.3.

Dynamical transfer functions, where the response of the population is not instantaneous, were addressed by Wilson and Cowan, who proposed a filtered, exponential decay dynamics towards the firing rate activity,

where is a time scale controlling the response time of the neuron. As mentioned at the beginning of this section, this exponential decay response has often been used interchangeably with the exponential decay dynamics presented in Eq. (4.6), although notice that the time constants and account for different biophysical quantities.

Using operator language with and inverse (a linear filter or convolution), this can be expressed as a non-linear filter

with the transfer functional is In this case, we talk of a transfer functional and not function: a nonlinear operator mapping total input currents to firing rate time series with delay governed by the time scale .

4.4 A simple self-coupled model and different limits

Let’s consider now a population of neurons of the same type that receive an external input and has recurrent connectivity given with a strength (see Figure Fig. 4.5). Then, we can assume that the total input received by the population is

where is an electrical admittance relating membrane potential perturbation and current. Notice that to derive the expression of the total input (4.18) we have assumed that the current flux generated by the recurrent coupling is proportional to the post-synaptic potential. This is an additional hypothesis that we will revisit later on.

Therefore, the dynamics of the average PSP of this population is given by

which we can rewrite in full operator form as the master neural mass model equation

with the transfer functional and the usual synaptic filter inverse operator characterized by some time scales — see Eqs. (4.17) and (4.10).

Equation 4.20 provides the minimal building blocks to derive a range of NMMs from simple principles. As we discuss next, taking or with first or second order synaptic operators, it covers the Wilson-Cowan (see Section §5), NMM1 (Section §6) and NMM2 (Section §7) formalisms. Finally, it also encompasses other logical variants — for example, a Wilson–Cowan style transfer function paired with second-order synapses. Although extensions of the Wilson–Cowan model do incorporate synaptic filtering or plasticity dynamics, to the best of our knowledge there exists no commonly adopted neural-mass model that preserves the original Wilson–Cowan architecture while explicitly modelling separate finite time constants for the membrane and the synaptic stages and, in particular, uses second-order synaptic operators — as proposed in our ‘master NMM’ formalism.

Limit of fast synapses (Wilson-Cowan model).

Taking the synaptic time constant to be much smaller than the membrane constant, , leads to the limit of fast synapses,

with the scaling constant associated with the operator. This is the stance taken by Wilson-Cowan (see Section §5), which we can rewrite in firing rate form as

with , or, equivalently in synapse potential (PSP) form as

In this limit, we can think of firing rate and synaptic PSPs as the same up to scale (there is no delay, only a rescaling). Thus, the Wilson-Cowan equations can be read as firing rate or synaptic equations, although the dynamics are due to membrane capacitance rather than synaptic currents.

Limit of fast membrane.

Conversely, if — the limit of fast membrane —, the firing rate ODE becomes a static equation,

This is the approach taken in the Jansen-Rit formalism (NMM1) when the synaptic equations are second-order, for example.

Delays and refractory period.

Finally, Equation (4.19) contains two time constants ( and ). However, there are two more time constants that play an important role in neural dynamics, and which we briefly mention here: time delays and refractory periods. Time delays introduce additional time scales that account for the latency of neurons or connecting circuit elements for reaching the threshold voltage and releasing the neurotransmitters. This effect might be due to different elements, for instance, the propagation of signals along the axon. This might be simply modeled by adding a delay in the synaptic dynamics

Notice that delay-differential equations have some special properties. For instance, their dynamics are infinite-dimensional, and thus, phenomena such as oscillatory activity can arise in instances where this was not previously considered.

On the other hand, refractory periods correspond to the inability of neurons to respond to external stimuli for a few milliseconds after an action potential. Whether refractory periods have a major role in the collective dynamics and function of neural populations is still a matter of debate. In the simple framework explained so far, there are two main ways to include a refractory period. On one hand, some researchers assume that transfer functions with a sigmoid function, i.e., a saturation for higher input, are already accounting for these refractory periods, hence the saturation. However, many researchers consider transfer functions derived from models, which do not saturate for large inputs. On the other hand, following Wilson and Cowan [136], the fraction of neurons susceptible to external inputs at any time is given by

where (ms) is the refractory period. Therefore, Eq. (4.16) should instead read

By considering that , the integral can be approximated as and the previous equation then reads

For the sake of completeness, we now reproduce Eq. (4.19) with both time delays and refractory period:

5 The Wilson–Cowan model (WILCO)

Following Eq. (4.21) (fast synapse limit, where firing rate input and PSP are related linearly without delay), or more concretely the form in Eq. (4.23), for a single population, we consider the model for a pair of coupled populations introduced by Wilson and Cowan [136],

with and representing synapse PSPs and are membrane time‑constants, the effective synaptic weights (positive for excitation, negative for inhibition), a tonic slowly varying drive that acts as the bifurcation (control) parameter, and the synapse gain term [48]. The sigmoid functions (see Eq. 4.15) are specific for each population.

In the form corresponding to Eq. (4.22), where represent firing rates, we have

Regardless of its synaptic or firing rate form, the original WILCO model describes dynamics stemming from transfer function delays. Despite its origins, the model is sometimes interpreted as describing first-order synapse dynamics. The form remains the same, but the interpretation changes (dynamics are induced by synaptic, not membrane delays).

A Taylor expansion shows that Eq. (5.2) contains the Hopf bifurcation machinery in the SL normal form plus additional bifurcations. Moreover, each coefficient is anchored to a biophysical mechanism (e.g. self‑excitation or synaptic delay ).

WILCO contains the Stuart–Landau oscillator as one of its parameter-tuned limits. When the cross-coupling is weak, the model supports a clean Hopf bubble (HB LC HB) with no fold or SNIC anywhere on the slice, and the center-manifold reduction at either end is the Stuart–Landau normal form. In this limit WILCO is SL. Once the cross-coupling product crosses a threshold set by the diagonal terms, a Bogdanov–Takens point enters the oscillatory window. The saddle-node, SNIC, and Hopf curves emerge from it as three branches of a single codimension-2 unfolding: the cusp brings bistability, the SNIC brings Class I excitability, and the system displays dynamics genuinely beyond SL’s reach. In this regime WILCO is more than SL. The two regimes are 1-parameter slices through the same equations, separated by exactly one codimension-2 point in parameter space. This is what gives WILCO its place in the Rosetta Stone: it specializes to SL — and to the PING/Jansen–Rit gamma rhythm — at one limit, extends to bistability and Class I excitability at the other, and the hinge between them is a single locatable point.

5.1 Effect of forcing

Shifting the tonic drive moves the mean input along the sigmoid with respect to the sigmoid threshold, i.e., the population operating point; because the slope changes with input, the effective linear gain and thus the distance to Hopf vary smoothly.[137]

Introducing an additional term to the excitatory equation in firing rate form reads

Near its bifurcation, the model becomes a forced SL oscillator[138]. may represent noise or a time-varying (e.g. sinusoidal, noisy, pulsed), or constant but transiently switched on and off forcing; in either case, it probes or entrains the dynamics without altering the underlying bifurcation point set by . Sinusoidal forcing carves out Arnol’d tongues of phase locking whose geometry parallels that of the analytical SL normal form, yet the effective tongue width is filtered by the sigmoid slope , a feature absent from purely polynomial descriptions.

5.2 Minimal ingredients for a Hopf bifurcation

Linearizing (5.2) around its fixed point and applying the Routh–Hurwitz criteria, a Hopf bifurcation occurs when with [139]. Physiologically, this translates into three rules:

1. Reciprocal E–I coupling () is indispensable; a single self‑exciting pool cannot realize a Hopf bifurcation.

2. Sigmoid slope. Either recurrent excitation or an external drive must position the operating point on the steep part of so that , recovering Ermentrout’s criterion for “sufficient self‑coupling or bias”.[140]

3. Cubic saturation. The curvature of supplies the cubic term that tames the growth of infinitesimal oscillations; piece‑wise linear gains suppress this mechanism.

Projecting the WILCO equations onto their two‑dimensional center manifold yields—after normal‑form transformation—the complex Stuart–Landau equation [141]

where the real coefficients and their imaginary counterparts are explicit functions of and the first two derivatives of ; their signs decide whether the Hopf is super‑ or sub‑critical. Conversely, the physiological variables follow , with the critical eigenvector. See Appendix App. G for more details.

The logistic transfer curve packages dendritic saturation, membrane noise, and threshold heterogeneity into one analytic function.[4]. Its first derivative shapes the linear gain, its second derivative injects the cubic non‑linearity that stabilizes the limit cycle—precisely the term in (5.4). In other words, the WILCO sigmoid is a built‑in normal‑form generator.

Two examples of WILCO bifurcation diagrams are provided in Fig. Fig. 5.6; they are 1-parameter slices through the same model equations, differing only in which external input parameter is varied ( or ) while the others remain fixed.

The model in panel (a) — -driven regime (see Table Table C.4 for parameters) — exhibits three different bifurcations as the external input to the excitatory population, , is increased. Initially, for low , the system has a single stable fixed point. As increases, the system undergoes a saddle-node (SN) bifurcation, where two unstable fixed points emerge. Further increasing leads to the collision and annihilation of a stable and an unstable fixed point, and, beyond this point, a limit cycle emerges. This transition corresponds to a saddle-node on an invariant circle (SNIC) bifurcation. Finally, the limit cycle vanishes in an HB bifurcation, and from that point the system only has a stable fixed point. Thus, near HB, WILCO reduces on its center manifold to a Stuart–Landau oscillator (Class II onset, finite frequency, vanishing amplitude); near SNIC it reduces to a theta-neuron / QIF normal form (Class I onset, finite amplitude, vanishing frequency); and near SN, to a saddle-node normal form (bistability, no oscillation). All three reductions are local pieces of the same Bogdanov–Takens (BT) unfolding, of which panel (a) is a 1-parameter slice.

The model in panel (b) — the driven regime-coupling (see Table Table C.4 for parameters) — exhibits a different bifurcation structure: a clean Hopf bubble bounded by HB at one end and HB at the other, with no fold and no SNIC anywhere on the slice. The cycle is born of small amplitude at HB, grows smoothly across the bubble, and dies of small amplitude at HB. The center-manifold reduction at both ends of the bubble is the Stuart–Landau normal form, and the cycle is locally Stuart–Landau across the whole window (Class II only, no Class I, no bistability).

By navigating this parameter space, WILCO effectively contains the Stuart–Landau oscillator as a parameter-tuned limit while adding the BT-organized cusp and SNIC structures as a functional extension. This allows the same model equations to switch between acting as a simple harmonic-like oscillator and a complex neural integrator.

5.3 Coupling

To build a whole-brain network of equal Wilson–Cowan E-I motifs, the excitatory populations are usually coupled via long-range glutamatergic projections [7]. Introducing both tonic/exogenous inputs () and dynamic drives () at the excitatory nodes, the coupled system reads (firing rate version)

Here denotes the (row-normalised) structural connectivity weight from node  to node , typically estimated from diffusion-MRI tractography.

Long–range excitation is integrated inside the sigmoid: distant axonal currents enter the same dendritic current pool through instantaneous synapses; the non-linearity then converts the total current into a bounded firing rate, preserving the physiological ceiling set by .

The network form in Eq. (5.5) has been the subject of extensive analysis in the dynamical-systems and computational-neuroscience literature, going well beyond the original two-population formulation of Wilson and Cowan. Foundational mathematical analyses of coupled WILCO units characterised the phase-plane structure, oscillation onset, and bifurcation skeleton of the canonical pair [142,143], and subsequent work has mapped the multistability, metastability, and pattern-forming regimes that arise when many WILCO nodes are coupled through delayed connectivity [144,145]. The treatment in the present section is therefore best understood as a unified-notation summary of an established line of work, repositioned to serve the cross-walk to the Stuart–Landau, NMM1, and NMM2 formalisms developed in subsequent sections.

Of particular relevance to the cross-walk this review develops, WILCO networks admit phase and phase–amplitude reductions back to generalised Kuramoto-type oscillator networks of the form discussed in §§2. Hlinka & Coombes [146] demonstrated that phase reductions of WILCO networks reproduce empirically observed structural–functional connectivity relationships in resting-state fMRI, providing an explicit bridge between the rate-based formulation and the phase-oscillator analyses of §§2. Daffertshofer & van Wijk [144] further analysed how the amplitude dynamics of WILCO modulate phase connectivity, clarifying the regime in which the phase-only description remains predictive. These reductions are the operational bridge that makes the WILCO,,Kuramoto leg of the cross-walk explicit.

Equations (5.5) implement additive coupling: each excitatory node receives the transduced firing rates of its peers weighted by . This choice reflects the physiology of long‐range excitatory fibers and avoids the homeostatic constraints of diffusive schemes, allowing the network to exhibit rich collective dynamics—from large‐scale synchrony and traveling waves to stimulus‐induced entrainment—while preserving the local Wilson–Cowan bifurcation structure.

Equation 5.5 is one among many possible combinations of populations, couplings, and forcings.

5.4 Applications

Historically, the Wilson–Cowan (WILCO) equations were introduced in two seminal papers in the early 1970s to capture the coarse‑grained dynamics of interacting excitatory and inhibitory neuronal populations [147,148]. By replacing spikes with smooth population activity and embedding saturating input–output nonlinearities, the framework established a tractable mean‑field language for multistability, oscillations, and pattern formation. Later syntheses clarified when the mean‑field approximation is valid, how fluctuations and delays can be incorporated, and how WC relates to other neural‑mass and field formalisms [149,150,151].

As a modeling workhorse, WILCO now spans scales—from local microcircuits to whole‑brain simulations coupled with connectomes derived from diffusion MRI—precisely because its parameters map cleanly to biology (excitatory/inhibitory gains, operating points, time constants, inputs) while retaining enough nonlinearity to express the canonical dynamical regimes. On the electrophysiology side, delay‑coupled WC networks fitted to human MEG reproduce band‑limited amplitude‑envelope correlations and phase‑locking across subjects when equipped with biologically plausible ingredients such as inhibitory synaptic plasticity and heterogeneous conduction delays; they naturally exhibit waxing–waning synchrony (metastability) and predict how changes in E/I balance or propagation speed reshape macroscopic spectra and coupling structure [152]. Related firing‑rate analyses show how excitatory and inhibitory feedback loops jointly regulate gamma rhythms—useful for interpreting resonance and entrainment bandwidths observed in M/EEG [153].

For fMRI, WILCO nodes embedded on the structural connectome have been used in multiscale pipelines that connect anatomy to resting‑state BOLD statistics and clinical phenotypes. Personalized WC-based models in major depressive disorder, for example, identified executive–limbic dysregulation consistent with empirical FC and symptom profiles [154].

Within The Virtual Brain ecosystem, WILCO is a standard regional model for synthesizing MEG/EEG/fMRI observables and probing how inter‑regional coupling, local gains, and delays generate subject‑specific variability in FC and its dynamics [155,156]. In task and method‑development contexts, WILCO local dynamics frequently serve as the generative “ground truth” for benchmarking estimators of frequency‑specific or task‑modulated connectivity from MEG/fMRI [157,158].

Because WILCO encodes excitation and inhibition explicitly, it offers a clean bridge to perturbation and inference. In Dynamic Causal Modeling (DCM) for fMRI, replacing the usual bilinear neuronal state equation with a WILCO‑type nonlinearity improves model evidence on multiple datasets while preserving physiological interpretability of effective connectivity and local transfer functions [159]. Clinically, introducing a non‑monotone (depolarization‑block) activation into a single WILCO microcircuit reproduces focal epileptiform activity and its spread, providing mechanistic handles for presurgical hypothesis testing [160]; conversely, disease‑specific applications at the whole‑brain scale use WILCO oscillators to explore how regional vulnerabilities and synaptic downscaling alter global connectivity and responsiveness [161]. Stochastic variants poised near critical points reproduce avalanche statistics and scaling of spontaneous activity, offering a testbed for hypotheses about critical brain dynamics and their departures under pathology [162,163,164].

Comparative studies situate WILCO among other whole‑brain models. Multi‑modal head‑to‑head work reports that WILCO and SL networks achieve broadly comparable fits to MEG and fMRI benchmarks once conduction delays and local E/I homeostasis are respected; WC often affords advantages on spatiotemporal measures (e.g., functional‑connectivity dynamics, the size distribution of transient oscillatory modes) thanks to its explicit rate saturation and E/I partition, whereas SL’s normal‑form compactness facilitates analytic reductions and turbulence‑style analyses [7]. Systematic benchmarking across cohorts further underscores that no single model dominates all metrics: for some summary statistics and parcellations, simpler linear baselines can rival or exceed nonlinear formalisms, and reliability/subject‑specificity depend strongly on the targets of fit [165]. In practice, WILCO is most compelling when questions hinge on mechanistic E/I balance, pharmacology, stimulation, seizure dynamics or nonequilibrium signatures; when phase‑only timing, graph‑spectral tractability or normal‑form universality are paramount, Kuramoto/SL (and linear surrogates) may be the better lens. In all cases, WILCO’s strengths and limitations are transparent. WILCO is a phenomenological rate model, originally formulated at the population level rather than derived from spiking dynamics: a compact E/I mean-field that is expressive enough to capture oscillations, multistability, and metastability while staying close to the biological levers experiments can manipulate.

6 NMM with second-order synapses (NMM1)

The equations in NMM1 are similar to the Wilson-Cowan, but introduce second-order derivatives. Here we present them without self-coupling for simplicity (unlike in WILCO, it is not needed for a stable limit cycle)—see Figure Fig. 6.7 (a),

which can be expressed in operator formalism as

with second order operator notation. For example, for the case of equal rising and decay times,

As in WILCO, the sigmoid function represents the integration carried out by the “soma" of the population. The argument of the sigmoid is therefore the total membrane perturbation caused by the PSPs from all synapses and other effects, such as electric fields. The s represent the coupling strength of the synapse — the synaptic gain.

Because the synaptic dynamics are governed by second-order operators (with rise and decay impulse response), the neuronal circuit can sustain oscillations even without explicit self-coupling. Specifically, the second-order filter introduces a frequency-dependent phase shift in the system response. According to the Barkhausen stability criterion, sustained oscillations occur if the total loop gain equals unity and the total loop phase shift reaches an integer multiple of 360. In a push-pull arrangement—excitatory coupled to inhibitory populations and back—this intrinsic phase shift, arising purely from synaptic kinetics (distinct rise and decay constants), can fulfill the Barkhausen criterion. Thus, unlike the Wilson-Cowan case, even in the absence of explicit self-feedback, second-order dynamics inherently provide the necessary conditions for sustained oscillations. We make this argument precise in the next paragraph.

Why second-order synapses sustain oscillations: the Barkhausen condition.

The defining dynamical feature of NMM1 is its capacity to sustain oscillations without explicit self-coupling, in contrast to WILCO. The mechanism is most cleanly understood through the classical Barkhausen criterion for closed-loop oscillation, which states that a feedback loop sustains a self-oscillation at angular frequency when (a) the loop gain equals unity and (b) the total loop phase shift is an integer multiple of . In the linear regime around a fixed point, the NMM1 push–pull motif (EIE) realises a closed loop whose loop transfer function reads

where is the sigmoid slope at the operating point and is the frequency response of the second-order synaptic operator (cf. Eq. (6.3). Two features distinguish NMM1 from WILCO here. First, the second-order filter contributes a frequency-dependent phase shift bounded by , so the two synapses together can supply the full required to close the loop, with the explicit minus sign of the inhibitory pathway providing the remaining . Second, the rise and decay constants control both the magnitude and the phase of , so synaptic kinetics alone determine where the loop-phase condition is met — and hence the resonance frequency. The loop-gain condition then selects the operating point at which oscillations onset through a Hopf bifurcation. The push–pull architecture is in this sense the minimal architecture that satisfies the Barkhausen criterion in NMM1: a single population with second-order synapses and no self-coupling cannot supply both the gain inversion and the loop phase. The full algebraic derivation of (6.4) and the explicit oscillation-onset condition are given in Appendix §I.3.

6.1 Forcing and coupling

Because of its realistic biological origins, accounting for the effects of coupling to other populations or of an electric field is straightforward: they produce additive voltage perturbation terms to the membrane potential, i.e., the argument of the corresponding sigmoid. So more generally, the equations with multiple synapses in a population reflect the additive combination of synaptic inputs.

The generalization of Equation 6.2 to multi-populations nodes and whole brain models consists of one for each synapse and neuron ,

(forcing can be represented through synaptic coupling). As before, the first equation transduces connectivity-weighted firing rate inputs into PSPs using the operator, the second sums all the PSPs and other perturbations affecting the neuron (e.g., an electric field in the equation), and the last one produces the firing rate output using the sigmoid.

The first equation links the input firing rate to its associated membrane perturbation (PSP). It can be read as the synapse equation for the input from neuron to neuron . The input may also reflect an input to the model from some external neuron, in which case for some function (constant, noise, etc.).

The second equation evaluates the membrane potential of neuron as the sum of the synaptic voltage perturbations plus an electrical field perturbation (if present) [166].

The last equation is a static transfer function: it evaluates the firing rate of a cell as a function of its total membrane perturbation.

The connectome in Equation 6.5 includes intra-parcel (defining, e.g., Jansen-Rit or LaNMM nodes) and inter-parcel connectivities in whole-brain models. The latter are usually derived from diffusion MRI or from Ising modeling of fMRI [167].

As a simple network example of inter-parcel coupling from excitatory to excitatory populations in the simple E-I motif model in Eq. 6.2 (with all parcels equal), we have

with as before, including noise or deterministic forcing.

As a further example, an electric field perturbation can be computed from the local electric field in the population. For example, in the case of weak electric fields at low frequencies (transcranial electrical stimulation), the perturbation is the dot product of the coupling constant and the electric field vector [35,36,38,37]. Adding this perturbation to the simple example leads to (see Equation 2.12)

As with the WILCO networks of §§5, the rate-based formulation here admits explicit phase reductions back to generalised Kuramoto-type oscillator networks. Forrester et al. [168] derive such a reduction for Jansen–Rit-type NMM1 networks and use it to show how local node dynamics shape the emergent functional connectivity patterns observed in whole-brain simulations, providing the explicit operational bridge between NMM1 and the phase-oscillator descriptions of §§2.

6.2 Applications

We review canonical models using the second-order synapse formalism (see Figure Fig. 6.7) and summarize their characteristic features and uses.

Simple push–pull motif (PING-like) model

The simplest second-order neural mass captures the push–pull loop between an excitatory and an inhibitory population, the canonical PING motif. Second-order synapses (distinct rise and decay) provide the phase lag needed to meet Barkhausen’s condition for sustained oscillations without explicit self-coupling. Minimal two-population models reproduce gamma-band rhythms and noise-sustained oscillations near Hopf; frequency depends mainly on inhibitory decay and EI gain and can be shifted by drive or kinetics. Applications include modeling high-frequency oscillations (HFOs) at seizure onset [173]; mechanistic context from spiking/mean-field work on PING/ING is reviewed in Buzsáki & Wang (2012)[174] and Tiesinga & Sejnowski (2009)[175], and synchronization analyses such as in Börgers & Kopell (2003)[176] and Whittington et al. (2000)[177].

The Jansen–Rit model

The Jansen–Rit (JR) model [120] comprises three populations (pyramidal, excitatory interneurons, inhibitory interneurons) interconnected with second-order synapses and a static transfer function. It generates alpha-band activity and realistic evoked responses; its regimes (fixed point, alpha, spike-like) are organized by Hopf and other bifurcations [178,179]. JR serves as the neuronal model in Dynamic Causal Modeling (DCM) for M/EEG steady-state and evoked responses [6,180,181], linking synaptic gains/time constants to observed spectra and ERPs.

The model dynamics can be written as

The bifurcation diagram of the Jansen-Rit model is shown in Fig. Fig. 6.8. Going from left to right, an SN bifurcation creates two unstable fixed points. As the external input is increased, the system exhibits a subcritical Hopf bifurcation (HB); beyond this point, the stability of the fixed point shifts and an unstable limit cycle (depicted in light blue) is created, and the system exhibits two stable fixed points after the HB bifurcation. This bistability persists as the system undergoes another Hopf bifurcation, this time a supercritical one (HB), where the upper fixed point loses stability and a stable limit cycle emerges. In this regime, the system can display either oscillatory or stationary behavior, depending on its initial conditions. Further increasing causes the two lower fixed points to collide in a saddle-node on invariant circle (SNIC) bifurcation. Beyond this point, the system exhibits two stable limit cycles, corresponding to oscillations with different frequencies. Eventually, a saddle-node of limit cycles (SNLC) bifurcation occurs, where an outer stable and an unstable limit cycle collide and annihilate, leaving a single remaining limit cycle. This final oscillatory state persists until a last supercritical Hopf bifurcation (HB), after which the system returns to a stationary regime.

The Wendling model and its extensions

Wendling’s CA1-inspired extension adds fast and slow inhibitory subpopulations (GABA, GABA) to the JR scaffold. By tuning inhibitory gains and kinetics, it reproduces background alpha, interictal spikes/spike–waves, and low-voltage fast activity; seizure onset emerges with impaired dendritic inhibition [182]. The model is widely used for interpreting stereo-electroencephalography (SEEG) / EEG patterns, exploring ictogenesis mechanisms, and assessing interventions [183], making it a standard computational tool in epilepsy. Recent extensions include chloride dynamics and laminar integration [184]. Whole-brain use: Wendling nodes have been used for resting‑state whole‑brain network modeling that matches empirical FC [185] and for patient‑specific whole‑brain simulations of interictal SEEG to aid clinical interpretation [186]. They have been embedded in realistic forward models to synthesize SEEG and personalize local epileptogenic dynamics [184]; the same group reports personalised whole-brain seizure-propagation models that integrate SEEG, MRI, and dMRI to simulate patient-specific responses to surgical, stimulation, and pharmacological interventions within a unified physiological framework [187].

As in the Jansen–Rit model, the core architecture is built around a pyramidal cell population interacting with both SS and SST neurons. The Wendling model extends this structure by introducing a parvalbumin-positive (PV) interneuron population, which both receives input from and projects back to the pyramidal cells. Incorporating this additional inhibitory loop results in the following system of equations:

The bifurcation diagram in Fig. Fig. 6.9 illustrates the steady-state membrane potential of the pyramidal population as a function of its external input . As increases, the system undergoes a supercritical Hopf bifurcation (HB) at approximately , giving rise to oscillations that persist until .

The laminar model (LaNMM)

The laminar neural mass model (LaNMM) [122,123] extends the NMM1 framework to capture the layered organisation of cortical microcircuits, addressing observables that single-population JR or Wendling models cannot reach. LaNMM couples two NMM1 generators with distinct anatomical and dynamical roles. A deep Jansen–Rit-like generator, situated in infragranular layers, produces alpha/theta-band rhythms through the standard three-population JR architecture. A superficial PING-like generator, situated in supragranular layers, produces gamma-band rhythms through the second-order push–pull motif of §§6. The two generators are coupled according to cortical laminar anatomy, with the deep generator modulating the superficial one through ascending projections that gate gamma amplitude on the slow phase, producing the empirically observed cross-frequency phase–amplitude coupling characteristic of laminar recordings. Because the model is anchored to specific cortical layers, the synthetic local field potential (LFP) and current source density (CSD) signals it generates inherit realistic depth-dependent polarity and phase relations through volume-conduction physics, allowing direct comparison with depth-electrode recordings rather than only with surface EEG/MEG.

Based on the connectivity scheme shown in Fig Fig. 6.7 (d), the model is governed by the following system of equations

The bifurcation diagram of the LaNMM is shown in Fig. Fig. 6.10. It depicts the steady-state membrane potential of the bottom pyramidal population () as a function of its external input (). The overall structure of the diagram closely resembles that of the Jansen–Rit (JR) model, as expected, since the LaNMM combines features of both the JR and PING models. As the external input increases, the system undergoes a saddle-node (SN) bifurcation where two unstable fixed points emerge. However, unlike the JR model, the LaNMM does not exhibit bistability. Oscillations then arise through a saddle-node on invariant circle (SNIC) bifurcation. Further increases in lead to a pair of fold-of-limit-cycle (FLC) bifurcations, where the oscillation amplitude decreases while the frequency increases (see [188] for details). This behavior continues until the system passes through a Hopf bifurcation (HB). Beyond this point, and in contrast to the JR model, the system undergoes an additional HB bifurcation, associated with the PING mechanism, giving rise to faster oscillatory activity. A systematic bifurcation analysis of the LaNMM identifies the parameter regimes that sustain coupled multifrequency oscillations when an external input is also introduced to the top pyramidal population, in addition to the previously considered input to the bottom population, and determines the conditions under which AD-like alterations disrupt these dynamics [188].

Applications span (i) reproducing laminar spectral peaks and CSD sinks/sources observed in primate and rodent cortex  [189,190,191,192]; (ii) explaining the cross-frequency coupling signatures observed in laminar recordings; (iii) capturing the alpha-slowing, gamma-loss, and disrupted theta–gamma coupling that characterise the oscillatory phenotype of Alzheimer’s disease [193], and the differential effects of psychedelic compounds in AD [194]; (iv) integration with stimulation physics for transcranial electrical stimulation (tES/tACS) studies through coupling to the local electric field as in Eq. (6.7); and (v) deployment as a regional model in connectome-scale simulations to study gamma-band coordination and cooperative/competitive network rhythms [195]. LaNMM has also been used as a generative model for predictive-coding interpretations of laminar dynamics [196].

In the disease-modeling space, modulation of fast-interneuron coupling within LaNMM reproduces the biphasic oscillatory progression of Alzheimer’s disease — early hyperexcitability with elevated alpha and gamma power, followed by oscillatory slowing and reduced spectral power consistent with empirical biomarkers [193] — and the corresponding modulation of layer-5 pyramidal excitability characteristic of 5-HT-receptor activation reproduces the spectral signatures of serotonergic psychedelics, suggesting a candidate therapeutic axis for the prodromal and early phases of Alzheimer’s disease [194].

7 Next-generation models (NMM2)

Classical NMMs are coarse‑grained descriptions: they compress the high‑dimensional, spike‑resolved dynamics of local circuits into a few mesoscale order parameters (typically, population firing rate and a filtered postsynaptic potential). Conceptually, this follows Kadanoff–Wilson coarse‑graining: integrate out fast, microscopic degrees of freedom, retain slow collective variables, and allow parameters (gains, time constants, noise) to be renormalized by the elimination of small scales [214,215]. Empirically, data‑driven coarse‑graining of population activity can approach non‑Gaussian fixed forms with static and dynamic scaling—evidence that mesoscale statistics can be approximately scale‑invariant near special operating points [216].

Early NMMs (e.g., Wilson–Cowan, Jansen–Rit) thus combine a biophysically grounded synaptic filter with a phenomenological static nonlinearity — the same split made explicit in the previous paragraph. Population‑density and mean‑field limits provide a more principled route from spiking to masses [217,218], and field‑theoretic expansions make explicit how fluctuations and finite‑size effects correct mean‑field behavior—precisely the kinds of corrections induced by coarse‑graining [219]. In large‑scale modeling, these ideas motivate dynamic mean‑field reductions of biophysical networks into mesoscale nodes with a few state variables [9].

The mathematical foundation for these exact reductions was laid by the analysis of theta-neuron networks via the Ott–Antonsen ansatz: Luke, Barreto & So [12] provided the complete classification of macroscopic behaviour for heterogeneous theta-neuron networks, and So, Luke & Barreto [11] extended this analysis to networks with non-trivial topologies, in both cases obtaining low-dimensional ODEs for collective phase coherence. Because the QIF and theta-neuron models are equivalent under the standard quadratic-to-theta change of variables , these reductions and the parallel reduction of QIF networks via the same Ott–Antonsen manifold capture the same macroscopic dynamics expressed in different coordinate systems. Montbri’o, Paz’o & Roxin [220] (henceforth MPR) reformulated the reduction directly in physiological coordinates — population firing rate and mean membrane potential — yielding the two macroscopic ODEs that anchor the present section. The QIF coordinate system is convenient because it expresses the macroscopic dynamics directly in terms of measurable observables (, ); the theta coordinate system is convenient for the underlying Ott–Antonsen analysis. NMM2 inherits results from both lines of work; see also subsequent reviews and extensions [221,222,223,212,213,224]. In the MPR framework, the steady‐state relationship between firing rate and mean voltage defines a sigmoid‐type input–output curve (i.e., a static transfer function). However, when the full two‐dimensional dynamical system is considered (with and evolving in time), the effective gain becomes state‐ and history‐dependent, rather than being a fixed static curve.

Positive self‑coupling () shifts fixed points to higher rates and can bring the system closer to resonant/oscillatory regimes, aligning with the intuition that denser local recurrence enhances responsiveness to weak, spatially uniform drives (e.g., uniform electric fields) [220,225]. Exact mean‑field extensions capture synaptic filtering, electrical coupling, and other biophysics [211,212,213], and have been leveraged to model working‑memory circuits with short‑term plasticity [226].

NMM2 can be read as a principled coarse‑grained synthesis: it replaces Freeman’s static nonlinearity by the MPR dynamic firing‑rate relation, and completes it with biophysical synaptic filters (e.g., second‑order ‑synapses) and exogenous field coupling. In renormalization group (RG), and are the relevant mesoscale variables; synaptic and coupling parameters are effective (renormalized) couplings that depend on the level of coarse‑graining and circuit state. This yields (i) a physics‑based “dynamic sigmoid” with inertia and state‑dependent gain, (ii) a transparent link between microparameters and mesoscale responsiveness, and (iii) a natural path to incorporate fluctuation corrections when needed (finite‑size, correlations) [219].

In their uniform, mean-field derivation for a population of quadratic integrate and fire (QIF) neurons, Montbri’o et al [10] start from the equations for the neuron membrane potential perturbation from baseline,

In this equation, the total input current in neuron is and includes a quenched noise constant component drawn from a Lorentzian (Cauchy) distribution, the input from other neurons per connection received (the mean synaptic activation) with uniform coupling , and a common input . The common input can represent both a common external input or the effect of an electric field, e.g.,

for weak electric fields. Here is an external uniform current, and is the dipole conductance term in the spherical harmonic expansion of the response of the neuron to an external, uniform electric field. This is a good approximation if the neuron is in its subthreshold, linear regime and can be computed using realistic compartment models of the (see, e.g., [227] and [228]).

The mean synaptic activation is given by

where is the arrival time of the th spike from the th neuron, and the synaptic activation function, e.g., . Note that we can write , where is the synapse coupling strength (charge delivered to the neuron per action potential at the synapse) of each synapse the cell receives from the network (there are of them in a fully connected architecture with neurons).

We assume here, for simplicity, that all neurons are equally oriented with respect to the electric field. If the electric field is constant, variations in orientation can be absorbed by the quenched noise term. The total input is thus homogeneous across the population (does not depend on the neuron).

Starting from these, Montbri’o et al derive an effective theory for a single population in the large limit (Eq 12 in [10]),

To this equation, we add the usual operator dynamics for the synapse activation,

(see Figure Fig. 6.7 (A)). Here and are the population mean membrane potential and firing rate, respectively. The new parameters and refer to the centre and half-width-at-half-maximum (HWHM) of the Lorentzian (Cauchy) distribution for the quenched noise input . Note that the Lorentzian distribution does not admit moments of order in the conventional sense (the integral defining the mean is divergent), so is a location parameter rather than an expectation; under the Cauchy principal value, formally equals the median. The analysis in [10] hinges on the assumptions of all-to-all uniform connectivity (with synaptic weight ) and common input .

These equations can be read as a transfer functional mapping input currents into an output firing rate, as in the master NMM Eq. 4.20,

The transfer functional is captured by Eq. 7.4. For a self-coupled inhibitory population (ING) with second-order synapses, the system of equations reads

where is the synaptic time constant, and where we explicitly introduce , the membrane constant of the population [213].

7.1 Forcing and coupling

Forcing and coupling are natural in this model, inheriting naturally from the single neuron biophysics in Equation 7.1 through the term . We converge here on the notation used in Equation 2.12, where represents the input external to the network.

E-I push-pull motif.

For two coupled populations (see Figure Fig. 7.11), we have the six equations, in dimensional (membrane-time) form

Here , ; are the membrane time constants of the excitatory and inhibitory populations, and are the second-order synaptic operators (with their own AMPA/GABA synaptic time constants). The membrane time enters as in the exact dimensional MPR reduction [10,213]: the rate term becomes , the quadratic curvature term , and the synaptic current is scaled by ; the external input and forcing are unscaled. The sign conventions follow the rest of the manuscript: denotes excitatory self-coupling within the E population, denotes inhibitory self-coupling within the I population (entering with an explicit minus sign in the equation), and , denote the magnitudes of the cross-couplings (EI positive, IE negative as a consequence of the explicit minus in the equation). The forcing term enters only the excitatory population, consistent with the convention used throughout this review. In operator form, following Equation 7.7, we express the complete set more compactly,

where is the dimensional MPR transfer functional—the pair above with membrane time —and, as usual, we drive only the excitatory population. The push-pull motif is present through the nonlinearity.

General equations for multiple coupled populations.

The general equations for multiple interacting populations become

where is the membrane time constant of population and the synaptic operators carry their own (synaptic) time constants, or

Equations for coupled E-I push-pull motif.

We can also write equations for multiple E-I motifs as in the previous sections, where each E-I motif is identical and we only connect excitatory to excitatory populations (with, e.g., and only two types of synapses),

The bifurcation diagram of Eq. (7.8) is shown in Fig. Fig. 7.12. It illustrates the steady-state firing rate as a function of the external input . As in the previous models, increasing the external input induces a transition from a stable fixed point to oscillatory activity via a supercritical Hopf bifurcation (HB). The resulting limit-cycle oscillations persist over a range of input values and are terminated at a second HB, with both bifurcation points indicated by the vertical dashed lines.

The bifurcation diagram of Eq. (7.9), which yields a PING model, is shown in Fig. Fig. 7.13. It depicts the steady-state firing rate as a function of the external input (I). From left to right, as the external input increases, the system first undergoes a pair of saddle-node (SN) bifurcations, indicating the creation and annihilation of fixed points. Following these, the system experiences a supercritical Hopf bifurcation (HB), giving rise to stable oscillatory dynamics. The resulting limit-cycle oscillations persist over a range of input values until the system passes through another HB bifurcation, at which the oscillations disappear and the system returns to a stable fixed point.

7.2 Applications

The NMM2 formalism has begun to power a range of concrete applications:

1. Emergence of brain rhythms: Several works have analyzed models based on the MPR theory to explore motifs allowing for the emergence of fast collective oscillations in one or two neural populations. The simplest instances include a single population of inhibitory neurons with synaptic delays[211,229,230,231], and excitatory-inhibitory population pairs [229,225,232,233,234].

2. Whole‑brain modeling: Next‑generation neural masses embedded on human connectomes offer an improved framework to analyze and reproduce brain dynamics captured by fMRI. Some works have started to explore this venue, focusing on capturing pathology and healthy states [235,236,237,238] and performing detailed mathematical analyses for the emergence of complex spatiotemporal behavior [239,240].

3. Non‑invasive stimulation theory: Because NMM2 replaces a static sigmoid with dynamic firing‑rate equations, it predicts state‑dependent field sensitivity. With second‑order synapses, it explains enhanced resonant responses to weak uniform AC electric fields (tACS‑like) and how self‑coupling and synaptic time constants tune this sensitivity [231]. Relatedly, NMM2 has been used to study population responses to brief exogenous pulses (TMS‑like) and ERD/ERS phenomena [225,241].

4. Cognition: Exact mean‑field models reproduce working‑memory operations and associated oscillatory signatures, providing a coarse‑grained yet mechanistic account [242,243].

5. Neuromodulation and DBS: Extensions that include adaptation/neuromodulatory variables have been used to explore mechanism‑level effects of deep brain stimulation and dopamine on network dynamics [244,245].

Additionally, several extensions to the MPR theory have extended its range of applicability by challenging some of the simplifying assumptions of the QIF formulation (7.1). Paired with synaptic dynamics, these extensions provide even further refined versions of NMM2:

1. Dynamic noise: Originally, the only source of microscopic disorder in the QIF formulation (7.1) was through the quenched variables . Nonetheless, the exact mean-field theory also applies if are considered to be dynamic Cauchy white noise [246]. Additionally, approximation theories exist for the case of Gaussian white noise [247,248].

2. Quenched noise distribution: The MPR theory requires to be Lorentzian-distributed to obtain a closed low-dimensional model for and . Recent work generalized the theory to include -Gaussian distributions, at the expense of increasing the number of independent variables in the resulting neural mass[249,250].

3. Gap-junctions: The MPR theory also applies for neurons communicating via electrical synapses mediated by gap-junctions and other diffusion-based coupling [44,251]. This allows for deriving NMM2 including electrical coupling, a major breakthrough considering that heuristic formulations cannot account for this microscopic effect.

4. Adaptation variables: In spite of its significance, the QIF model is a simplified model for neuron dynamics. Further biophysical mechanisms in the single unit dynamics, such as spike-frequency adaptation or ion channel dynamics, require the inclusion of additional dynamical variables for which the exact mean-field theory might not apply. To address this problem, some works propose including a mean-field adaptation as an approximation to networks with single unit adaptation[252,253]. A recent approach, instead, proposes including a suitable form of spike-frequency adaptation for which the MPR still holds exactly[254]. Beyond the QIF/theta line, Nicola & Campbell [255] derived an approximate mean-field reduction for heterogeneous networks of two-dimensional integrate-and-fire neurons.

5. Full low-dimensional theory: The MPR theory, which is strongly related to the Ott-Antonsen ansatz for the Kuramoto model[256], poses some challenging theoretical questions. The most relevant of them asks whether the low-dimensional manifold is attracting in the microscopic representation. A rigorous theoretical framework has been proposed recently, demonstrating that this is indeed the case provided the single units are not all identical ()[257].

6. Finite and asymmetric spikes: The original theory requires QIF neurons with a spike apex and reset at . This poses a limitation resulting neural mass model, especially when gap-junctions are included. Montbrió and Pazó [258] introduced a model accounting for spike asymmetry in order to properly capture these effects in networks with gap junctions. Recently, Cestnik [259] extended the exact low-dimensional reduction to a two-phase model capturing the dynamics for neurons with finite spikes ().

8 Summary and Outlook

We have treated neural mass models as points on a continuous ladder connecting microscopic spiking descriptions to macroscopic whole-brain dynamics. Griffiths, Bastiaens & Kaboodvand [260]organise the same space along an orthogonal axis — spatial scale (cells circuits networks), whereas the present review uses the derivation–representation axes summarised in Fig. Fig. 8.14. The two organisations are complementary.

At the lowest rung, the harmonic and Stuart–Landau (SL) oscillators provide the normal-form language for local rhythms and their phase–amplitude structure. Wilson–Cowan (WILCO) models make the underlying excitatory–inhibitory push–pull architecture explicit and introduce sigmoidal transfer functionals as coarse-grained summaries of neuronal input–output relations. Second-order synaptic neural mass models (NMM1) add realistic synaptic filters with distinct rise and decay constants, thereby separating synapse from soma dynamics and producing realistic postsynaptic potentials, delays and phase shifts. Finally, next-generation models (NMM2) derived exactly from quadratic integrate-and-fire (QIF) networks close the loop from spikes to masses by replacing static sigmoids with analytically derived dynamic transfer equations, showing how population firing rates and mean membrane potentials evolve jointly under recurrent and external drive.[220,221,222,213]

A main message is that all these formalisms share a simple core: a push–pull motif between two effective degrees of freedom, coupled through linear filters and nonlinear transfer functions, and embedded in a network via structured coupling and delays. In the linear limit, this core reduces to damped resonators or complex Ornstein–Uhlenbeck processes, for which covariances and spectra can be obtained in closed form and mapped onto empirical functional connectivity.[86,79] Close to Hopf, SL-like normal forms capture the onset of self-sustained oscillations, metastability, and turbulence-like dynamics under connectome coupling.[85,204,102] WILCO and NMM1 add biophysical levers (synaptic gains, time constants, E/I balance, laminar circuits) without losing this dynamical structure, and NMM2 shows that, at least for QIF networks, these macroscopic descriptions can be derived exactly rather than postulated.

Looking ahead, one natural direction is to extend exact mean-field reductions beyond QIF neurons. Current next-generation neural mass models exploit specific analytic properties of the QIF nonlinearity and Lorentzian excitability distributions.[220,221] A major open problem is to obtain similarly low-dimensional closures for more realistic single-cell models, such as exponential integrate-and-fire or multi-variable Hodgkin–Huxley-type neurons, possibly via systematic approximations (e.g. moment closures, population-density expansions, or renormalization-group inspired coarse graining).[217,219,247] Such derivations would clarify which aspects of spiking physiology (e.g. spike-frequency adaptation, active dendrites, NMDA currents) survive coarse graining and how they renormalize effective gains, time constants, and non-Gaussian noise at the neural-mass level.

A second axis concerns more realistic local circuits. Laminar neural mass models already capture distinct deep and superficial generator loops and their contribution to laminar LFP, CSD, and cross-frequency coupling.[8,122] Next steps include: (i) explicitly modeling multiple inhibitory subtypes (PV, somatostatin (SOM), vasoactive intestinal peptide (VIP)) with layer-specific projections; (ii) incorporating both chemical and electrical synapses in a unified mean-field description; and (iii) embedding short-term synaptic plasticity and intrinsic adaptation into next-generation masses in an exact or controlled approximate way [211,226,252]. These refinements would allow laminar models to speak more directly to cell-type and layer-resolved data, calcium imaging, and perturbation experiments, and to test hypotheses about how microcircuit motifs implement canonical computations (e.g., gain control, winner–take–all, and gating).

At the whole-brain scale, three modeling ingredients become increasingly important. First, heterogeneity across regions: empirical work shows that cortical areas differ systematically in local time scales, recurrent excitation, and laminar architecture, which in turn shape their role in large-scale hierarchies.[261,262] Future models should assign region-specific parameters (e.g., SL working points, WILCO gains, NMM1/NMM2 synaptic kinetics) informed by cytoarchitectonics, gene expression, and receptor densities, rather than using identical nodes everywhere. Second, directed and laminar-specific coupling: predictive-processing accounts suggest that feedforward and feedback pathways are implemented via distinct layers, with frequency-specific channels (gamma/beta) carrying prediction errors and predictions.[263,264] Embedding laminar neural masses in directed structural connectomes with frequency-dependent effective connectivity is a natural route to test such proposals against MEG/EEG and laminar recordings. Third, more refined parcellations and subcortical circuits: incorporating high-resolution cortical atlases (e.g. Glasser multimodal 360-area parcellation and its extensions)[265] and explicit thalamic, basal ganglia, and cerebellar masses[266,267,262] will be essential to capture cortico–subcortical loops that shape rhythms, state transitions, and neuromodulatory control.

Neuromodulatory systems provide a complementary, low‑dimensional control axis over the same large‑scale circuits. Gradients of receptor densities and gene expression co‑localize with the structural and timescale hierarchies described above, suggesting that neuromodulatory tone shapes regional gain, effective time constants, and long‑range coupling in a spatially specific way [268,269]. From a modeling standpoint, even when neuromodulatory nuclei are not explicitly represented, their effects can be approximated by treating external drives and gain parameters as proxies for neuromodulatory axes—for example, using global or projection‑specific changes in background input, SL working points, or WILCO gains to implement neuromodulation‑like shifts in E/I balance, integration time and noise statistics, in line with classical circuit‑level neuromodulation work [270]. Recent whole‑brain modeling incorporates neurotransmission‑weighted connectivity to account for a wide repertoire of task‑evoked states, showing that neuromodulator‑dependent changes in effective coupling can reproduce diverse cognitive configurations within a single structural scaffold [271]. These findings motivate configuring “external inputs” in large‑scale models using receptor and gene‑expression maps, or fitting low‑dimensional neuromodulatory control variables to match state transitions induced by pharmacological, arousal, or task manipulations.

We see two further connections, with statistical inference and with information-theoretic analyses of state. Linear and SL-based whole-brain models already support perturbation-based measures of nonequilibrium, susceptibility, and information routing across states (wake, sleep, anesthesia, psychedelics).[85,102,104] NMM1 and NMM2 provide richer, mechanistically grounded arenas in which to define and compute such quantities, including algorithmic-information or compression-based characterizations of oscillatory dynamics. Coupling these models to modern inference frameworks (variational data assimilation, Bayesian model comparison, active inference) can turn them into forward models for multimodal data (fMRI, M/EEG, SEEG, laminar probes) and for in silico perturbation experiments (TMS, DBS, tES).[9,159,80,272]

Acknowledgments

GR and FC have received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement No 855109, GALVANI) and from FET under the European Union’s Horizon 2020 research and innovation programme (grant agreement No 101017716, NEUROTWIN). PC has received financial support from the grant PID2024-155942NB-I00 funded by MCIN/AEI/10.13039/501100011033. Raul de Palma Aristides and Jordi Garcia-Ojalvo are supported by the European Commission under European Union’s Horizon 2020 research and innovation programme Grant Number 101017716 (NEUROTWIN). Jordi Garcia-Ojalvo was also financially supported by the European Research Council (ERC) under the Synergy grant 101167121 (CeLEARN), by the Spanish Ministry of Science and Innovation and FEDER under project PID2024-160263NB-I00), and by the ICREA Academia program.

Declaration of generative AI and AI-assisted technologies in the manuscript preparation process. During the preparation of this work the author(s) used ChatGPT in order to streamline the narrative. After using this tool/service, the author(s) reviewed and edited the content as needed and take(s) full responsibility for the content of the published article.

Continue to Appendix (Terminology, Coordinate Transformations, Bifurcations, Kuramoto, and more) →