Rosetta Stone of Neural Mass Models

Rosetta Stone of Neural Mass Models Francesca Castaldo∗∗, Raul de Palma Aristides, Pau Clusella, Jordi Garcia-Ojalvo, Giulio Ruffini∗ July 2025 Abstract Brain dynamics dominate every level of neural organization—from single-neuron spiking to the macroscopic waves captured by fMRI, MEG, and EEG—yet the mathematical tools used to interrogate those dynamics remain scattered across a patchwork of traditions. Neural mass models (NMMs) (aggregate neural models) provide one of the most popular gateways into this landscape, but their sheer variety—spanning lumped parameter models, firing-rate equations, and multi-layer generators— demands a unifying framework that situates diverse architectures along a continuum of abstraction and biological detail. Here, we start from the idea that oscillations originate from a simple push-pull interaction between two or more neural populations. We build from the undamped harmonic oscillator and, guided by a simple push–pull motif between excitatory and inhibitory populations, climb a systematic ladder of detail. Each rung is presented first in isolation, next under forcing, and then within a coupled network, reflecting the progression from singlenode to whole-brain modeling. By transforming a repertoire of disparate formalisms into a navigable ladder, we hope to turn NMM choice from a subjective act into a principled design decision, helping both theorists and experimentalists translate between scales, modalities, and interventions. In doing so, we offer a Rosetta Stone for brain oscillation models—one that lets the field speak a common dynamical language while preserving the dialectical richness that fuels discovery. ∗[email protected], [email protected] (equal contribution) 1 CONTENTS 2 Contents 1 Introduction 4 1.1 Key Concepts .................................. 5 2 Harmonic Oscillator 9 2.1 Undamped (Phase-Only) Oscillator ...................... 10 2.2 Damped Harmonic Oscillator (DHO) ..................... 15 3 Stuart–Landau Oscillator (SL) 18 3.1 Effect of Forcing ................................ 19 3.2 Coupling ..................................... 19 3.3 Applications ................................... 21 4 Interlude: Synapses and Transfer Functionals 22 4.1 Synaptic Dynamics and operator formalism .................. 23 4.2 Transfer functional ............................... 25 4.3 A simple self-coupled model and different limits ............... 27 5 The Wilson–Cowan model (WILCO) 30 5.1 Effect of forcing ................................. 30 5.2 Minimal ingredients for a Hopf bifurcation .................. 31 5.3 Coupling ..................................... 31 5.4 Applications ................................... 33 6 NMM with second-order synapses (NMM1) 34 6.1 Forcing and coupling .............................. 35 6.2 Applications ................................... 37 7 Next-generation models (NMM2) 39 7.1 Forcing and coupling .............................. 42 7.2 Applications ................................... 44 8 Summary and Outlook 46 A Terminology and Scope of Aggregate Neuronal Models 68 B Coordinate Transformations 69 B.1 Undamped Oscillator: Polar ←→ Cartesian ←→ Complex ......... 69 B.2 Damped Oscillator: Polar ←→ Cartesian ←→ Complex ........... 70 B.3 WILCO and second order equations ...................... 71 C Linear Stability and Bifurcation Analysis 72 C.1 Linear stability analysis ............................ 72 C.2 Bifurcation diagrams .............................. 74 D Phase dynamics and Kuramoto model 81 D.1 Phase of a perturbed oscillator ......................... 81 D.2 From Winfree to Kuramoto .......................... 82 D.3 Phase reduction of a Stuart-Landau system .................. 84 CONTENTS 3 E System of Linear Coupled Complex Oscillators 86 E.1 Constant Forcing and the Affine Shift ..................... 87 E.2 From linear complex networks to Kuramoto phases ............. 87 F From Wilson-Cowan to the Damped Harmonic Oscillator 90 G From Wilson-Cowan to Stuart-Landau 93 H Stuart–Landau Oscillator: Parameters and Geometry 95 H.1 Geometry .................................... 95 H.2 Parameter Redundancy and Scaling in the Stuart–Landau Normal Form . 96 H.3 Alternative Form Emphasizing the Limit-Cycle Radius ........... 98 H.4 Stuart–Landau in Push/Pull Form ...................... 99 H.5 DC-Shifted Formulation of the Oscillator ................... 99 I Oscillations, Topology and Simplicity 103 I.1 Algorithmic-information definition .......................105 I.2 U1 and the topology of the Stuart-Landau equation .............107 J Linear Operators, Green–Laplace Tools, and E-I Oscillations: a Pedagogical View 110 J.1 Two complementary tools: Laplace and Green ................110 J.2 The synapse as a forced harmonic oscillator .................111 J.3 E–I motifs, Barkhausen conditions, and where the phase lag comes from . . 112 1 INTRODUCTION 4 1 Introduction “Understanding is the ability to see one thing in many ways.” —R. P. Feynman Modern neuroscience confronts us with an extraordinary diversity of dynamical phenomena, unfolding across an intricate hierarchy of spatial and temporal scales. Oscillations permeate neuronal systems, from subthreshold membrane resonances to the macroscopic rhythms observed in MEG/EEG or fMRI recordings, each reflecting underlying computational roles or pathological signatures. As experimental data accumulate, the theoretical neuroscientist faces a critical challenge: bridging these empirical observations with an equally vast and heterogeneous theoretical landscape. Mathematical neural mass models, also referred to as Firing Rate equations, offer a coarse-grained, biophysically informed description of population dynamics that bridges microcircuit mechanisms with macroscopic signals and enables principled inference and prediction. They range from minimalist phase oscillators, through amplitude-modulated systems, to elaborate firing-rate and lumped-parameter descriptions. Each formulation carries distinct assumptions, parameters, and interpretative frameworks, complicating efforts to unify insights or rigorously justify model selection. The field currently lacks a principled theoretical bridge—akin to a Rosetta Stone — that not only translates smoothly between these mathematical dialects but also guides objective model choice. Such a bridge should expose the biophysical correspondences among formulations, standardize how exogenous inputs (“forcing”) and inter-areal coupling are encoded at the node level, and provide clear recipes for assembling network models that can be interrogated by perturbations, including sensory drive and brain stimulation. Our ambition is to show how a simple linear oscillator—augmented systematically by damping, forcing, and nonlinearity—leads naturally to Wilson–Cowan firing-rate dynamics and then to layered neural-mass formulations, making explicit the assumptions introduced at each step. The present paper aims precisely to construct and elucidate this fundamental concordance, transforming a fragmented theoretical landscape into a coherent and navigable ladder. 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 1 INTRODUCTION 5 respond realistically to physiological and pharmacological interventions. Crucially, at each step we uncover a shared push–pull dynamical structure —the unifying mechanism across these different models. Finally, we discuss the first neural mass model rigorously derived from first principles (the quadratic integrate-and-fire neuron model, QIF), termed NMM2. In NMM2, the transfer function and synaptic couplings are both dynamical. Finally, we highlight the original push-pull core motif underlying the oscillatory behavior in all these models. The conceptual framework we develop yields important theoretical and practical benefits. By distilling the dynamical core of neural mass models to a minimal yet powerful set of parameters, we enable systematic transformations from one formalism into another. Thus, seemingly distinct modeling approaches become recognizable as points along a coherent ladder of abstraction and biological detail. Moreover, this unified view offers a principled rationale for model selection. Lastly, our approach emphasizes clarity and pedagogy: we provide explicit, step-by-step derivations throughout, ensuring the mathematical logic is transparent even to those new to the field. Our ambition is for the theoretical mapping presented here to equip both experimentalists and theorists with an intuitive yet rigorous language, enabling a fluent navigation of the modeling landscape. 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 wholebrain modeling. By following this progression, readers will gain insight not only into how these models interrelate mathematically, but also into why these interrelations are critically relevant for experimental design, clinical interventions, and the fundamental interpretation of neural data. 1.1 Key Concepts For clarity and ease of reference, we present here a concise set of key definitions and foundational concepts that recur throughout the paper. These will facilitate the subsequent mathematical development and make explicit the modeling assumptions connecting elementary oscillatory systems to neural mass formalisms. 1. Computational modeling concepts (a) A neural mass collapses the electrical activity of thousands of similar neurons into a handful of population-averaged state variables governed by lowdimensional ordinary differential equations, with interactions between populations mediated by effective synaptic couplings. More generally, we treat a neural mass as any variable pair (x, y) whose reciprocal interaction generates a collective mode: xpushes the state away from equilibrium, while ysupplies a restoring (or damping) pull. (b) We loosely call an oscillator a dynamical system whose state is approximately periodic in time. Similarly, an oscillation is a signal that approximately repeats. More formally, a dynamical variable is said to oscillate when it exhibits sustained, approximately periodic departures around a reference value such that the system returns to a similar state after a characteristic interval T (its period), or, equiv- 1 INTRODUCTION 6 alently, at a dominant frequency f = 1/T. The repetition can be exact (strictly periodic) or approximate - quasi-periodic, weakly modulated, chaotic (i.e., the Rossler chaotic attractor oscillates irregularly close to a given frequency), or noise-jittered. 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 I.2 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.1 Limit-cycle oscillators are nonlinear and dissipative; after transients, they settle onto a stable closed orbit. (c) 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. (d) Push–pull (E/I) motif. Whether cast as Eversus I, displacement versus momentum, or real versus imaginary coordinates, every oscillator can be decomposed into a push variable xthat drives the system forward and a pull variable ythat drags it back. The canonical Wilson-Cowan loop captures this antagonism, exactly as a harmonic oscillator’s restoring force balances inertia or a Stuart–Landau oscillator’s nonlinear damping balances growth. The same motif powers oscillations in the more complex models. (e) Forcing is any external drive that breaks the autonomy of a node, such as an input from another node in the network (which we call coupling, see next point), electrical stimulation, pharmacological modulation, sensory pulses, or broadband synaptic noise from other sources not explicitly in the model, for example. Weak forcing entrains phase via the phase-response curve;2strong forcing can add or destroy dynamics altogether.3 (f) Coupling allows one neural mass to serve as a dynamic forcing for another. Coupling can yield synchrony, phase slips, amplitude death, or chimera states, depending on delay and gain.4,5 2. Mathematical tools (a) Fixed point: Given a dynamical system ˙x=f(x), a fixed point is a state x∗ that does not change over time; that is f(x∗) = 0. The stability of a fixed point is determined by the eigenvalues of the Jacobian matrix Df. A stable fixed point attracts nearby states, while an unstable one repels nearby states.6 (b) Hopf bifurcation (sometimes referred to as Hopf-Andronov) (HA): From a mathematical perspective, there are 4 different ways in which a two-dimensional system begins, or ceases, to oscillate—for a detailed analysis on how oscillations can arise, we invite the reader to study more detailed and complete analyses.6–8 Here, we will only focus on HA. This bifurcation consists of a fixed point turning stable or unstable through a pair of complex conjugate eigenvalues of Df crossing the imaginary axis. 1 INTRODUCTION 7 (c) Differential Linear operator (ˆ L): Every synapse in the model behaves like a linear filter: it receives an incoming firing-rate trace r(t) and converts it into a post-synaptic potential (PSP) x(t). Written explicitly, this is a linear operation (convolution, see Appendix J) x(t) = ˆ K[r(t)],(1.1) or, equivalently, r(t) = ˆ L[x(t)],(1.2) where ˆ Kis the synapse’s impulse-response kernel and ˆ Lthe inverse operator of ˆ K,ˆ L−1=ˆ K. Framing the dynamics in terms of ˆ Lmakes the filter’s key properties—gain, decay time and delay—visible in a handful of coefficients and shows immediately how the population will amplify, attenuate or phase-shift small perturbations.6In the simplest push–pull oscillator, there are two such operators, one for the excitatory synapse and one for the inhibitory synapse; extending the network merely adds one ˆ Lrow per additional connections between populations. (d) Nonlinearity: Nonlinear functions play a key role in the models beyond the simplest case (the harmonic oscillator). 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. 1 INTRODUCTION 8 Coupled Whole-Brain Models Summary of Neural Mass Model Equations Model Coupled node equation Phase-only (Kuramoto / undamped HO) ˙ θi(t) = ωi+G N X j=1 Cij sinθjt−τij −θi(t)−ιij +ˆ Fe;i(t), i = 1,...,N. Linear damped oscillator network ˙zi(t) = αi+i ωizi(t) + G N X j=1 Cij zjt−τij −zi(t)+ˆ Fe;i(t), i = 1,...,N. Stuart–Landau (Hopf) network ˙zi(t)=(α+i ωi)zi(t)−(γ+i β)|zi(t)|2zi(t) +G N X j=1 Cij zjt−τij −zi(t)+ˆ Fe;i(t), i = 1,...,N. Wilson–Cowan (WILCO) E–I rate network τx˙xi+xi=σx wxx xi−wxy yi+Px,i +ˆ Fe;i(t) + X j=i Cij xj, τy˙yi+yi=σywyx xi−wyy yi, i = 1,...,N. NMM1 (second-order synapses, E–I motif) Lx[xi] = σx−wxy yi+ˆ Fe;i(t) + X j=i Cij xj(t−τij ), Ly[yi] = σywyx xi, i = 1,...,N, where Lα=1 γατ2 α d2 dt2+ 2τα d dt + 1. NMM2 (nextgeneration / QIF-based E–I motif) r(i) x=ˆ ΦxhCxx s(i) x−Cxy s(i) y+ˆ F(i) e(t) + X j=i Cij s(j) xi, s(i) x=ˆ Kxr(i) x, r(i) y=ˆ Φyh−Cyy s(i) y+Cyx s(i) xi, s(i) y=ˆ Kyr(i) y, i = 1,...,N. 2 HARMONIC OSCILLATOR 9 Figure 1.1: Oscillatory dynamics across increasing model complexity. Top row: phase-space portraits. Bottom row: qualitative one-parameter bifurcation diagrams. The horizontal arrow indicates increasing model complexity from left to right. Blue traces highlight attracting (or neutrally stable) periodic orbits and their stable cycle branches; gray traces denote unstable invariant sets and illustrative transients. Panels: (a) phaseonly oscillator ( no amplitude dynamics); (b) damped oscillator with a globally attracting fixed point; (c) Stuart–Landau oscillator (SL) where a Hopf bifurcation creates a stable limit cycle; (d) Wilson–Cowan (WILCO) E–Imodel with coexistence of equilibria and oscillations; (e–f) two neural-mass models (NMM1, NMM2) showing parameter-dependent onset and growth of oscillations and possible bistability. Sketches are schematic and not to scale. 2 Harmonic Oscillator The harmonic oscillator lies at the heart of many rhythmic phenomena in neural systems, providing a mathematically transparent yet conceptually rich framework for understanding how neurons and networks generate, sustain, and modulate oscillations. 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—, and (ii) a mechanism to rotate or cycle continuously through phase space. The simplest realization of these requirements is the undamped (phase-only) oscillator, in which a single phase variable advances uniformly and the amplitude remains fixed. Such a phase reduction captures the essence of lossless oscillatory dynamics and serves as a pedagogical starting point.9,10 Real neural circuits, however, are neither lossless nor isolated: they exhibit intrinsic decay or amplification and receive time-dependent inputs (“forcing”). Accordingly, we proceed in two passes. First, we study the undamped, phase-only oscillator; then, we add external forcing and extend to phase coupling (i.e., Kuramoto paradigm) to examine entrainment and synchronization. Next, we introduce a damping coefficient to allow decay or growth and repeat the sequence—again adding forcing and coupling—to characterize the linear damped oscillator as both a driven resonator and a network element. In both passes, we highlight applications in computational neuroscience. 2 HARMONIC OSCILLATOR 16 From the linear stability analysis, trajectories spiral into (for α < 0) or out of (for α > 0) the off-center equilibrium z∗(see details in Appendix C). In polar coordinates (see Eqs. (2.21)–(2.22)), Re[F e−iθ] can exactly balance the leakage term αr. A constant real Fxensures that x(t) does not decay to zero, mirroring how a steady current injection holds a neuron at a depolarized potential. More generally, time-dependent F(t) 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 Noscillators that, when uncoupled, each satisfy ˙z= (α+i ω)z+ˆ F(t), where α(real) is the common damping (α < 0) or growth (α > 0) rate, ωis the natural frequency, and ˆ F(t) is any external (complex-valued) input. Introducing all-to-all diffusive coupling among these units yields the network equation ˙zi= (α+i ωi)zi+G N N X j=1 zj−zi+ˆ Fe;i(t), i = 1, . . . , N. (2.28) Here, each ωidenotes the intrinsic frequency of node i,G≥0 is the global coupling strength, and the term G NPj(zj−zi) represents a diffusive interaction that pulls each oscillator toward the network centroid 1 NPjzj. 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: ˙zi(t) = αi+i ωizi(t) + G N X j=1 Cij zjt−tij−zi(t)+ˆ Fe;i(t), i = 1, . . . , N. (2.29) Here zi=xi+i yiis the complex state of node i;αi<0 sets linear damping (decay time −1/αi) and ωiis the intrinsic angular frequency; Gis an optional global coupling gain sometimes used to fit models; Cij ≥0 are (possibly directed) connectome weights; τij ≥0 are propagation delays (set τij ≡0 if ignored) and ˆ Fe;i(t) is the external forcing; Setting τij ≡0 and taking Cij =1/N recovers the all-to-all form in Eq. (2.28). 2 HARMONIC OSCILLATOR 17 Coupled Damped oscillator (Eq. 2.29) Whole-brain Simulations Parameters & Physiological Meaning Node iNeural mass (parcel/column or nucleus) represented as a linear damped resonator. NNumber of nodes (parcels). zi(t) = xi+i yiComplex mesoscopic activity; x, y are quadrature “push–pull” components (in-phase / quadrature of a local E/I loop). αiLinear damping (leak/gain). αi<0 gives a stable focus with decay time −1/αi; reflects local E/I balance, membrane and synaptic dissipation. ωiIntrinsic angular frequency (preferred local timescale); set by effective synaptic/membrane constants. fi=ωi/(2π). GGlobal coupling gain controlling the overall strength of inter-areal drive. Cij Connectivity from region jto i(SC from tractography by default; FC/EC alternatives or synthetic graphs if needed). τij Propagation delay on pathway j→i(axonal conduction & synaptic latency); induces frequency-dependent phase offsets. ˆ Fe;i(t) Exogenous drive (external forcing from nodes or elements outside the network or an electric field, see Equation 2.12). Complex or real valued. When to Use It Use this when you want brain-scale resonance with decay (stable focus): fitting FC/covariances, lag structures, and power spectra; probing linear entrainment and susceptibility. Assumptions small fluctuations around equilibrium (αi<0), additive noise and/or weak drive; linear coupling over the connectome (optionally with delays). Best for closed-form second-order statistics (covariances, cross-spectra), rapid parameter sweeps, effective-connectivity estimation, graph-spectral analyses. Avoid if self-sustained oscillations, limit cycles, or multistability are essential—use full Stuart–Landau/Hopf bifurcation dynamics instead. How to Use It Provide C(SC default; FC/EC or synthetic graphs if SC unavailable), G,αi<0, ωi(centered on target band), optional τij ,Fi(t), and ˆηi(t). Defaults normalize C←C/λmax(W); pick αi=−1/τiwith τiin a plausible range; Euler–Maruyama with ∆t≤1/(100 fmax); set τij = 0 initially. Readouts stationary covariance/lagged covariance, PSD and cross-spectra; impulse/transfer functions for linear response; reconstruct xi(t) = Re zi(t) if a real observable is needed. 2.2.3 Applications Equation (2.25) 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 α < 0 sets the decay time back to equilibrium while ωfixes the resonance frequency. Crucially, linearity means that second-order statistics—instantaneous and lagged covariances, cross-spectra and power spectral densities—admit closed-form expressions in terms of the Jacobian and the connectome (with optional delays), so one can predict functional connectivity and spectra without heavy 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.51 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 3 STUART–LANDAU OSCILLATOR (SL) 18 statistics and fast grid searches,51 MOU-based estimation of directed effective connectivity from fMRI,52 and analytic graph-neural-field or spectral-graph models that link Laplacian eigenmodes to MEG/EEG spectra and spatial patterns.53–55 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.56 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 F(t) 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. These principles have been used to interpret phase-specific responses to electrical stimulation—such as phase-locked DBS explained by linearization around a stable focus—and they provide a clean bridge to personalized in-silico perturbations.57 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.58 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.59–63 3 Stuart–Landau Oscillator (SL) Beyond the linear damped oscillator lies the Stuart–Landau (SL) oscillator, a canonical nonlinear model that elegantly captures the transition from decaying or diverging dynamics to robust, self-sustained limit cycles. Purely linear neural mass models cannot adequately describe experimentally observed neuronal oscillations: without nonlinear regulation, excitation would either dominate, leading to runaway growth in firing rates, or inhibition would prevail, extinguishing neural activity altogether. To address this limitation, neural models incorporate biophysical nonlinearities—such as synaptic depression, spike-frequency adaptation, or homeostatic mechanisms—that naturally suppress excessive neuronal activity. Qualitatively, these nonlinear feedback processes act as intrinsic stabilizers, pulling the system back whenever its amplitude begins to exceed physiologically plausible limits. Here we introduce the SL oscillator, exploring its unforced form in Cartesian, polar, and complex representations. Then we discuss how external forcing shapes its dynamics, and introduce its coupled form. We conclude by highlighting SL applications in the context of whole-brain computational models. 3 STUART–LANDAU OSCILLATOR (SL) 19 In polar form writing z=r eiθ yields ˙r=α r −γ r3,(3.1) ˙ θ=ω−β r2.(3.2) Here, the linear terms (α, ω) generate growth and rotation, while the amplitude-dependent nonlinearities (γ, β) self-regulate both amplitude and frequency, stabilizing the limit cycle at r=r∗. More specifically, α∈Ris the linear growth (α > 0) or decay (α < 0) rate, ω > 0 is the intrinsic oscillation frequency, γ > 0 governs nonlinear amplitude saturation, β∈Rintroduces nonlinear frequency modulation, coupling amplitude and phase. Equations ((3.1), (3.2)) represent the normal form of a Hopf bifurcation,64 which gives rise to the sustained oscillations characteristic of the SL (or sometimes called Hopf) model (see Appendix Cfor details). The amplitude equation drives rtoward r∗=rα γ,(3.3) while the phase evolves at an amplitude-dependent rate ˙ θ=ω−β r2. In Cartesian form with x=rcos θand y=rsin θ, one obtains ˙x=α x −ω y −γ(x2+y2)x+β(x2+y2)y, (3.4) ˙y=α y +ω x −γ(x2+y2)y−β(x2+y2)x, (3.5) which makes explicit once again the dominating, push-pull motif in the first-order terms. In complex form, in turn, we have: ˙z= (α+i ω)z−(γ+i β)|z|2z. (3.6) 3.1 Effect of Forcing To model constant or time-varying synaptic inputs, we add a complex forcing term F(t), ˙z= (α+i ω)z−(γ+i β)|z|2z+F(t).(3.7) A constant Fdisplaces the limit cycle to a new fixed point z∗, implicitly given by (α+iω)z∗−(γ+iβ)|z∗|2z∗+F= 0, generally yielding |z∗| =pα/γ. Thus, F(t) 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). 3.2 Coupling Extending to a network of Ndiffusively coupled SL oscillators, each node iobeys ˙zi= (α+i ωi)zi−(γ+i β)|zi|2zi+G N N X j=1 zj−zi+ˆ Fe;i(t),(3.8) 3 STUART–LANDAU OSCILLATOR (SL) 20 where ωiis the intrinsic frequency of node i,Kthe coupling strength, and Fi(t) 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 rich 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: ˙zi(t)=(α+i ωi)zi(t)−(γ+i β)|zi(t)|2zi(t) +G N X j=1 Cij zjt−tij−zi(t)+ˆ Fe;i(t), i = 1, . . . , N. (3.9) where zi=xi+i yiis the complex state; α, γ, β ∈Rand ωi∈Rare the local SL parameters; G≥0 is a global coupling gain; C= [Cij]∈RN×N ≥0is the (possibly directed) connectivity matrix (set Cij = 1/N for all-to-all if no connectome is used); τij ≥0 are propagation delays (set τij ≡0 if ignored) and ˆ Fe;i(t) is the external forcing. We can rewrite these equations in polar coordinates to reconnect with the Kuramoto model65 (see Appendix D). Stuart-Landau Network (Eq. 3.9) Whole-brain Simulations Parameters & Physiological Meaning Node iNeural mass (parcel/column or nucleus) with nonlinear self–limiting dynamics. NNumber of nodes (parcels). zi(t) = xi+i yiComplex mesoscopic state; x, y are quadrature “push–pull” components (E/I-like in-phase and quadrature). ri, θiAmplitude and phase: zi=rieiθi. Isolated steady amplitude (if αi>0): r∗ i=pαi/γi; mean rotation Ωi=ωi−βir∗2 i. αiLinear growth/decay (Hopf bifurcation parameter). αi>0: self-sustained rhythm; αi<0: decay to rest. Physiologically: net local gain/E–I balance, neuromodulatory tone. ωiIntrinsic angular frequency (preferred local timescale) set by effective membrane/synaptic constants. γi>0 Nonlinear amplitude saturation (self-limiting gain); prevents runaway activity, sets r∗ i. βiAmplitude–phase coupling (“shear”): amplitude changes shift instantaneous frequency; captures nonlinear dispersion/adaptation effects. GGlobal coupling gain scaling inter-areal drive over the connectome. Cij Connectivity weight from region jto i(SC by default; FC/EC or synthetic graphs when needed). τij Propagation delay j→i(conduction + synaptic latency); introduces frequency-dependent phase offsets. ˆ Fe;i(t) Exogenous drive (external forcing from nodes or elements outside the network or an electric field, see Equation 2.12). Complex or real valued. When to Use It Use this when you need self-sustained rhythms with amplitude regulation and phase–amplitude coupling; to place nodes near a Hopf bifurcaiton point and study metastability, envelopes, and input-dependent modulation. Coupled form Diffusively coupled SL nodes on a connectome; nonlinear saturation keeps amplitudes bounded and supports synchrony, clustering, chimeras, waves, and amplitude death (all-to-all in (3.8); connectome form in (3.9)). How forcing enters Add a complex input per node to shift the working point, entrain, or reshape amplitude and instantaneous frequency; constant biases displace the attractor, periodic drives yield nonlinear locking (single-node in (3.1)–(3.2); forcing in (3.7)). Strengths Minimal nonlinear phase–amplitude model; explicit control of oscillation onset via α; captures band-limited power, envelope FC/FCD, metastability/turbulence; reduces to linear analytics near the stable focus. Limitations Abstract (few biophysical knobs); parameter identifiability harder far from the Hopf bifurcation; sensitive to delays/coupling very close to the bifurcation; multi-timescale adaptation needs extensions. Typical analyses Fit FC/FCD and spectra; map working point vs. coupling/delays; quantify metastability/turbulence; assess entrainment and phase–amplitude metrics; graph-eigenmode/wave analyses. Good defaults Set γ > 0; place α≈0 (slightly negative for subcritical “fluctuating” regime, slightly positive for limit cycles); modest ωheterogeneity; empirical W(optional delays τij ); small noise σ; tune βto match frequency–amplitude shifts. Report Working point (α), coupling/delay model, frequency distribution, input/noise statistics, and fitted readouts (PSD, FC/FCD, envelopes, metastability). 3 STUART–LANDAU OSCILLATOR (SL) 21 How to Use It Model (drop-in): Provide C(SC default; FC/EC or synthetic graphs if SC unavailable), G, local (αi, ωi, γi>0, βi), optional τij ,Fi(t), and ˆηi(t). Defaults normalize C←Cs/λmax(W); choose small αi(tune distance to Hopf bifurcation: <0 damped, >0 limit cycle), set γi=1 for scaling, βi≈0 initially; initialize zi(0) with small random amplitude; Euler–Maruyama or RK methods with ∆t≤1/(200 fmax). Readouts amplitudes ri(t), phases θi(t), PSD/cross-spectra, spatial/temporal FC, order parameters for phase and amplitude; compare rito r∗ i=pαi/γiwhen αi>0. 3.3 Applications The Stuart–Landau (SL) oscillator traces its roots to Lev Landau’s 1944 phenomenological model of the laminar–turbulent transition, where the complex amplitude Aof an incipient instability obeys ˙ A=σA −β|A|2A, coupling linear growth to cubic saturation so that amplitudes cannot grow without bound. John Trevor Stuart (1958–1960) gave the rigorous weakly nonlinear derivation for parallel shear flows—first in general and then explicitly for plane Poiseuille flow1—showing how amplitude and phase dynamics arise near Hopf bifurcation onset within the Navier–Stokes framework. The resulting SL equation is now understood as the normal form of a Hopf bifurcation: the minimal phase–amplitude model that transitions from damped or growing behavior to a finite-amplitude limit cycle. Its enduring appeal lies in this universality; whenever sustained oscillations emerge from an instability with weak nonlinearity, SL provides a faithful macroscopic scaffold.66–68 In computational neuroscience, these same properties make SL a natural generative model for whole-brain dynamics. At the node level, the bifurcation parameter acontrols the local working point—negative agives noise-driven, damped fluctuations; a≈0 yields large, susceptible excursions; and positive aproduces 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.69 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.70 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.71 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.70 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.72 1Plane Poiseuille flow is a classical fluid–mechanical configuration describing laminar motion of a viscous fluid between two infinitely long, parallel plates. 4 INTERLUDE: SYNAPSES AND TRANSFER FUNCTIONALS 22 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.62,73 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.74 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.75,76 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. 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.77 This system can exhibit a wide array of behaviors—including fixed points and oscillations —depending on the choice of coupling parameters and external inputs. Similar to the phase-only oscillator discussed in Section 2, Wilson–Cowan dynamics unfold in a twodimensional state space. Here, however, the dimensions are the firing rates (x, y) of the two subpopulations, rather than amplitude and phase or Cartesian coordinates. This reflects a more biological perspective on oscillatory phenomena, wherein excitatory and inhibitory pools generate emergent rhythms through their recurrent interactions. 4 Interlude: Synapses and Transfer Functionals Before we embark on describing more biologically realistic models, it is worth discussing in more detail two key elements: the synapse and the transfer functional. The reader will note that the discussion will sometimes refer to neurotransmitter release, ionic currents, or membrane potential perturbations. Insofar as all these quantities are related linearly (by linear operators), we may loosely interchange concepts and equations. That is, the impact of firing rate on a PSP in the receiving cell is mediated by neurotransmitter release, ion channel conduction, and electrical properties of the systems involved — a transformation ladder from an idealized delta function action potential to the dynamics of conductance, synaptic current, and voltage or chloride concentration in the receiving cell. Mathematically, this amounts to the composition of linear operators, which is itself a linear operator (see Appendix Jfor more details on the operator formalism). The reader will see this chain of mechanisms cast as a single first or second-order differential equation. 4 INTERLUDE: SYNAPSES AND TRANSFER FUNCTIONALS 23 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.78 The second is bound to the existence of gap-junctions between nearby cells.79 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, 80–82). Thus, earlier work on NMMs focused on chemical coupling only. Communication through chemical synapses consists on the release of neurotransmitters by the presynaptic neurons whenever they fire. These molecules might bind, then, to the receptors of postsynaptic neurons, which will cause certain ligand-gated ion channels to either open or close, depending on the nature of the neurotransmitter. The flux of ions caused by the synaptic transmission modifies the membrane potential of the postsynaptic neuron. The time scales in which these changes occur vary depending on the neurotransmitters and channels involved, but they generally range between a few to hundreds of milliseconds. This temporal range is similar to that of the neuron’s internal action potential and, therefore, synaptic dynamics play a fundamental role in neural communication. For a single synaptic connection, the dynamics of neurotransmitter binding can be modelled as kinetic reactions.83–85 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 sj(t) (mV) of a single, isolated neuron is usually modeled by the second order ordinary differential equation85,86 τrτd¨sj=γr(t)−(τr+τd) ˙sj−sj(4.1) where τrand τdare the rise and decay times (ms), γis the amplitude of the PSP (mV/kHz), and r(t) is the input term, modeling the arrival of presynaptic spikes. While Eq. (4.1) corresponds to that of a single neuron j, its linearity allows us to quantify the mean synaptic activity of Nneurons receiving and reacting identically to the same inputs as s(t) = 1 NPN j=1 sj(t). Then, we obtain the same equation for the mean synaptic activity, τrτd¨s(t) = γ r(t)−(τr+τd) ˙s−s . (4.2) 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 (r(t) = δ(t)), the solution to (4.2) is s(t) = γ τd−τre−t/τd−e−t/τr.(4.3) Usually, different types of neurotransmitters will result in different timescales. Table 4.1 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, τ=τr=τd. Then Eq. 4.2 reads τ2¨s(t) = γ r(t)−2τ˙s−s . (4.4) 4 INTERLUDE: SYNAPSES AND TRANSFER FUNCTIONALS 24 Neurotransmitter τr[ms] τd[ms] AMPA 0-2 2-5 NMDA 3-15 40-100 GABAA0-2 6-20 GABAB25-50 100-300 Table 4.1: Generic ranges for the values of τrand τd, extracted from86 and references therein. The solution of this equation for a single pulse at time zero (r(t) = δ(t)) reads s(t) = γ t τ2e−t/τ (4.5) which is sometimes referred to as alpha synapse. On the other hand, if one considers that the rise time is almost instantaneous, τr→0, then Eq. 4.2 reads τd˙s=γr(t)−s(t) (4.6) which is just an equation for exponential decay, thus a single pulse at r(t) = δ(t) gives the solution s(t) = γ τd e−t/τd.(4.7) Other types of neurotransmitter kinetics might result in more complex synaptic dynamics.83,84 Of particular importance is the glutamate NMDA receptors, whose dynamics also depend on the voltage of the postsynaptic neuron.85,86 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 ˆ Ls[s(t)] = 1 γτrτd d2 dt + (τr+τd)d dt + 1s(t).(4.8) Then Eq. (4.2) can be written succinctly as ˆ Ls[s(t)] = r(t).(4.9) The synapse operator ˆ Lcan 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 studies87–90 (see Appendix Jfor details). The synapse operation can be cast as a causal filter. Since ˆ Lis a linear differential operator, using proper boundary conditions, we can define its inverse ˆ L−1=ˆ K(another linear operator), s(t) = ˆ Ks[r(t)].(4.10) 4 INTERLUDE: SYNAPSES AND TRANSFER FUNCTIONALS 25 Figure 4.1: Post synaptic potentials in the first-order or the more realistic second-order synapse models. These plots are the response to an “impulse” at time zero, i.e., the solutions to ˆ L[x] = δ(t). 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, ˆ Ls[s(t)] = 1 γ[τdd/dt + 1]s(t).(4.11) From this lens, we view synapses as 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, ˆ L[x] = δ(t−t0). In the case of WilsonCowan, where the operator is first order, the response is simply a decaying exponential (see Fig. 4.1). In more realistic models, such as Jansen-Rit, Wendling, or LaNMM, 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 (τrise) and decay (τ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 meaningful temporal scale, crucial for correctly modeling neuronal integration and timing-dependent neural computations Synaptic dynamics depend on the dynamics of the firing rate of the presynaptic population, which we now turn to. 4.2 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 f-I curve, used to characterize the dynamics of single neurons, to a neural population. 5 THE WILSON–COWAN MODEL (WILCO) 32 Here Cij denotes the (row-normalised) structural connectivity weight from node jto node i, typically estimated from diffusion-MRI tractography. Long–range excitation PjCijxjis 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 σx. Equations (5.5) implement additive coupling: each excitatory node ireceives the transduced firing rates xjof its peers weighted by Cij. 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. It is important to note that Equation 5.5 is one among many potential combination of populations, couplings, and forcings. Wilson-Cowan Network (Eq. 5.5) Whole-brain Simulations Parameters & Physiological Meaning Node iCoarse-grained cortical/subcortical region represented by two interacting subpopulations: excitatory (xi) and inhibitory (yi). NNumber of regions (E–I motifs) in the network. xi(t), yi(t) Population activities (firing rates or PSP proxies, depending on the form used); the E/I push–pull loop produces oscillations and multistability (cf. (5.1), (5.2)). τx, τyMembrane/synaptic time constants (ms) setting local response speeds and the E–I timescale separation. σx(·), σy(·) Static input–output (sigmoid) of each population; typical logistic with slope ραand threshold θα controlling gain and saturation. wxx, wxy, wyx, wyy Local coupling weights (E→E, I→E, E→I, I→I). Signs follow (5.2): −wxy yiimplements inhibition of E by I; −wyyyiself-inhibition of I. Px,i (Py,i)Tonic drives (bias currents) shifting the operating point along the sigmoids and hence the effective linear gains. GGlobal coupling gain scaling long-range excitation from other regions into xi. Cij Long-range (row-normalized) connectivity weight from region jto i(typically structural; functional/effective or synthetic graphs are alternatives). Enters additively inside σxin (5.5). τij Propagation delay on pathway j→i(conduction + synaptic latencies); optional. ˆ Fe;i(t) Exogenous drive (external forcing from nodes or elements outside the network or an electric field, see Equation 2.12). When to Use It Use this when explicit excitation–inhibition, saturating gains, and operating-point control are central (Hopf onset, multistability, stimulus responses, seizure-like dynamics). Assumptions mean-field population description; first-order E/I kinetics; static sigmoids capturing dendritic saturation and threshold dispersion; long-range inputs summed into E. Best for whole-brain simulations with biophysical levers (E/I balance, gains, delays) and macroscopic readouts (FC/FCD, spectra, waves, metastability). Avoid if only phase relations matter (prefer Kuramoto) or small-fluctuation linear analytics suffice (prefer linear damped resonators). How to Use It Provide C(SC by default; FC/EC or synthetic if SC absent), G, local (τx, τy), (w··), σαparameters (ρα, θα), biases Px,i (and optionally Py,i), plus optional τij ,Fi(t), noise. Defaults normalize C←C/λmax(C); choose τx>τyfor gamma-like loops or comparable for alpha/theta; start near Hopf by tuning Px,i and wxx to place the operating point on the steep part of σx; Euler/Heun or RK with ∆t≤1/(100 fmax). Readouts xi(t), yi(t) time series; PSD/cross-spectra; FC/FCD; phase–amplitude metrics; wave/cluster structure; operating-point maps (via linearization around (x0, y0)). 5 THE WILSON–COWAN MODEL (WILCO) 33 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.109,110 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.111–113 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.114 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.115 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.116 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.117,118 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.119,120 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.121 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;122 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.123 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.124–126 Comparative studies situate WILCO among other whole-brain models. Multi-modal head-to-head work reports that WILCO and SL networks achieve broadly compara- 6 NMM WITH SECOND-ORDER SYNAPSES (NMM1) 34 ble 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.127 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.128 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: it trades spiking detail for 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) In the Jansen-Rit or NMM1 formalisma, the level of biological realism is significantly enhanced. NMM1 enriches the Wilson–Cowan and Stuart–Landau E-I motif formalisms by explicitly distinguishing synaptic and somatic stages with secondorder synaptic operators. Rather than collapsing synaptic dynamics into a single firing–rate equation, NMM1 treats each post-synaptic potential x(t) and y(t) as the output of a biologically grounded linear filter ˆ Lτ(with separate rise and decay time constants), sums these to form the membrane perturbation v, and then applies a sigmoidal transfer S(v) to yield the population firing rate. This separation of synapse (impulse-response filters) and soma (sigmoid) provides a clear physiological connection and endows the model with proper delay time scales and intrinsic phase shifts that can sustain oscillations without the need for self-coupling. The biological ontology includes synapses, post-synaptic potentials (PSPs), the population membrane potential (actually the perturbation from its baseline), and the transfer function from membrane potential to firing rate output of the population (the sigmoid). Accordingly, the variables in the equations include the postsynaptic potentials or PSPs (xand yin the E-I model we will discuss), the membrane potential v, the firing rates rxand ry, and the sigmoid S. As usual, all dynamical quantities refer to “mean” population averages. a(We use the term here NMM1 to avoid confusion with the specific model of three populations of Jansen and Rit). 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 6.1 (a), τ2 x¨x+ 2 τx˙x+x=γxσx−wxy y+ˆ Fe, τ2 y¨y+ 2 τy˙y+y=γyσywyx x(6.1) 6 NMM WITH SECOND-ORDER SYNAPSES (NMM1) 35 which can be expressed in operator formalism as Lx[x] = σx−wxy y+ˆ Fe, Ly[y] = σywyx x(6.2) with second order operator notation. For example, for the case of equal rising and decay times, Lα=1 γατ2 α d2 dt2+ 2τα d dt + 1(6.3) 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 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. 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 (m, n) and neuron m, ˆ Lm←num←n=Cm←nrn vm= Λ[ ˆ E] + X n:Cm←n=0 um←n rm=σmvm (6.4) (forcing can be represented through synaptic coupling). As before, the first equation transduces connectivity-weighted firing rate inputs into PSPs using the ˆ Loperator, 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 rnto its associated membrane perturbation (PSP). It can be read as the synapse equation for the input from neuron nto neuron m. The input rnmay also reflect an input to the model from some external neuron, in which case rn=f(t) for some function (constant, noise, etc.). 6 NMM WITH SECOND-ORDER SYNAPSES (NMM1) 36 The second equation evaluates the membrane potential vmof neuron mas the sum of the synaptic voltage perturbations plus an electrical field Eperturbation (if present).129 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.4 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.130 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 Lx[xi] = σx−wxy yi+ˆ Fe;i(t) + X j=i Cij xj(t−τj), Ly[yi] = σywyx xi(6.5) with ˆ Fe(t) 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  E.11–14 Adding this perturbation to the simple example leads to (see Equation 2.12) Lx[xi] = σx−wxy yi+fi(t) +  λi· Ei(t) + ˆηi(t) + N X j=i Cij xj(t−τj), Ly[yi] = σywyx xi (6.6) NMM1 Network (Eq. 6.5) Whole-brain Simulations Parameters & Physiological Meaning Node iNeural mass with explicit synapse and soma: PSP states drive a membrane perturbation vi, which passes through a sigmoid to yield a firing rate ri. xi(t), yi(t) Excitatory and inhibitory postsynaptic potentials (PSPs). They are the outputs of second-order synaptic filters (rise/decay), not directly the rates. vi(t) Membrane potential perturbation (sum of PSPs and exogenous terms), i.e., the input to the soma/nonlinearity. rx,i(t), ry,i(t) Excitatory/inhibitory firing rates (soma outputs). These feed other synapses locally and across the network. σx(·), σy(·) Static sigmoids (e.g., logistic) mapping vto firing rate; slope controls effective gain; saturation bounds activity. Lx, LySecond-order synaptic operators (cf. (6.3)): implement biophysical PSP kinetics with rise/decay; supply intrinsic phase lags that can sustain oscillations. τα,(τα,r, τα,d) Synaptic time constants (single or separate rise/decay); set resonance frequency and phase lag of each synapse. γαSynaptic gains (PSP amplitude scale). wyx, wxy Local E→I and I→E coupling (push–pull loop); self-coupling often omitted in NMM1 because synaptic phase lags can close the Barkhausen loop. ˆ Fe;i(t) Exogenous drive (external forcing from nodes or elements outside the network or an electric field)— see Equation 2.12). Cij Long-range connectivity from node jto i(SC default; FC/EC or synthetic graphs if needed); typically targets the excitatory pathway. NNumber of nodes (parcels). 6 NMM WITH SECOND-ORDER SYNAPSES (NMM1) 37 Figure 6.1: Four models using the second-order formalism: a) PING-like push–pull motif, b) Jansen–Rit model,131,132 c) Wendling model,133 d) Laminar model.134,135 Synapses are shown as arrowheads (excitatory) or buttons (inhibitory). Noise or external inputs are indicated on pyramidal cells, although other targets (µ) are also possible. When to Use It Use this when explicit synaptic kinetics and their phase lags matter (evoked responses, resonance, photic entrainment, band-limited power and envelope dynamics) or when linking parameters to PSP amplitudes/time-constants is essential. Assumptions linear synaptic filters (second order) feeding a static nonlinearity; E/I pathways combined at the soma; long-range excitation enters via excitatory synapses. Best for EEG/MEG/fMRI generative modeling (spectra and ERPs), seizure phenomenology (fast/slow inhibition variants), and whole-brain simulations where synaptic time constants set rhythms. Relations near a stable focus it reduces to linear resonators; near Hopf it displays SL-like amplitude–phase dynamics but with biophysical PSP knobs (gains/time-constants). Avoid if only relative phase is of interest (Kuramoto) or if closed-form linear statistics suffice (damped linear network). How to Use It Model (drop-in) use the operator form (6.4) or the E–I pair with coupling (6.5)–(6.6). Provide synaptic operators Lα(choose ταor (τα,r, τα,d)) and gains γα; local couplings wyx, wxy; soma sigmoids σα (gain, midpoint, max rate); drives Fi(t) (and optionally  λ· Ei). Network supply C(SC default; FC/EC or synthetic if SC is absent), optional delays τij ; long-range input enters the excitatory pathway. Defaults normalize C←C/λmax(C); start with standard JR-style kinetics (faster E than I or vice-versa depending on band); small noise; Euler–Maruyama/RK with ∆t≤1/(200 fmax). Readouts PSPs (xi, yi), membrane vi, rates ri; PSD/cross-spectra and ERPs; envelope/FC/FCD; assess resonance by scanning τ’s and gains. 6.2 Applications We review canonical models using the second-order synapse formalism (see Figure 6.1) 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 E→I gain and can be shifted by drive or kinetics. Applications include modeling high-frequency oscillations (HFOs) at seizure onset;136 mechanistic context from spiking/mean-field work on PING/ING is reviewed in Buzs´aki & Wang (2012)137 and Tiesinga & Sejnowski (2009),138 and synchronization analyses such as in B¨orgers & Kopell (2003)139 and Whittington et al. (2000).140 6 NMM WITH SECOND-ORDER SYNAPSES (NMM1) 38 The Jansen–Rit model The Jansen–Rit (JR) model92 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.141,142 JR serves as the neuronal model in Dynamic Causal Modeling (DCM) for M/EEG steadystate and evoked responses,143–145 linking synaptic gains/time constants to observed spectra and ERPs. The Wendling model and its extensions Wendling’s CA1-inspired extension adds fast and slow inhibitory subpopulations (GABAA, GABAB) 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.146 The model is widely used for interpreting SEEG/EEG patterns, exploring ictogenesis mechanisms, and assessing interventions,147 making it a standard computational tool in epilepsy. Recent extensions include chloride dynamics and laminar integration.148 Whole-brain use: Wendling nodes have been used for resting-state whole-brain network modeling that matches empirical FC149 and for patient-specific whole-brain simulations of interictal SEEG to aid clinical interpretation.150 They have been embedded in realistic forward models to synthesize SEEG and personalize local epileptogenic dynamics;148 the same group reports personalized whole-brain seizure-propagation models integrating SEEG, MRI and dMRI.151 The laminar model (LaNMM) LaNMM embeds laminar-specific projections and volume-conduction physics to connect mesoscopic generators to depth LFP/CSD. Practically, it combines a deep JR-like slow generator (alpha/theta) with a superficial PING-like fast generator (gamma), coupled according to cortical laminar anatomy, yielding coexisting slow/fast rhythms with realistic depth-dependent polarity and phase relations. It reproduces laminar spectral peaks and CSD sinks/sources and accounts for cross-frequency interactions observed in laminar recordings152,153 and the oscillatory features of Alzheimer’s disease.154 The framework also supports integration with stimulation physics for tES/tACS studies and has been used to study predictive coding155 mechanisms and the effects of psychedelics in AD.156 It has also been deployed at the connectome scale to study gamma-band coordination and cooperative/competitive network rhythms.157 7 NEXT-GENERATION MODELS (NMM2) 39 7 Next-generation models (NMM2) The neural mass models we have discussed so far, such as the Jansen-Rit, Wendling systems or LaNMM, provide a practical framework for representing and interpreting electrophysiological activity in both local and global brain models.102,131,152,154,158–166 However, they are only partly derived from first principles. While the post-synaptic potential (PSPs) dynamics are inferred from data and can be grounded on diffusion physics,106,167,168 Freeman’s “wave to pulse” sigmoid function,169–171 used to transduce mean population membrane potential into firing rate, rests on a weaker theoretical standing. More importantly, it is far from obvious how the activity of a large number of neurons can be coarse-grained into mesoscale models with many fewer degrees of freedom. Recently, Montbri´o et al172 derived an exact mean-field theory (MPR) for a population of quadratic integrate-and-fire neurons under some simplifying assumptions, thereby connecting microscale neural mechanisms and meso/macroscopic phenomena. The MPR model can be seen to replace Freeman’s sigmoid function with a pair of differential equations for the mean membrane potential and firing rate variables—a dynamical relation between firing rate and membrane potential—, providing a more fundamental interpretation of the semi-empirical NMM sigmoid parameters. In doing so, it sheds light on the mechanisms behind enhanced network response to weak but uniform perturbations. In the exact mean-field theory, intrinsic population connectivity modulates the steady-state firing rate sigmoid relation in a monotonic manner, with increasing excitatory self-connectivity leading to higher firing rates. This provides a plausible mechanism for the enhanced response of densely connected networks to weak, uniform inputs such as the electric fields produced by non-invasive brain stimulation. This new, “dynamic sigmoid” also endows the neural mass model with a form of “inertia”, an intrinsic delay to external inputs that depends on, e.g., self-coupling strength and state of the system. Models resulting from the MPR mean-field theory can be completed by adding the first or second-order equations for delayed post-synaptic currents and the coupling term with an external electric field,173–175 bringing together the MPR and the usual NMM formalisms into a unified exact mean-field theory (NMM2, for short) displaying rich dynamical features. In the single population model, we show that the resonant sensitivity to a weak alternating electric field is enhanced by increased self-connectivity and slow synapses. Thus, the NMM2 framework further elevates biological grounding by deriving the equations from first principles. It replaces the static “wave-to-pulse” sigmoid with an exact mean-field firing-rate dynamics derived from quadratic integrate-and-fire (QIF) neurons. 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.176,177 Empirically, data-driven coarse-graining of population activity can approach non-Gaussian fixed forms with static 7 NEXT-GENERATION MODELS (NMM2) 40 and dynamic scaling—evidence that mesoscale statistics can be approximately scale-invariant near special operating points.178 Early NMMs (e.g., Wilson–Cowan, Jansen–Rit) are phenomenological: the linear synaptic filter is biophysically grounded, but the wave-to-pulse static nonlinearity is heuristic. Population-density and mean-field limits provide a more principled route from spiking to masses,179,180 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.181 In large-scale modeling, these ideas motivate dynamic mean-field reductions of biophysical networks into mesoscale nodes with a few state variables.182 A major step forward is the exact reduction of all-to-all QIF networks with heterogeneous excitabilities to two macroscopic ODEs for population firing rate r(t) and mean membrane potential v(t) (via the QIF–θtransform and the Ott–Antonsen manifold). This is the Montbri´o–Paz´o–Roxin (MPR) theory;183 see also reviews and extensions.174,175,184–187 In the MPR framework, the steady-state relationship between firing rate rand mean voltage vdefines a sigmoid-type input–output curve (i.e., a static transfer function). However, when the full two-dimensional dynamical system is considered (with r(t) and v(t) evolving in time), the effective gain becomes stateand history-dependent, rather than being a fixed static curve. Positive self-coupling (J↑) 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).183,188 Exact mean-field extensions capture synaptic filtering, electrical coupling, and other biophysics,173–175 and have been leveraged to model working-memory circuits with short-term plasticity.189 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 RG language, rand vare 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 (η, ∆, J) and mesoscale responsiveness, and (iii) a natural path to incorporate fluctuation corrections when needed (finite-size, correlations).181 In this sense, NMM2 is a next-generation neural mass that is both biophysically anchored and explicitly multiscale.182,187 In their uniform, mean-field derivation for a population of Nquadratic integrate and fire (QIF) neurons, Monbri´o et al172 start from the equations for the neuron membrane potential perturbation from baseline, ˙ Vj=V2 j+ηj+Js(t) + I(t),if Vj≥Vp,then Vr←Vj(7.1) In this equation, the total input current in neuron jis Ij=ηj+Js(t) + I(t) and includes a quenched noise constant component ηjdrawn from a Lorentzian (Cauchy) distribution, the input from other neurons s(t) per connection received (the mean synaptic activation) with uniform coupling J, and a common input I(t). The common input I(t) can represent both a common external input or the effect of an electric field, e.g., I(t) = p(t) +  λ· E(t) (7.2) 7 NEXT-GENERATION MODELS (NMM2) 41 for weak electric fields. Here p(t) 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.,190 and191). The mean synaptic activation is given by s(t) = 1 N N X j=1 X k|tk j<t Zt −∞ dt′aτ(t−t′)δ(t′−tk j) (7.3) where tk jis the arrival time of the kth spike from the jth neuron, and a(t) the synaptic activation function, e.g., a(t) = e−t/τ /τ. Note that we can write J=j·N, where jis 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 Nof them in a fully connected architecture with Nneurons). 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 p(t) +  λ· E(t) + J s(t) 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 Nlimit (Eq 12 in172), ˙r= ∆/π + 2rv (7.4) ˙v=v2−π2r2+J s + ¯η+I(t) (7.5) To this equation, we add the usual operator dynamics for the synapse activation, ˆ Ls[s] = r, or, equivalently, s=ˆ Ks[r] (7.6) (see Figure 6.1 (A)). Here vand rare the population mean membrane potential and firing rate, respectively. The new parameters ¯ηand ∆ refer to the mean and half-width of the Lorentzian distribution for the quenched noise input ηj. The analysis in172 hinges on the assumptions of all-to-all uniform connectivity (with synaptic weight J) and common input I(t). These equations can be read as a transfer functional mapping input currents into an output firing rate, as in the master NMM Equation 4.19, r=ˆ Φm[J s +I(t)] , s=ˆ Ks[r].(7.7) The transfer functional is captured by Equations 7.4. REFERENCES 48 to NMM1 and NMM2—offers a modular blueprint for building, interpreting, and extending neural mass models. Each rung distills a distinct set of assumptions about synapses, transfer functions, and coupling while remaining dynamically and conceptually connected to the others. Future work will likely move up and down this ladder: deriving more realistic masses from increasingly detailed neuron models; enriching local circuits with laminar, cell-type and synaptic complexity; and scaling up to heterogeneous, directed whole-brain networks that incorporate subcortical structures and predictive processing. Our hope is that such a unified dynamical language will make it easier to compare models, design experiments and interventions, and eventually bridge microcircuit physiology, whole-brain activity, and cognition within a single coherent framework. 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 GarciaOjalvo was also financially supported by the 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. References [1] Steven H. Strogatz. Nonlinear Dynamics and Chaos. Addison–Wesley, 1994. [2] Bard Ermentrout and Nancy Kopell. Multiple traveling waves in neuronal networks. SIAM Journal on Applied Mathematics, 51:179–194, 1991. [3] Eugene M. Izhikevich. Dynamical Systems in Neuroscience: The Geometry of Excitability and Bursting. MIT Press, 2007. [4] Arkady Pikovsky, Michael Rosenblum, and J¨urgen Kurths. Synchronization: A universal concept in nonlinear sciences. Cambridge Nonlinear Science Series, 12, 2003. [5] Daniel M. Abrams and Steven H. Strogatz. Chimera states for coupled oscillators. Physical Review Letters, 93:174102, 2004. [6] Steven H Strogatz. Nonlinear dynamics and chaos with student solutions manual. (No Title), 2018. REFERENCES 49 [7] Frank C Hoppensteadt and Eugene M Izhikevich. Weakly connected neural networks, volume 126. Springer Science & Business Media, 2012. [8] Yuri A. Kuznetsov. Elements of Applied Bifurcation Theory. Applied Mathematical Sciences. Springer Cham, 4 edition, 2023. [9] Steven H. Strogatz. Nonlinear Dynamics and Chaos: With Applications to Physics, Biology, Chemistry, and Engineering. Westview Press, Boulder, CO, 2nd edition, 2018. [10] Arthur T. Winfree. The Geometry of Biological Time. Springer, New York, NY, 2nd edition, 2001. [11] Giulio Ruffini, Fabrice Wendling, Isabelle Merlet, Behnam Molaee-Ardekani, Abeye Mekonnen, Ricardo Salvador, Aureli Soria-Frisch, Carles Grau, Stephen Dunne, and Pedro C. Miranda. Transcranial Current Brain Stimulation (tCS): Models and Technologies. IEEE Transactions on Neural Systems and Rehabilitation Engineering, 21(3):333–345, May 2013. Conference Name: IEEE Transactions on Neural Systems and Rehabilitation Engineering. [12] Giulio Ruffini, Michael D Fox, Oscar Ripolles, Pedro Cavaleiro Miranda, and Alvaro Pascual-Leone. Optimization of multifocal transcranial current stimulation for weighted cortical pattern targeting from realistic modeling of electric fields. Neuroimage, 89:216–225, 2014. [13] Adri`a Galan-Gadea, Ricardo Salvador, Fabrice Bartolomei, Fabrice Wendling, and Giulio Ruffini. Spherical harmonics representation of the steady-state membrane potential shift induced by tDCS in realistic neuron models. Journal of Neural Engineering, 20(2), March 2023. [14] Aman S. Aberra, Angel V. Peterchev, and Warren M. Grill. Biophysically realistic neuron models for simulation of cortical stimulation. Journal of Neural Engineering, 15(6):066023, December 2018. [15] Aman S. Aberra, Boshuo Wang, Warren M. Grill, and Angel V. Peterchev. Simulation of transcranial magnetic stimulation in head model with morphologicallyrealistic cortical neurons. Brain Stimulation, 13(1):175–189, January 2020. [16] David C. Roberts. Linear reformulation of the kuramoto model of self-synchronizing oscillators. Physical Review E, 77(3):031114, 2008. [17] Laurie Conteville and Elena Panteley. Linear reformulation of the kuramoto model: Asymptotic mapping and stability properties. In Proceedings of the 2013 European Control Conference (ECC), pages 821–826, 2013. [18] Yoshiki Kuramoto. Chemical Oscillations, Waves, and Turbulence, volume 19 of Springer Series in Synergetics. Springer, Berlin, Heidelberg, 1984. [19] Hiroya Nakao. Phase reduction approach to synchronization of nonlinear oscillators. Contemporary Physics, 57(2):188–214, 2016. [20] Bastian Pietras, Federico Devalle, Alex Roxin, Andreas Daffertshofer, and Ernest Montbri´o. Exact firing rate model reveals the differential effects of chemical versus electrical synapses in spiking networks. Phys. Rev. E, 100:042412, Oct 2019. REFERENCES 50 [21] Yoshiki Kuramoto. Self-entrainment of a population of coupled non-linear oscillators. In Huzihiro Araki, editor, International Symposium on Mathematical Problems in Theoretical Physics, pages 420–422, Berlin, Heidelberg, 1975. Springer. [22] Steven H. Strogatz. Nonlinear dynamics and chaos: with applications to physics, biology, chemistry, and engineering. Boulder, CO: Westview Press, a member of the Perseus Books Group, 2015. [23] Juan A. Acebr´on, L. L. Bonilla, Conrad J. P´erez Vicente, F´elix Ritort, and Renato Spigler. The Kuramoto model: A simple paradigm for synchronization phenomena. Reviews of Modern Physics, 77(1):137–185, April 2005. [24] G. Bard Ermentrout and David H. Terman. Mathematical Foundations of Neuroscience. Springer, New York, NY, 2010. [25] Steven H. Strogatz. From Kuramoto to Crawford: Exploring the onset of synchronization in populations of coupled oscillators. Physica D: Nonlinear Phenomena, 143(1):1–20, September 2000. [26] Arthur T. Winfree. Biological rhythms and the behavior of populations of coupled oscillators. Journal of Theoretical Biology, 16(1):15–42, 1967. [27] Yoshiki Kuramoto. Lecture notes in physics, international symposium on mathematical problems in theoretical physics. In Lecture Notes in Physics, volume 30 of Lecture Notes in Physics. 1975. [28] Steven H. Strogatz. From kuramoto to crawford: Exploring the onset of synchronization in populations of coupled oscillators. Physica D: Nonlinear Phenomena, 143(1–4):1–20, 2000. [29] J. A. Acebr´on, L. L. Bonilla, C. J. P´erez Vicente, F. Ritort, and R. Spigler. The kuramoto model: A simple paradigm for synchronization phenomena. Reviews of Modern Physics, 77(1):137–185, 2005. [30] Michael Breakspear, Stewart Heitmann, and Andreas Daffertshofer. Generative models of cortical oscillations: neurobiological implications of the kuramoto model. Frontiers in Human Neuroscience, 4:190, 2010. [31] Joana Cabral, Etienne Hugues, Olaf Sporns, and Gustavo Deco. Role of local network oscillations in resting-state functional connectivity. NeuroImage, 57(1):130– 139, 2011. [32] Joana Cabral, Henry Luckhoo, Mark Woolrich, Morten Joensson, Hamid Mohseni, Adam Baker, Morten L. Kringelbach, and Gustavo Deco. Exploring mechanisms of spontaneous functional connectivity in meg: How delayed network interactions lead to structured amplitude envelopes of band-pass filtered oscillations. NeuroImage, 90:423–435, 2014. [33] J. Cabral, E. Hugues, M. L. Kringelbach, and G. Deco. Role of local network oscillations in resting-state functional connectivity. NeuroImage, 57:130–139, 2014. [34] Gy¨orgy Buzs´aki. Rhythms of the Brain. Oxford University Press, Oxford, UK, 2006. REFERENCES 51 [35] Joana Cabral, Emmanuel Hugues, Olaf Sporns, and Gustavo Deco. Role of local network oscillations in resting-state functional connectivity. NeuroImage, 57(1):130– 139, 2011. [36] Ruben Schmidt, Karl J. R. LaFleur, Marcel A. de Reus, Leonard H. van den Berg, and Martijn P. van den Heuvel. Kuramoto model simulation of neural hubs and dynamic synchrony in the human cerebral connectome. BMC Neuroscience, 16:54, 2015. [37] Murray Shanahan. Metastable chimera states in community-structured oscillator networks. Chaos, 20(1):013108, 2010. [38] Johanne Hizanidis, Nikos E. Kouvaris, Gorka Zamora-L´opez, Albert D´ıaz-Guilera, and Chris G. Antonopoulos. Chimera-like states in modular neural networks. Scientific Reports, 6:19845, 2016. [39] Joana Cabral, Morten L. Kringelbach, and Gustavo Deco. Exploring the network dynamics underlying brain activity during rest. Progress in Neurobiology, 114:102– 131, 2014. [40] Changgui Gu, Zonghua Liu, William J. Schwartz, and Premananda Indic. Photic desynchronization of two subgroups of circadian oscillators in a network model of the suprachiasmatic nucleus with dispersed coupling strengths. PLoS ONE, 7(5):e36900, 2012. [41] Spase Petkoski, Andreas Spiegler, Timoth´ee Proix, Parham Aram, Jean-Jacques Temprado, and Viktor K. Jirsa. Heterogeneity of time delays determines synchronization of coupled oscillators. Physical Review E, 94(1):012209, 2016. [42] Spase Petkoski, J. Matias Palva, and Viktor K. Jirsa. Phase-lags in large-scale brain synchronization: Methodological considerations and in-silico analysis. PLoS Computational Biology, 14(7):e1006160, 2018. [43] Spase Petkoski, Andreas Spiegler, Timoth´ee Proix, Parham Aram, Jean-Jacques Temprado, and Viktor K. Jirsa. Transmission time delays organize the brain network synchronization. Philosophical Transactions of the Royal Society A, 377(2153):20180132, 2019. [44] P. Fries. Rhythms for Cognition: Communication through Coherence. Neuron, 88(1):220–35, October 2015. [45] Peter A. Tass. Desynchronization by means of a coordinated reset of neural subpopulations. Progress of Theoretical Physics Supplement, 150:281–296, 2003. [46] Oleksandr V. Popovych, Serhiy Yanchuk, and Peter A. Tass. Multisite delayed feedback for electrical brain stimulation. Frontiers in Neuroscience, 12:644, 2018. [47] Gihan Weerasinghe, Hayriye Cagnan, Peter Brown, and Rafal Bogacz. Predicting the effects of deep brain stimulation using a reduced coupled oscillator model. PLoS Computational Biology, 15(7):e1006575, 2019. [48] Gihan Weerasinghe, Hayriye Cagnan, Peter Brown, and Rafal Bogacz. Optimal closed-loop deep brain stimulation using multiple independently controlled contacts. PLoS Computational Biology, 17(8):e1009281, 2021. REFERENCES 52 [49] Maximilian Sadilek and Stefan Thurner. Physiologically motivated multiplex kuramoto model describes phase diagram of cortical activity. Scientific Reports, 5:10015, 2015. [50] Lukas G. Bauer, Fabian Hirsch, Craig Jones, Michael Hollander, Philipp Grohs, Amit Anand, Chris Plant, and Andreas Wohlschl¨ager. Quantification of kuramoto coupling between intrinsic brain networks applied to fmri data in major depressive disorder. Frontiers in Computational Neuroscience, 16:729556, 2022. [51] Adri´an Ponce-Alvarez and Gustavo Deco. The hopf whole-brain model and its linear approximation. Scientific Reports, 14:2615, 2024. [52] Matthieu Gilson, Gustavo Deco, Katsumi Friston, Patric Hagmann, Dante Mantini, Michael L. L. M. van Mieghem, D. J. O. de Pasquale, Viktor Jirsa, Andrea J. R. Mess´e, and H. R. W. R. Senden. Estimation of directed effective connectivity from fmri functional connectivity hints at asymmetries of cortical connectome. PLOS Computational Biology, 12(1):e1004762, 2016. [53] Marco Aqil, Selen Atasoy, Morten L. Kringelbach, and Rikkert Hindriks. Graph neural fields: A framework for spatiotemporal dynamical models on the human connectome. PLOS Computational Biology, 17(1):e1008310, 2021. [54] Ashish Raj, Chang Cai, Xihe Xie, Eva Palacios, Julia Owen, Pratik Mukherjee, and Srikantan S. Nagarajan. Spectral graph theory of brain oscillations. Human Brain Mapping, 41(11):2980–2998, 2020. [55] Parul Verma, Srikantan S. Nagarajan, and Ashish Raj. Spectral graph theory of brain oscillations—revisited and improved. NeuroImage, 249:118919, 2022. [56] Erfan Nozari, Maxwell A. Bertolero, Jennifer Stiso, Lorenzo Caciagli, Eli J. Cornblath, Xiaosong He, Arun S. Mahadevan, George J. Pappas, and Dani S. Bassett. Macroscopic resting-state brain dynamics are best described by linear models. Nature Biomedical Engineering, 8:68–84, 2024. [57] Benoit Duchet and Rafal Bogacz. How to design optimal brain stimulation to modulate phase-amplitude coupling? Journal of neural engineering, 21(4):10.1088/1741– 2552/ad5b1a, July 2024. [58] Felix Effenberger, Pedro Carvalho, Igor Dubinin, and Wolf Singer. The functional role of oscillatory dynamics in neocortical circuits: A computational perspective. Proceedings of the National Academy of Sciences, 122(4):e2412830122, January 2025. Publisher: Proceedings of the National Academy of Sciences. [59] Gustavo Deco, Christopher W. Lynn, Yonatan Sanz Perl, and Morten L. Kringelbach. Violations of the fluctuation–dissipation theorem reveal distinct nonequilibrium dynamics of brain states. Physical Review E, 108(6):064410, 2023. [60] Joana Cabral, Gustavo Deco, and Morten L. Kringelbach. Remote synchronization in the human cortex. Scientific Reports, 14(1):5678, 2024. [61] Gustavo Deco, Juan Cruzat, Joana Cabral, Enzo Tagliazucchi, Helmut Laufs, Nikos K. Logothetis, and Morten L. Kringelbach. Modeling the impact of lsd on whole-brain dynamics using the hopf model. NeuroImage, 218:116871, 2020. REFERENCES 53 [62] Gustavo Deco and Morten L. Kringelbach. Turbulent-like dynamics in the human brain. Cell Reports, 33(10):108471, December 2020. Publisher: Elsevier BV. [63] Adri´an Ponce-Alvarez and Gustavo Deco. The Hopf whole-brain model and its linear approximation. Scientific Reports, 14(1):2615, January 2024. Publisher: Nature Publishing Group. [64] Jurij A. Kuznecov. Elements of applied bifurcation theory. Number 112 in Applied mathematical sciences. Springer, New York Berlin Heidelberg, 1995. [65] Ana P. Mill´an, David Poyato, David N. Reynolds, and Francesco Tudisco. Synchronization of coupled Stuart-Landau oscillators: How heterogeneity can facilitate synchronization, October 2025. [66] L. D. Landau. On the problem of turbulence. Dokl. Akad. Nauk SSSR, 44(8):339– 349, 1944. [67] J. T. Stuart. On the non-linear mechanics of hydrodynamic stability. Journal of Fluid Mechanics, 4(1):1–21, 1958. [68] J. T. Stuart. On the non-linear mechanics of wave disturbances in stable and unstable parallel flows. part 1. the basic behaviour in plane poiseuille flow. Journal of Fluid Mechanics, 9(3):353–370, 1960. [69] Gustavo Deco, Morten L. Kringelbach, Viktor K. Jirsa, and Petra Ritter. The dynamics of resting fluctuations in the brain: metastability and its dynamical cortical core. Scientific Reports, 7(1):3095, 2017. [70] Joana Cabral, Francesca Castaldo, Jakub Vohryzek, Vladimir Litvak, Christian Bick, Renaud Lambiotte, Karl Friston, Morten L. Kringelbach, and Gustavo Deco. Metastable oscillatory modes emerge from synchronization in the brain spacetime connectome. Communications Physics, 5:184, 2022. [71] Gustavo Deco, Joana Cabral, Mark W. Woolrich, Angus B. A. Stevner, Tim J. van Hartevelt, and Morten L. Kringelbach. Single or multiple frequency generators in on-going brain activity: A mechanistic whole-brain model of empirical meg data. NeuroImage, 152:538–550, 2017. [72] Francesca Castaldo, Francisco P´ascoa Dos Santos, Ryan C Timms, Joana Cabral, Jakub Vohryzek, Gustavo Deco, Mark Woolrich, Karl Friston, Paul Verschure, and Vladimir Litvak. Multi-modal and multi-model interrogation of large-scale functional brain networks. NeuroImage, 277:120236, 2023. [73] Yonatan Sanz Perl, Anira Escrichs, Enzo Tagliazucchi, Morten L. Kringelbach, and Gustavo Deco. Strength-dependent perturbation of whole-brain model working in different regimes reveals the role of fluctuations in brain dynamics. PLOS Computational Biology, 18(11):e1010662, 2022. [74] Anira Escrichs, Yonatan Sanz Perl, Camila Uribe, et al. Unifying turbulent dynamics framework distinguishes different brain states. Communications Biology, 5:638, 2022. [75] Beatrice M. Jobst, Selen Atasoy, Adri´an Ponce-Alvarez, Ana Sanju´an, Leor Roseman, Mendel Kaelen, Robin Carhart-Harris, Morten L. Kringelbach, and Gustavo REFERENCES 54 Deco. Increased sensitivity to strong perturbations in a whole-brain model of lsd. NeuroImage, 230:117809, 2021. [76] Josephine Cruzat, Yonatan Sanz Perl, Anira Escrichs, Jakub Vohryzek, Christopher Timmermann, Leor Roseman, Andrea I. Luppi, Agustin Iba˜nez, David Nutt, Robin Carhart-Harris, Enzo Tagliazucchi, Gustavo Deco, and Morten L. Kringelbach. Effects of classic psychedelic drugs on turbulent signatures in brain dynamics. Network Neuroscience (Cambridge, Mass.), 6(4):1104–1124, 2022. [77] Adri´an Ponce-Alvarez and Gustavo Deco. The hopf whole-brain model and its linear approximation. Scientific Reports, 14:2615, 2024. [78] D. Hillis, H.C. Heller, S.D. Hacker, D. Hall, D. Sadava, and M. Laskowski. Life: The Science of Biology. Macmillan Learning, 2020. [79] Daniel A. Goodenough and David L. Paul. Gap Junctions. Cold Spring Harbor Perspectives in Biology, 1(1):a002576, July 2009. [80] Alberto E. Pereda. Electrical synapses and their functional interactions with chemical synapses. Nature Reviews Neuroscience, 15(4):250–263, April 2014. [81] Pepe Alcam´ı and Alberto E. Pereda. Beyond plasticity: The dynamic impact of electrical synapses on neural circuits. Nature Reviews Neuroscience, 20(5):253–271, May 2019. [82] Mitchell J. Vaughn and Julie S. Haas. On the Diverse Functions of Electrical Synapses. Frontiers in Cellular Neuroscience, 16, 2022. [83] Alain Destexhe, Zachary F. Mainen, and Terrence J. Sejnowski. Synthesis of models for excitable membranes, synaptic transmission and neuromodulation using a common kinetic formalism. Journal of Computational Neuroscience, 1(3):195–230, August 1994. [84] A. Destexhe, Z. F. Mainen, and T. J Sejnowski. Kinetic Models of Synaptic Transmission. In Methods in Neuronal Modeling: From Ions to Networks, pages 1–25. MIT Press, Cambridge, MA, 2nd edition, 1998. [85] Peter Dayan and L. F. Abbott. Theoretical Neuroscience: Computational and Mathematical Modeling of Neural Systems. Computational Neuroscience. Massachusetts Institute of Technology Press, Cambridge, Mass, 2001. [86] Wulfram Gerstner, Werner M. Kistler, Richard Naud, and Liam Paninski. Neuronal Dynamics: From Single Neurons to Networks and Models of Cognition. Cambridge University Press, first edition, July 2014. [87] W. Rall, R. E. Burke, T. G. Smith, P. G. Nelson, and K. Frank. Dendritic location of synapses and possible mechanisms for the monosynaptic EPSP in motoneurons. Journal of Neurophysiology, 30(5):1169–1193, September 1967. [88] R. E. Burke. Composite nature of the monosynaptic excitatory postsynaptic potential. Journal of Neurophysiology, 30(5):1114–1137, September 1967. [89] P. G. Nelson and K. Frank. Anomalous rectification in cat spinal motoneurons and effect of polarizing currents on excitatory postsynaptic potential. Journal of Neurophysiology, 30(5):1097–1113, September 1967. REFERENCES 55 [90] W. Rall. Distinguishing theoretical synaptic potentials computed for different somadendritic distributions of synaptic input. Journal of Neurophysiology, 30(5):1138– 1168, September 1967. [91] Hugh R. Wilson and Jack D. Cowan. Excitatory and Inhibitory Interactions in Localized Populations of Model Neurons. Biophysical Journal, 12(1):1–24, 1972. [92] Ben H. Jansen and Vincent G. Rit. Electroencephalogram and visual evoked potential generation in a mathematical model of coupled cortical columns. Biological Cybernetics, 73(4):357–366, September 1995. [93] Arnold J. F. Siegert. On the first passage time probability problem. Phys. Rev., 81:617–623, Feb 1951. [94] D. Amit. Model of global spontaneous activity and local structured activity during delay periods in the cerebral cortex. Cerebral Cortex, 7(3):237–252, April 1997. [95] Nicolas Brunel and Vincent Hakim. Fast Global Oscillations in Networks of Integrate-and-Fire Neurons with Low Firing Rates. Neural Computation, 11(7):1621–1671, October 1999. [96] Nicolas Fourcaud-Trocm´e, David Hansel, Carl van Vreeswijk, and Nicolas Brunel. How Spike Generation Mechanisms Determine the Neuronal Response to Fluctuating Inputs. Journal of Neuroscience, 23(37):11628–11640, December 2003. [97] Nicolas Brunel and Peter E. Latham. Firing Rate of the Noisy Quadratic Integrateand-Fire Neuron. Neural Computation, 15(10):2281–2306, October 2003. [98] Ernest Montbri´o, Diego Paz´o, and Alex Roxin. Macroscopic Description for Networks of Spiking Neurons. Physical Review X, 5(2):021028, June 2015. [99] Pau Clusella and Ernest Montbri´o. Regular and sparse neuronal synchronization are described by identical mean field dynamics, 2022. [100] Walter J. Freeman. Mass Action in the Nervous System. Elsevier, 1975. [101] Walter J. Freeman. Nonlinear gain mediating cortical stimulus-response relations. Biological Cybernetics, 33(4):237–247, 1979. [102] Hugh R Wilson and Jack D Cowan. Excitatory and Inhibitory interactions in localized populations of model neurons. Biophysical Journal, 12(1):1–24, 1972. [103] Michael Breakspear. Dynamic models of large-scale brain activity. Nature Neuroscience, 20(3):340–352, 2017. [104] Lea Fredrickson-Hemsing, Seung Ji, Robijn Bruinsma, and Dolores Bozovic. Modelocking dynamics of hair cells of the inner ear. Physical Review E, 86(2):021915, August 2012. Publisher: American Physical Society. [105] John Guckenheimer and Philip Holmes. Nonlinear oscillations, dynamical systems, and bifurcations of vector fields. Number 42 in Applied mathematical sciences. Springer Science+Business Media, New York, NY, corrected seventh printing edition, 2002. [106] G. Ermentrout and David H. Bard, Terman. Mathematical Foundations of Neuroscience. Springer-Verlag New York, 2010. REFERENCES 56 [107] Bastian Pietras and Andreas Daffertshofer. Network dynamics of coupled oscillators and phase reduction techniques. Physics Reports, 819:1–105, 2019. [108] Gustavo Deco and et al. The dynamic brain: from spiking neurons to neural masses and cortical fields. PLOS Computational Biology, 4(8), August 2008. [109] Hugh R. Wilson and Jack D. Cowan. Excitatory and inhibitory interactions in localized populations of model neurons. Biophysical Journal, 12(1):1–24, 1972. [110] Hugh R. Wilson and Jack D. Cowan. A mathematical theory of the functional dynamics of nervous tissue. Kybernetik, 13:55–80, 1973. [111] Alain Destexhe and Terrence J. Sejnowski. The wilson–cowan model, 36 years later. Biological Cybernetics, 101(1):1–2, 2009. [112] Jack D. Cowan, Jeremy Neuman, and Wim van Drongelen. Wilson–cowan equations for neocortical dynamics. Journal of Mathematical Neuroscience, 6(1):1, 2016. [113] Carson C. Chow and Yasaman Karimipanah. Before and beyond the wilson–cowan equations. Journal of Mathematical Neuroscience, 10(1):61, 2020. [114] Romesh G. Abeysuriya, Jonathan Hadida, Stamatios N. Sotiropoulos, Saad Jbabdi, Robert Becker, Benjamin A. E. Hunt, Matthew J. Brookes, and Mark W. Woolrich. A biophysical model of dynamic balancing of excitation and inhibition in fast oscillatory large-scale networks. PLOS Computational Biology, 14(2):e1006007, 2018. [115] S. L. Keeley et al. Firing rate models for brain rhythms. Journal of Neurophysiology, 122(4):1561–1588, 2019. [116] Guoshi Li, Yujie Liu, Yanting Zheng, Ye Wu, Danian Li, Xinyu Liang, Yaoping Chen, Ying Cui, Pew-Thian Yap, Shijun Qiu, Han Zhang, and Dinggang Shen. Multiscale neural modeling of resting-state fmri reveals executive–limbic malfunction as a core mechanism in major depressive disorder. NeuroImage: Clinical, 31:102758, 2021. [117] Paula Sanz Leon, Stuart A. Knock, Andreas Spiegler, and Viktor K. Jirsa. The virtual brain: a simulator of primate brain network dynamics. Frontiers in Neuroinformatics, 7:10, 2013. [118] Paula Sanz-Le´on, Stuart A. Knock, Andreas Spiegler, and Viktor K. Jirsa. Mathematical framework for large-scale brain network modeling in The Virtual Brain. NeuroImage, 111:385–430, 2015. [119] Thomas F. Varley, Prejaas Tewarie, et al. Non-reversibility as a quantitative marker of neurophysiological change from meg. NeuroImage, 273:120002, 2023. Uses WC local dynamics as a generative benchmark. [120] Rustam Masharipov, Dongha Shin, Qisheng Zhao, Quinn K. Telesford, Paula SanzLeon, et al. Task-modulated functional connectivity: a comprehensive evaluation of state-of-the-art methods. Communications Biology, 7(1):???, 2024. Includes WC generative simulations. [121] Sadjad Sadeghi, Daniela Mier, Martin F. Gerchen, Susanne N. L. Schmidt, and Joachim Hass. Dynamic causal modeling for fmri with wilson–cowan–based neuronal equations. Frontiers in Neuroscience, 14:593867, 2020. REFERENCES 57 [122] Hil G. E. Meijer, Tahra L. Eissa, Bert Kiewiet, Jeremy F. Neuman, Catherine A. Schevon, Robert G. Emerson, Robert R. Goodman, Guy M. McKhann, Charles J. Marcuccilli, Andrew K. Tryba, Jack D. Cowan, Stephan A. van Gils, and Wim van Drongelen. Modeling focal epileptic activity in the wilson–cowan model with depolarization block. Journal of Mathematical Neuroscience, 5(1):7, 2015. [123] Laura M. S´anchez-Rodr´ıguez, Miriam Vila-Vidal, Perrine Berroir, Arnau M´endez, Gustavo Deco, Morten L. Kringelbach, and Joana Cabral. Bridging scales in alzheimer’s disease by whole-brain modelling: from the braak staging scheme to functional connectivity. Communications Biology, 7(1):???, 2024. [124] Antonio de Candia, Alessandro Sarracino, Ilenia Apicella, and Lucilla de Arcangelis. Critical behaviour of the stochastic wilson–cowan model. PLOS Computational Biology, 17(8):e1008884, 2021. [125] Ilenia Apicella, Antonio de Candia, Daniele Conte, Lucilla de Arcangelis, and Alessandro Sarracino. Power spectrum and critical exponents in the 2d stochastic wilson–cowan model. Scientific Reports, 12(1):22612, 2022. [126] Hamed Alvankar Golpayegan, Andrea Cavagna, Serena di Santo, ´ Edgar Rold´an, Massimiliano Viale, Alessio Attanasi, and Massimiliano Viale. Bistability and criticality in the stochastic wilson–cowan model. Physical Review E, 107(3):034404, 2023. [127] Francesca Castaldo, Maarten De Vos, Jakub Vohryzek, Joana Cabral, Gustavo Deco, Morten L. Kringelbach, et al. Multi-modal and multi-model interrogation of large-scale functional brain networks. NeuroImage, 275:120236, 2023. [128] J. W. M. Domhof, S. B. Eickhoff, and O. V. Popovych. Reliability and subject specificity of personalized whole-brain dynamical models. NeuroImage, 257:119321, 2022. [129] Giulio Ruffini, Michael D. Fox, Oscar Ripolles, Pedro Cavaleiro Miranda, and Alvaro Pascual-Leone. Optimization of multifocal transcranial current stimulation for weighted cortical pattern targeting from realistic modeling of electric fields. NeuroImage, 89:216 – 225, 2014. [130] Borja Mercadal, Maria Guasch-Morgades, Lucia Mencarelli, Giacomo Koch, and Giulio Ruffini. Bridging local and global dynamics: a biologically grounded model for cooperative and competitive interactions in the brain. bioRxiv, pages 2025–07, 2025. Publisher: Cold Spring Harbor Laboratory. [131] B. H. Jansen, G. Zouridakis, and M. E. Brandt. A neurophysiologically-based mathematical model of flash visual evoked potentials. Biol Cybern, 68(3):275–83, 1993. [132] F. Grimbert and Olivier Faugeras. Analysis of Jansen’s model of a single cortical column. INRIA, RR-5597:34, June 2006. [133] Fabrice Wendling, Pascal Benquet, Fabrice Bartolomei, and Viktor Jirsa. Computational models of epileptiform activity. Journal of Neuroscience Methods, 260:233 – 251, 2016. REFERENCES 64 [213] Liang Chen and Sue Ann Campbell. Exact mean-field models for spiking neural networks with adaptation. Journal of Computational Neuroscience, 50(4):445–469, November 2022. [214] Alberto Ferrara, David Angulo-Garcia, Alessandro Torcini, and Simona Olmi. Population spiking and bursting in next-generation neural masses with spike-frequency adaptation. Physical Review E, 107(2):024311, February 2023. [215] Bastian Pietras, Pau Clusella, and Ernest Montbri´o. Low-dimensional model for adaptive networks of spiking neurons. Phys. Rev. E, 111:014422, Jan 2025. [216] Edward Ott and Thomas M. Antonsen. Low dimensional behavior of large systems of globally coupled oscillators. Chaos: An Interdisciplinary Journal of Nonlinear Science, 18(3):037113, 2008. [217] Bastian Pietras, Rok Cestnik, and Arkady Pikovsky. Exact finite-dimensional description for networks of globally coupled spiking neurons. Phys. Rev. E, 107:024315, Feb 2023. [218] Jorge F. Mejias, John D. Murray, Henry Kennedy, and Xiao-Jing Wang. Feedforward and feedback frequency-dependent interactions in a large-scale laminar network of the primate cortex. Science Advances, 2(11):e1601335, November 2016. [219] Ryan V Raut, Anish Mitra, Scott Marek, Mario Ortega, Abraham Z Snyder, Aaron Tanenbaum, Timothy O Laumann, Nico U F Dosenbach, and Marcus E Raichle. Organization of Propagated Intrinsic Brain Activity in Individual Humans. 30(3):1716– 1734. [220] Andre M. Bastos, W. Martin Usrey, Rick A. Adams, George R. Mangun, Pascal Fries, and Karl J. Friston. Canonical microcircuits for predictive coding. Neuron, 76(4):695, November 2012. [221] Karl Friston. The free-energy principle: a unified brain theory? Nature Reviews Neuroscience, 11(2):127–138, February 2010. Publisher: Springer Science and Business Media LLC. [222] Matthew F. Glasser, Timothy S. Coalson, Emma C. Robinson, Carl D. Hacker, John Harwell, Essa Yacoub, Kamil Ugurbil, Jesper Andersson, Christian F. Beckmann, Mark Jenkinson, Stephen M. Smith, and David C. Van Essen. A multi-modal parcellation of human cerebral cortex. Nature, 536(7615):171–178, August 2016. [223] Morten L Kringelbach, Josephine Cruzat, Joana Cabral, Gitte Moos Knudsen, Robin Carhart-Harris, Peter C Whybrow, Nikos K Logothetis, and Gustavo Deco. Dynamic coupling of whole-brain neuronal and neurotransmitter systems. Proceedings of the National Academy of Sciences, 117(17):9566–9576, 2020. Publisher: National Acad Sciences. [224] Rodrigo Cofr´e, Rub´en Herzog, Pedro A.M. Mediano, Juan Piccinini, Fernando E. Rosas, Yonatan Sanz Perl, and Enzo Tagliazucchi. Whole-Brain Models to Explore Altered States of Consciousness from the Bottom Up. Brain Sciences, 10(9):626, September 2020. [225] Murat Demirta¸s, Joshua B. Burt, Markus Helmer, Jie Lisa Ji, Brendan D. Adkinson, Matthew F. Glasser, David C. Van Essen, Stamatios N. Sotiropoulos, Alan REFERENCES 65 Anticevic, and John D. Murray. Hierarchical Heterogeneity across Human Cortex Shapes Large-Scale Neural Dynamics. Neuron, 101(6):1181–1194.e13, March 2019. Publisher: Elsevier. [226] Joshua B. Burt, Murat Demirta¸s, William J. Eckner, Natasha M. Navejar, Jie Lisa Ji, William J. Martin, Alberto Bernacchia, Alan Anticevic, and John D. Murray. Hierarchy of transcriptomic specialization across human cortex captured by structural neuroimaging topography. Nature Neuroscience, 21(9):1251–1259, September 2018. Publisher: Nature Publishing Group. [227] Eve Marder. Neuromodulation of Neuronal Circuits: Back to the Future. Neuron, 76(1):1–11, October 2012. Publisher: Elsevier. [228] Gustavo Deco, Yonatan Sanz Perl, Jakub Vohryzek, Andrea Luppi, and Morten L. Kringelbach. Evolution’s boldest trick: Neurotransmission modulated whole-brain computation captures full task repertoire, June 2025. Pages: 2025.06.02.657368 Section: New Results. [229] George Datseris and Ulrich Parlitz. Nonlinear Dynamics: A Concise Introduction Interlaced with Code. Undergraduate Lecture Notes in Physics. Springer International Publishing, Cham, 2022. [230] E. J. Doedel and others. AUTO-07P: Continuation and bifurcation software for ordinary differential equations. manual, Concorida University, Montreal, Canada, 2007. [231] Raul de Palma Aristides, Pau Clusella, Roser Sanchez-Todo, Giulio Ruffini, and Jordi Garcia-Ojalvo. Emergence of multifrequency activity in a laminar neural mass model. arxiv, 2025. tex.affiliation: Department of Medicine and Life Sciences, Universitat Pompeu Fabra, Barcelona, Spain; Department of Mathematics, Universitat Polit`ecnica de Catalunya, Manresa, Spain; Center of Brain and Cognition, Universitat Pompeu Fabra, Barcelona, Spain; Brain Modeling Department, Neuroelectrics, Barcelona, Spain. [232] Frank C. Hoppensteadt and Eugene M. Izhikevich. Weakly Connected Neural Networks, volume 126 of Applied Mathematical Sciences. Springer New York, New York, NY, 1997. [233] Louis M. Pecora and Thomas L. Carroll. Master stability functions for synchronized coupled systems. Phys. Rev. Lett., 80(10):2109–2112, 1998. [234] H. Fujisaka and T. Yamada. Stability theory of synchronized motion in coupledoscillator systems. Prog. Theor. Phys., 69(1):32–47, 1983. [235] M. Barahona and L. M. Pecora. Synchronization in small-world systems. Phys. Rev. Lett., 89(5):054101, 2002. [236] R. M. May. Will a large complex system be stable? Nature, 238:413–414, 1972. [237] D. G. Aronson, G. B. Ermentrout, and N. Kopell. Amplitude response of coupled oscillators. Physica D, 41(3):403–449, 1990. [238] D. V. Ramana Reddy, Abhijit Sen, and G. L. Johnston. Time delay induced death in coupled limit cycle oscillators. Phys. Rev. Lett., 80(23):5109–5112, 1998. REFERENCES 66 [239] Arkady Pikovsky, Michael Rosenblum, and J¨urgen Kurths. Synchronization: A Universal Concept in Nonlinear Sciences. Cambridge University Press, 2001. [240] Giulio Ruffini, Francesca Castaldo, and Jakub Vohryzek. Structured dynamics in the algorithmic agent. Entropy, 27(1):90, 2025. [241] Jack Carr. Applications of Centre Manifold Theory, volume 35 of Applied Mathematical Sciences. Springer-Verlag, New York, 1981. [242] Stephen Wiggins. Introduction to Applied Nonlinear Dynamical Systems and Chaos, volume 2 of Texts in Applied Mathematics. Springer, New York, 2003. [243] Michael J. Field. Dynamics and Symmetry, volume 3 of ICP Advanced Texts in Mathematics. Imperial College Press, London, 2007. [244] Viktor Jirsa and Hiba Sheheitli. Entropy, free energy, symmetry and dynamics in the brain. Journal of Physics: Complexity, 3(1):015007, February 2022. Publisher: IOP Publishing. [245] Iv´an Le´on and Hiroya Nakao. Analytical Phase Reduction for Weakly Nonlinear Oscillators. Chaos, Solitons & Fractals, 176:114117, November 2023. arXiv:2308.02105 [nlin]. [246] Erwin Kreyszig, Herbert Kreyszig, and Edward Norminton. Advanced Engineering Mathematics. John Wiley & Sons, Hoboken, NJ, 10 edition, 2011. [247] Erol Ba¸sar. Brain oscillations in neuropsychiatric disease. Dialogues in Clinical Neuroscience, 15(3):291–300, 2013. [248] Thomas Donoghue, Mikael Haller, Ethan J. Peterson, Pankaj Varma, Patrick Sebastian, Richard Gao, Tory Noto, Andrea H. Lara, Jonathan D. Wallis, Robert T. Knight, Andrey Shestyuk, and Bradley Voytek. Parameterizing neural power spectra into periodic and aperiodic components. Nature Neuroscience, 23(12):1655–1665, 2020. [249] Steven H. Strogatz. Nonlinear Dynamics and Chaos : With Applications to Physics, Biology, Chemistry, and Engineering. Addison–Wesley, Reading, MA, 1994. [250] B. O. Koopman. Hamiltonian systems and transformation in hilbert space. Proceedings of the National Academy of Sciences, 17(5):315–318, 1931. [251] T. Andrew Whitten, Andrew M. Hughes, Courtney T. Dickson, and Jeremy B. Caplan. A better oscillation detection method robustly extracts eeg rhythms across brain state changes: The human alpha rhythm as a test case. NeuroImage, 54(2):860–874, 2011. [252] Thomas M. Cover and Joy A. Thomas. Elements of information theory. John Wiley & sons, 2 edition, 2006. [253] B. van der Pol. On “relaxation-oscillations”. The London, Edinburgh, and Dublin Philosophical Magazine and Journal of Science, 2(11):978–992, 1926. [254] J. T. Stuart. On the non-linear mechanics of hydrodynamic stability. Journal of Fluid Mechanics, 4:1–21, 1958. REFERENCES 67 [255] Richard FitzHugh. Impulses and physiological states in theoretical models of nerve membrane. Biophysical Journal, 1(6):445–466, 1961. [256] Catherine Morris and Harold Lecar. Voltage oscillations in the barnacle giant muscle fiber. Biophysical Journal, 35(1):193–213, 1981. [257] Alan J. McKane and Timothy J. Newman. Predator–prey cycles from resonant amplification of demographic stochasticity. Physical Review Letters, 94(21):218102, 2005. [258] Ming Li and Paul M. B. Vitanyi. Applications of algorithmic information theory. Scholarpedia, 2(5):2658, 2007. [259] Giulio Ruffini and Edmundo Lopez-Sola. AIT foundations of structured experience. Journal of Artificial Intelligence and Consciousness, 9(2):153–191, September 2022. Publisher: World Scientific Pub Co Pte Ltd. [260] Giulio Ruffini. Structured dynamics in the algorithmic agent, December 2023. Pages: 2023.12.12.571311 Section: New Results. [261] John Guckenheimer and Philip Holmes. Nonlinear oscillations, dynamical systems, and bifurcations of vector fields, volume 42. Springer Science & Business Media, 2013. A TERMINOLOGY AND SCOPE OF AGGREGATE NEURONAL MODELS 68 A Terminology and Scope of Aggregate Neuronal Models Under the umbrella of aggregate neuronal modeling (also called neural mass,population, or lumped models), we encompass a spectrum of formalisms that trade biological detail for analytical or computational simplicity. This hierarchy can be traversed in both directions—adding realism to derive more mechanistic descriptions or stripping back complexity to reveal core dynamical principles: •Phase oscillator: describes the evolution of a single phase variable θ(t), capturing limit-cycle interactions at the most abstract level. •Damped oscillator: introduces amplitude relaxation, e.g., a second-order linear system with friction. •Stuart–Landau (SL) normal form: a generic nonlinear oscillator near a Hopf bifurcation, used phenomenologically to model neural rhythms. •Wilson–Cowan (WILCO): the prototypical two-population rate (population/lumped/neuralmass) model, describing mean excitatory and inhibitory firing with coupled ODEs. •Neural Mass Model 1 (NMM1): adds biophysical filtering of post-synaptic potentials (PSPs) and explicit conversion from mean membrane potentials to firing rates. •Neural Mass Model 2 (NMM2 or MPR): extends NMM1 by incorporating dynamics in the transfer function. When we take the limit of infinitely many, infinitesimally small populations (or let the spatial coupling kernel become continuous), this lumped description generalizes to integrodifferential or partial-differential equations known as neural field models. Aggregate neuronal models (phase oscillators, SL, WILCO, NMM1, NMM2, etc.) are fundamentally statistical constructs that provide macro-level, effective descriptions of complex networks of many neurons. All share these core features: 1. Mean variables: each model tracks ensemble averages—phase and amplitude (phase oscillator, SL), firing rate (WILCO), membrane potential and PSP (NMM1), plus synchronization metrics (NMM2). 2. Aggregation by type: neurons are grouped into populations based on anatomical and functional characteristics, each treated as a single “lumped” unit. 3. Effective parameters: time-constants, gains and connectivities summarize average behavior; they do not map one-to-one onto single-cell properties. 4. Scale invariance: the dynamical equations predict the same trajectories regardless of neuron count, provided effective parameters (e.g. mean synapses per neuron) remain fixed. 5. Bridge to measurements: to relate outputs (phase, rate, or voltage) to LFP, currentsource density, or BOLD signals, morphological and density-based scaling factors must be applied post-hoc. B COORDINATE TRANSFORMATIONS 69 B Coordinate Transformations This appendix provides the full mathematical derivations of how one obtains Cartesian and complex representations from polar form (and vice versa) for both the undamped and damped oscillators. B.1 Undamped Oscillator: Polar ←→ Cartesian ←→ Complex Polar to Cartesian. Starting from the polar equations ˙ θ=ω, (B.1) ˙r= 0,(B.2) define x=rcos θ, y =rsin θ. (B.3) Differentiate xand ywith respect to time. Since r=r0is constant, ˙x=d dtrcos θ=−rsin θ˙ θ=−rsin θ(ω) = −ω(rsin θ) = −ω y, ˙y=d dtrsin θ=rcos θ˙ θ=rcos θ(ω) = ω(rcos θ) = ω x. Hence, ˙x=−ω y, ˙y=ω x, (B.4) which recovers the undamped Cartesian form (2.1)–(2.2). Cartesian to Complex. Given x(t) and y(t) satisfying (B.4), set z=x+i y. (B.5) Then ˙z= ˙x+i˙y= (−ω y) + i(ω x) = i ω (x+i y) = i ω z, recovering the undamped complex form (2.8). Conversely, writing z=r ei θ and differentiating yields ˙z= ˙r ei θ +r ei θ i˙ θ= ( ˙r+i r ˙ θ)ei θ. By matching ˙z=i ω z =i ω r ei θ, one finds ˙r= 0,˙ θ=ω, which recovers (B.2)–(B.1). B COORDINATE TRANSFORMATIONS 70 B.2 Damped Oscillator: Polar ←→ Cartesian ←→ Complex Polar to Cartesian. Begin with the damped polar equations: ˙r=α r +Fr(t),(B.6) ˙ θ=ω+Fθ(t) r.(B.7) Define x=rcos θ, y =rsin θ, (B.8) and differentiate: ˙x= ˙rcos θ−rsin θ˙ θ, ˙y= ˙rsin θ+rcos θ˙ θ. Substitute ˙rand ˙ θfrom (B.6)–(B.7): ˙x=α r +Fr(t)cos θ−rsin θω+Fθ(t) r =α(rcos θ) + Fr(t) cos θ−ω(rsin θ)−sin θ Fθ(t) =α x −ω y +Fr(t) cos θ−Fθ(t) sin θ. Similarly, ˙y=α r +Fr(t)sin θ+rcos θω+Fθ(t) r =α(rsin θ) + Fr(t) sin θ+ω(rcos θ) + cos θ Fθ(t) =α y +ω x +Fr(t) sin θ+Fθ(t) cos θ. Define Ix(t) = Fr(t) cos θ−Fθ(t) sin θ, Iy(t) = Fr(t) sin θ+Fθ(t) cos θ. (B.9) Then ˙x=−ω y +α x +Ix(t),(B.10) ˙y=ω x +α y +Iy(t),(B.11) which recover the damped Cartesian form (2.23)–(2.24). Cartesian to Complex. From (B.10)–(B.11), let z=x+i y and Iz(t) = Ix(t)+i Iy(t). Then ˙z= ˙x+i˙y =−ω y +α x +Ix(t)+iω x +α y +Iy(t) =α(x+i y) + i ω (x+i y)+(Ix+i Iy) = (α+i ω)z+Iz(t), yielding the damped-and-forced complex form (2.25). B COORDINATE TRANSFORMATIONS 71 Conversely, writing z=r ei θ, ˙z= ˙r ei θ +i r ˙ θ ei θ =˙r+i r ˙ θei θ. Matching ˙z= (α+i ω)z+Iz(t)=(α+i ω)r ei θ +Iz(t), one identifies ˙r=α r + Ree−i θIz(t), r˙ θ=ω r + Ime−i θIz(t), so that, by defining Fr(t) = Re(e−i θIz(t)) and Fθ(t) = Im(e−i θIz(t)), one recovers the polar form (B.6)–(B.7). Notation and Labels in Appendix: •The external inputs Ix(t), Iy(t) in Cartesian form are related to radial/tangential forcings Fr(t), Fθ(t) via (B.9). •The complex input Iz(t) is Ix(t) + i Iy(t). B.3 WILCO and second order equations Unlike phase–amplitude models that admit a natural polar form, the Wilson-Cowan system is most naturally expressed in a 2D (x, y) phase plane. We can interpret X= x y!,˙ X= ˙x ˙y!, and rewrite the Wilson-Cowan equations in vector form: τx˙x τy˙y= −x+Sx(wxx x−wxy y+Ix) −y+Sy(wyx x−wyy y+Iy)!. Phase-plane analysis (nullclines, fixed points, and limit cycles) is used to study how X evolves in R2. Second order equations. Because each equation is second-order, one can rewrite them as pairs of first–order ODEs. For instance, let x1=xand x2= ˙x. Then:      ˙x1=x2, ˙x2=1 τ2 xh−2τxx2+KxSx···−x1i. Similarly for (y1, y2). Alternatively, one may keep the original form and perform phase– space analysis in four dimensions (x, ˙x, y, ˙y). The −xand −yterms (right–hand side) act as a linear decay, while the second–order operator allows for resonant or damped oscillatory responses that are further modulated by the nonlinear saturations Sxand Sy. C LINEAR STABILITY AND BIFURCATION ANALYSIS 72 C Linear Stability and Bifurcation Analysis C.1 Linear stability analysis A feature that all neural mass models discussed in main text have in common, except for the harmonic oscillator, is that they are all nonlinear models. Although most nonlinear systems are impractical to solve analytically and the use of numerical methods are often required, we can still gain insights in the systems behaviors through linearization methods. As the literature in this topic is very rich, we will briefly review main concepts that will help the reader understand the dynamical behavior of the neural mass models discussed in the main text. First, is it worth to recall that a N-dimensional linear system, which can be written as ˙ x(t) =      a11 a12 ... a1N a21 a22 ... a2N . . .. . ..... . . aN1aN2... aNN           x1(t) x2(t) . . . xN(t)      =Ax(t) (C.1) Notice that at each instant tthe matrix Alinearly transforms the vector x. To obtain the solution of this systems we can rely on the fact that Acan be diagonalized, hence it can be expressed as: A=PΛP−1(C.2) where Pis the matrix whose columns are the eigenvectors of A, and Λ is the diagonal matrix whose entries are the eigenvalues {λi}of A. Substituting we have ˙x =PΛP−1x(t).(C.3) Introducing z(t) = P−1x(t), we have ˙ z(t) = Λz(t) =      λ1z1(t) λ2z2(t) . . . λNzN(t)      .(C.4) And have a solution of the form zi(t) = cieλit, thus z(t) =      c1eλ1t c2eλ2t . . . cNeλNt      .(C.5) Finally, using that x(t) = Pz(t), we get x(t) = c1eλ1tv1+c2eλ2tv2+. . . +cNeλNtvN(C.6) where v1,v2,...,vNare the corresponding eigenvectors of Aand c1, c2, . . . , cNare constants determined by the initial conditions of the system. In this form, it is clear that the solution and hence the dynamics of the system will be dictated by its eigenvalues, which C LINEAR STABILITY AND BIFURCATION ANALYSIS 73 can be real or complex. For example, if we take a two-dimensional system with ci>0, we have x(t) = c1eλ1tv1+c2eλ2tv2(C.7) If both λ1and λ2are positive, the system diverges as both exponentials grow to infinity; on the other hand, if both are negative, the exponentials tend to zero. Notice how the magnitude of the eigenvalues will determine the speed with which the system evolves in each direction. Interestingly, if the eigenvalues form a complex conjugate, λ1,2=α±iβ which allow solutions of the type x(t) = c1e(α1±iβ1)t+c2e(α2±iβ2)t=c1eα1tcos(β1t)v1+c2eα2tcos(β2t)v2(C.8) and hence oscillatory behavior. Of special interest to us is the behavior of the system around fixed points x∗. Fixed points are defined by solutions to ˙ x(t) = F(x) = 0 (C.9) which represent equilibrium solutions, since if the system starts at x∗, it will stay there if it is kept unperturbed. Now, if we take a nonlinear system ˙x =F(x) (C.10) where Frepresents a nonlinear transformation of the vector xand usually cannot be diagonalized, we can not obtain a solution with form of Eq. (C.6). Fortunately, we can linearize nonlinear systems around their fixed points, which at first order yield us ˙x ≈F(x∗) + DF(x)|x∗(x−x∗) + O2.(C.11) And by using F(x∗) = 0, discarding higher order terms we get: ˙x ≈DF(x)|x∗(x−x∗).(C.12) This linear version of our original system, allow us to use the linear analysis tools to get an insight on the dynamical properties of the nonlinear system: the eigenvalues λof the Jacobian matrix DFwill encoded the behavior of system around the fixed points. Below, we summarize the different types of fixed points in Table C.1. While we focus on two-dimensional systems for clarity, the reasoning applies to N-dimensional systems: the eigenvalues of the Jacobian still determine the local stability and type of equilibrium. The classification in higher dimensions, though, involves more combinations of stable, unstable, and center directions. Readers interested in the general case can consult standard references on nonlinear dynamics such as.8,22,229 While fixed points correspond to states where the system remains constant over time, some nonlinear systems exhibit limit cycles, in which the system evolves along a closed, repeating trajectory in phase space. Limit cycles represent self-sustained oscillations and can be either stable, attracting nearby trajectories, or unstable, repelling them. In this way, limit cycles extend the concept of equilibrium from a single point to a periodic orbit. Their stability are encoded by Floquet exponents. Beyond understanding the local behavior around fixed points and limit cycles, it is equally important to study how changes in the system’s parameters affect its qualitative dynamics. Such transitions occur when fixed points (the solutions of Eq. (C.9)) or limit cycles C LINEAR STABILITY AND BIFURCATION ANALYSIS 80 Figure C.7: Bifurcation diagram of the PING (NMM2) model with bifurcation parameter I. To summarize, bifurcation diagrams are particularly valuable for neural mass models because they: •Reveal the onset of oscillations and other dynamical regimes. They identify where fixed points lose stability and give rise to limit cycles, as well as where bistability or other qualitative transitions occur, providing a clear picture of the model’s possible behaviors. •Guide parameter selection and sensitivity analysis. By showing how solutions depend on parameters such as coupling strengths or synaptic gains, bifurcation diagrams help identify which parameters critically shape the dynamics and which regimes are physiologically plausible. •Inform control and intervention strategies. In applications ranging from neuromodulation to pharmacology, knowing how small parameter changes can shift the system between regimes is essential. Bifurcation diagrams make these transitions explicit by highlighting where stability changes occur. D PHASE DYNAMICS AND KURAMOTO MODEL 81 D Phase dynamics and Kuramoto model In the main text, we begin with a pure phase oscillator model and gradually introduce increasing mathematical complexity. Nonetheless, coupled phase oscillator models, particularly the Kuramoto model, were originally derived from arbitrary oscillator systems subject to weak perturbations by employing rigorous mathematical techniques. Such approach, usually known as phase reduction, was pioneered by Winfree and Kuramoto, and has been widely applied in mathematical neuroscience to capture the phase dynamics of both single-neuron models and neural population models. For a comprehensive treatment of these methods, readers may refer to.4,10,18,232 Drawing on this literature, in this Appendix we provide a concise introduction to the phase reduction technique in a general setting, and then illustrate it using the Stuart–Landau oscillator, which naturally leads to the Kuramoto model. D.1 Phase of a perturbed oscillator Let us consider an m-dimensional variable x∈Rm, whose time evolution is ruled by the differential equation ˙ x=F(x) (D.1) where Fis a smooth function. Let us assume that equation (D.1) has an attracting limit-cycle with period T,xLC(t) = xLC(t+T). Thus, the curve xis, by definition, parameterized by time, t∈[0, T). We define the phase of a point x(t) of xLC as ϕ(x(t)) := t2π T∈[0,2π).(D.2) For simplicity, we denote ω:= 2π T. The concept of phase is well defined for any point lying exactly on the limit cycle. Nevertheless, the phase reduction approach is based on extending the notion of phase to a neighborhood G⊂Rmof xLC(t). Let y0=y(0) ∈Rmbe a point in the state space close to, but not exactly on the periodic attractor xLC. Since the limit cycle is stable, the trajectory y(t) will end up arbitrarily close to xLC as t→ ∞, thus will have a well defined phase, ϕ(y(t)). Therefore, even if y0is not on the invariant manifold, one can find a point, x0∈xLC such that, at long-term, ϕ(y(t)) = ϕ(x(t)). One can then define the phase of y0as ϕ(y0) = ϕ(x0). All the points in the neighborhood of xLC associated to the same phase form a isochrone curve, i.e., the isochrones are the level curves of the function ϕ. This extension of the notion of phase to the vicinity of the periodic orbit preserves a constant time evolution: ˙ ϕ(x) = ω , x∈G. (D.3) Applying the chain rule and using Eq. D.1, one obtains ˙ ϕ(x) = ∇xϕ·F(x).(D.4) Combining the two previous equations provides ∇xϕ·F(x) = ω . (D.5) D PHASE DYNAMICS AND KURAMOTO MODEL 82 This equation expresses then, that the variation of the phase under the effect of the field Fis constant. We want to study oscillatory systems that are interacting with a certain environment, so, usually, equation (D.1) is not enough for our purposes. One needs to account for the effect that external perturbations have on the system. Let us assume a perturbation on the evolution of x,p(t), with a small amplitude parameter 0 < ε ≪1, ˙ x=F(x) + εp(t).(D.6) In order to keep arguing in terms of phases, one needs to find which is the impact of such perturbation on the phase of the oscillator. Applying again the chain rule to the time evolution of ϕ, and using Eq. (D.5) one obtains ˙ ϕ(x) = ∇xϕ·(F(x) + εp(t)) (D.7) =ω+ε∇xϕ·p(t) (D.8) Thus, as one would expect, the linear response of the phase ϕto a perturbation of the trajectory is given by its gradient. For this reason, the gradient of the phase along a certain unitary direction δx is usually referred to as the Phase Response Curve (PRC), Z(ϕ) := ∇xϕ·δx. Following the definition of gradient, given a infinitesimal perturbation εδx of a point on the limit cycle x, one can effectively compute the phase response curve as the normalized difference between the phase of x, and that of x+εδx, Z(ϕ) = ϕ(x)−ϕ(x+εδx) ε.(D.9) Thus, a positive (resp. negative) PRC indicates that the perturbation has the effect of advancing (resp. delaying) the phase of the oscillator. A zero PRC means that the perturbed point lies exactly on the isochrone of the original point. For simplicity, on the following we assume that the perturbation p(t) has direction δx and modulus p(t) so that ∇xϕ·p(t) = Z(ϕ)p(t). Thus, Eq (D.8) now reads ˙ ϕ(x) = ω+εZ(ϕ)p(t).(D.10) Systems with this form are known as Winfree-type models of phase oscillators,10 where the PRC, Z(ϕ), and the shape of the forcing p(t) are defined independently one from the other. D.2 From Winfree to Kuramoto Let’s assume now the phase dynamics of two interacting oscillators of Winfree type: ˙ ϕ1=ω1+εZ(ϕ1)g(ϕ2) (D.11) ˙ ϕ2=ω2+εZ(ϕ2)g(ϕ1).(D.12) Notice that now the forcing term p(t) has been replaced by a 2π-periodic interaction function g(ϕ). It is possible to further simplify this system by separating time scales on the dynamics of the phase model. Since ε≪1 in Eq. (D.11), one expects that the changes in the evolution of ϕjmostly happen at a slow time scale, since the fast dynamics is driven D PHASE DYNAMICS AND KURAMOTO MODEL 83 by the natural frequency of oscillation ωj. Let’s study then the time evolution of the slow variables φj=ϕj−ωjt. With this change of variables, equation (D.11) now reads ˙φj=εZ(φj−ωjt)g(φk−ωkt).(D.13) In order to retain the terms relevant for the evolution of φjone can expand the previous equation in Fourier series, ˙φj=εZ(φj−ωjt)g(φk−ωkt) = ε 4π2X m,n ˜ Z(m)˜p(n)e−i[m(φj−ωjt)+n(φk−ωkt)] (D.14) =ε 4π2X m,n ˜ Z(m)˜p(n)e−im(φj−φk)ei(mωj+nωk)t.(D.15) The last expression provides a separation of the fast and slow terms in the series: the terms with mωj+nωk≃0 (D.16) are resonant, and, therefore, lead to slow dynamics. The other terms, instead, have explicit dependence on t, and, therefore, lead to rapid changes. Thus, in order to keep only the relevant dynamics of the system, one can average out the non-resonant terms. Under the assumption that the natural frequencies of the two oscillator are comparable, i.e.,ωj≃ωk, the terms of Eq. (D.14) relevant for the dynamics are those with m=−n. This leads to ˙φj== ε 4π2 ∞ X n=−∞ ˜ Z(−n)˜p(n)e−in(φk−φj)=H(φj−φk) (D.17) where the coupling function Hcan be computed as H(ψ) := 1 4π2 ∞ X n=−∞ ˜ Z(−n)˜p(−n)e−inψ =1 2πZ2π 0 Z(t)p(t+ψ)dt. (D.18) Writing now system (D.11) into the original variables, we finally obtain: ˙ ϕ1=ω1+εH(ϕ2−ϕ1) (D.19) ˙ ϕ2=ω2+εH(ϕ1−ϕ2).(D.20) Phase oscillators ruled by equations of this type are known as Kuramoto-Daido type models. Despite their simplicity, such systems provide a good framework to study emerging phenomena on networks of oscillators. The simplest of these models corresponds to H being composed of a single harmonic, i.e., H(ψ) = sin(ψ−α), leading to the famous Kuramoto-Sakaguchi model. Nonetheless, the phase reduction of most nonlinear oscillators contain more harmonics in their coupling function, which usually also increase the complexity of the model dynamics. Here, we provided a simple intuitive arguments on how the time-scale separation is being performed. Formally, the coupling function Hemerges from applying averaging theory to system (D.11). These rigorous derivations also allow to extend this approach to more complex situations, including an arbitrary number of oscillators, cases in which the interaction function gdepends on both, pre and post-synaptic units, g(ϕ1, ϕ2), and cases with other resonant relation between the units. Nonetheless, here we want to emphasize the two main assumptions that allow to obtain Kuramoto-Daido models from arbitrary limit cycles: weak coupling (ε≪1), and weak frequency heterogeneities (ω1≈ω2). D PHASE DYNAMICS AND KURAMOTO MODEL 84 D.3 Phase reduction of a Stuart-Landau system In general, performing the phase reduction of a nonlinear oscillator requires a numerical approach. Due to its simplicity and symmetry properties, the Stuart-Landau oscillator is a notable exception. Here we reproduce a well-known result that a system of coupled Stuart-Landau oscillators leads to the Kuramoto-Sakaguchi model: ˙zi= (α+i ωi)zi−(γ+i β)|zi|2zi+ε N N X j=1 zj−zi.(D.21) Let’s first determine the phase variable of an uncoupled oscillator. We work using the polar representation: ˙r=α r −γ r3,(D.22) ˙ θ=ω−β r2.(D.23) The system contains an (unstable) fixed point at r= 0. For any other initial condition, r(0) = r0= 0 and θ(0) = θ0, the solution of the system reads: r(t) = sα γ+α r2 0−γe−2αt θ(t) = ω−αβ γ+θ0−β 2γln hγ αr2 0+e−2αt 1−γ αr2 0i. (D.24) As t→ ∞ the equation evolves towards a limit cycle with amplitude r∗=pα/γ and frequency Ω = ω−αβ γ. If the system is initialized exactly at the limit-cycle, then the phase is simply ϕ(r, θ) = θ. For initial conditions outside the invariant manifold, the system converges to the limit cycle with phase ϕ(r, θ) = θ−β 2γln γ αr2 0. From this explicit expression we can observe the isochrones of the Stuart-Landau oscillator are given by setting ϕ=k, where k∈[0,2π) is constant. We can also check that this definition of phase produces a uniform rotation from Eq. (D.5): ∇ϕ(r, θ)·F=ω−αβ γ. Let’s consider now the perturbed system: ˙zi= (α+i ωi)zi−(γ+i β)|zi|2zi+εδz. (D.25) As explained before, the PRC of the system can be computed as the phase gradient at a given perturbation direction, in this case Z(ϕ) = ∇zϕ·δz. For simplicity, we proceed in Cartesian coordinates: ϕ(x, y) = arctan y x−β 2γln γ α(x2+y2)(D.26) D PHASE DYNAMICS AND KURAMOTO MODEL 85 thus ∇zϕ=−1 x2+y2β γx+y,1 x2+y2x−β γy (D.27) =1 √α(−βcos(ϕ)−γsin(ϕ), γ cos(ϕ)−βsin(ϕ)) .(D.28) Therefore, the phase dynamics of a single, weakly perturbed Stuart-Landau oscillator are given by ˙ ϕ= Ω + ε √α[(−βcos(ϕ)−γsin(ϕ)) δx (γcos(ϕ)−βsin(ϕ)) δy] where δz = (δx, δy). In the coupled system, the perturbation comes from the other oscillators of the network, and thus chages in time. In particular δxi=1 N N X j=1 (xj−xi) = 1 Nrα γ N X j=1 (cos(ϕj)−cos(ϕi)) (D.29) δyi=1 N N X j=1 (yj−yi) = 1 Nrα γ N X j=1 (sin(ϕj)−sin(ϕi)).(D.30) Inserting these expressions into the PRC and simplifying the terms using trigonometric identities, we finally obtain ϕi= Ωi+εK β γ+εK N N X j=1 sin(ϕj−ϕi−α) (D.31) where Ωi=ωi−αβ γ, K =s1 + β2 γ2,and α= arctan β γ. In this case, no averaging is needed to obtain a Kuramoto-Daido type of model, since we have already obtained a coupling function that depends only on the phase differences. Moreover, the coupling is given by a single harmonic, thus ultimately providing a Kuramoto-Sakaguchi model. We emphasize, again, that this is specific for the StuartLandau oscillator, and other models can lead to coupling functions with several harmonics. E SYSTEM OF LINEAR COUPLED COMPLEX OSCILLATORS 86 E System of Linear Coupled Complex Oscillators We consider a network of Nlinearly coupled complex oscillators governed by ˙zi= (a+iω)zi+ N X j=1 cij zj, i = 1, . . . , N, (E.1) where a∈Ris the self–damping coefficient, ω∈Rthe intrinsic frequency, and C= [cij]∈ CN×Nthe coupling matrix. In vector form this reads ˙z=aI +iωI +Cz. If Cis diagonalizable with eigenpairs {(λk, vk)}N k=1, the change of variables z=V w decouples into modes ˙wk=a+λk+iωwk, k = 1, . . . , N. (E.2) Asymptotic stability (synchronization to z= 0) requires max kRea+λk(C)<0.(E.3) In particular, if Cis real symmetric with largest eigenvalue λmax, the simple sufficient condition is a < −λmax. More generally, the Master Stability Function (MSF) formalism233 characterizes the region in the complex plane α=a+λwhere perturbations decay. Early Lyapunov–matrix approaches appear in,234 and applications to small-world and complex topologies were developed in.235 Key Stability Criterion. max kRea+λk(C)<0⇐⇒ all trajectories decay to 0. Although the Master Stability Function (MSF) framework of Pecora & Carroll233 treats coupled dynamics in full generality, the purely linear network ˙zi= (a+iω)zi+X j cij zj has also been studied in its own right: •Large random systems. May analyzed the eigenvalue spectrum of aI +Cfor large random C, showing a sharp transition to instability when √N σ > |a|.236 •Amplitude death. Early stability analyses of coupled Stuart–Landau oscillators derive exactly the same linear condition against the trivial fixed point.237,238 •Hopf-normal-form variational equation. When linearizing two or more Hopf normal-form oscillators (Stuart–Landau) around the zero solution, one recovers the model above.239 In each case, asymptotic decay to zero reduces to max kRea+λk(C)<0. E SYSTEM OF LINEAR COUPLED COMPLEX OSCILLATORS 87 E.1 Constant Forcing and the Affine Shift Adding constant terms bidoes not alter the spectral–stability condition; it simply translates the asymptotic resting state from the origin to z∗=−A−1b. All transient decay (or growth) rates are exactly those of the original linear network. Consider appending a constant complex bias bito every node in (E.1): ˙zi= (a+iω)zi+bi,+ N X j=1 cij zji= 1, . . . , N. (E.4) Collecting the states z= (z1, . . . , zN)⊤and biases b= (b1, . . . , bN)⊤gives the affine system ˙z=Az +b, A ≡aI +iωI +C. Equilibrium. If Ais nonsingular—equivalently 0 /∈σ(A)—there is a unique fixed point z∗=−A−1b. (E.5) When Ais singular, a constant input may generate a line (or plane) of equilibria or even unbounded drift along ker Awhen b /∈Im A. Shift–of–origin reduction. Define the deviation w=z−z∗. Substituting (E.5) into (E.4) yields the homogeneous system ˙w=Aw, showing that constant forcing merely translates the flow; all eigenvalues, Jordan blocks, and Master-Stability-Function curves coincide with those of the original network. Hence the stability criterion max kRea+λk(C)<0 continues to be necessary and sufficient—now guaranteeing exponential convergence toward z∗instead of the origin. Dynamical consequences. •Stable spectrum (Re λk<0∀k): every trajectory decays at the same rates as before but settles on the static pattern (E.5). •Marginal spectrum (Re λk= 0 for some k): the constant drive excites those neutral modes, producing bounded oscillations (purely imaginary eigenvalues) or linear growth (λk= 0). •Unstable spectrum: divergence persists; bonly offsets the blow-up. E.2 From linear complex networks to Kuramoto phases Starting from the undamped linear network ˙zi(t) = i ωizi(t) + N X j=i Cijzj(t), i = 1, . . . , N, (E.6) E SYSTEM OF LINEAR COUPLED COMPLEX OSCILLATORS 88 we introduce polar coordinates zi(t) = ri(t)eiθi(t), ri(t)≥0, θi(t)∈R.(E.7) Substituting into (E.6), multiplying by e−iθi, and separating real and imaginary parts gives ˙ri= Re"X j=i Cij rjei(θj−θi)#,(E.8) ri˙ θi=ωiri+ Im"X j=i Cij rjei(θj−θi)#.(E.9) Writing Cij =Kijeiαij with Kij ≥0 yields ˙ θi=ωi+1 riX j=i Kij rjsin θj−θi+αij,(E.10) a Kuramoto–Sakaguchi–type equation with time–dependent amplitudes. To understand how a true phase model emerges, it is useful to rewrite (E.6) in vector form and add a uniform decay term, ˙ z= (A−γI)z, Aij =iωiδij +Cij,(E.11) as in the linear reformulations of Kuramoto.16,17 The dynamics are then determined by the eigenvalues λkand eigenvectors v(k)of A−γI. For a critical value of γ, one can arrange that Re λ1= 0,Re λk<0 (k≥2),(E.12) so that all but one mode decay. Decomposing the initial condition as z(0) = X k akv(k),(E.13) the solution is z(t) = X k akv(k)eλkt−−−→ t→∞ a1v(1)eλ1t,(E.14) provided a1= 0. Writing λ1=iΩcoll and v(1) = (v1, . . . , vN)⊤, we obtain zi(t)∼a1vieiΩcollt⇒ri(t)−−−→ t→∞ r∗ i:= |a1vi|, θi(t)≈Ωcollt+ arg(a1vi).(E.15) In the synchronized or collective–oscillation regime of the linear system, the amplitudes therefore do not decay to zero; rather, they converge to a fixed spatial profile r∗ i>0 determined by the dominant eigenvector. All other directions in state space are exponentially damped. E SYSTEM OF LINEAR COUPLED COMPLEX OSCILLATORS 89 Role of U(1) symmetry and collapse onto a reduced manifold. The network (E.11) is equivariant under global phase rotations, zi7→ eiϕzi∀i, ϕ ∈R/2πZ,(E.16) i.e. if z(t) is a solution then so is eiϕz(t). This defines a smooth action of the compact Lie group U(1) on the phase space CN, and the dynamics commute with that action. In general, such a continuous symmetry partitions the solution set into group orbits: each trajectory has a whole U(1) family of symmetry–related copies. Under standard uniqueness assumptions, this induces a conserved “group label” (a map from phase space to the group or its Lie algebra that is constant along trajectories), providing a direct link between continuous symmetries, conserved quantities and reduced manifolds.240 In local canonical coordinates, the conserved quantities appear as cyclic variables that do not enter the equations of motion except through their derivatives; trajectories then lie in a lower–dimensional manifold parameterized by these invariants.240 From the linear–systems side, the spectral picture above says that the state space R2N∼ = CNdecomposes into a one–complex–dimensional center eigenspace (spanned by v(1)) with Re λ1= 0 and a (2N−2)–dimensional stable subspace with Re λk<0. Center manifold theory then guarantees the existence of a locally invariant center manifold Wctangent to the center eigenspace at the origin, such that all nearby trajectories are exponentially attracted to Wcand the reduced dynamics on Wccapture the long–time behavior of the full system.241,242 In our case, Wcis generated by the center eigenvector and the U(1) action: up to small nonlinear corrections (e.g. saturating terms that fix the overall amplitude), the manifold is the set Wc≈av(1)eiϕ :a∈R+, ϕ ∈[0,2π).(E.17) Dissipation collapses all transverse directions onto this family, while the U(1) symmetry guarantees a neutrally stable direction along the global phase. If a weak nonlinearity or normalization fixes the amplitude a(or if we quotient out the trivial overall scaling), the effective attractor becomes a one–dimensional invariant circle S1generated by global phase rotations. This is the precise sense in which the combination of (i) U(1) symmetry and (ii) a spectral gap (all other eigenvalues with negative real part) leads to a collapse of the dynamics onto a reduced manifold, as emphasized in symmetry–based analyses of reduced manifolds in neural dynamics.240,243,244 Projecting onto this neutrally stable manifold and parametrizing the state by phases alone, the frozen amplitudes r∗ irender the couplings in (E.10) effectively constant, ˙ θi=ωi+X j=i ˜ Kij sin θj−θi+αij,˜ Kij =Kij r∗ j r∗ i ,(E.18) which is of Kuramoto–Sakaguchi type.16,17 The reduced phase model describes the dynamics on (or very close to) the U(1)–generated center manifold where all non–symmetric directions have been damped out. This construction should be viewed as one realization of Kuramoto dynamics, rather than an equivalence in the opposite direction. The effective couplings ˜ Kij in (E.18) are constrained by the spectrum of Aand the associated eigenvector profile r∗ i. In contrast, the Kuramoto model is typically introduced as a phenomenological phase reduction for H STUART–LANDAU OSCILLATOR: PARAMETERS AND GEOMETRY 96 Hence x(t)2+y(t)2=r2 ∞for all t, so the trajectory in the (x, y)-plane is a perfect circle of radius r∞. The parameter β appears only in the phase: ˙r=µ r −g r3,˙ θ=ω−β r2. Since βdoes not enter ˙r, the equilibrium amplitude r∞is independent of β. Therefore, β does not distort the circular shape into an ellipse. Instead, βshifts the frequency of rotation once the amplitude has saturated. Concretely, a nonzero βproduces an amplitudedependent frequency: as rgrows, ˙ θdecreases by β r2. On the limit cycle r=r∞, the constant frequency is ω−β(µ/g). Finally, the (x, y)-trajectory remains a uniform circular orbit at this shifted frequency; there is no ellipticity. In conclusion, while µsets the radius of the circle and ω(together with β) fixes its angular speed, only the pair (µ, g) determines the geometric shape (the radius) of the limit cycle. The parameter βinfluences when around the circle the oscillator moves (phase), but not how it traces out space (shape). H.2 Parameter Redundancy and Scaling in the Stuart–Landau Normal Form Must all four parameters appear, or can some be scaled away in the normal-form reduction to simplify analysis of the dynamics? H.2.1 Linear part: µand ω Near the Hopf bifurcation of a real dynamical system, one obtains a conjugate-pair of eigenvalues λ1,2(ν) = α(ν)±iΩ(ν), where νis the original system’s bifurcation parameter. By definition, at the bifurcation point ν=νc,α(νc) = 0 and Ω(νc)= 0. In the SL reduction µis chosen so that µ(νc) = 0 and µ≈α′(νc) (ν−νc) measures the distance from the Hopf point, and ωis (to leading order) the Hopf frequency Ω(νc). Because the fixed-time-units normal form must preserve the linear spectral center (i.e. the imaginary part at the bifurcation), both µand ωare in general essential parameters. However, one can make a rotating-frame transformation z(t) = u(t)ei ω t to eliminate ωentirely if one is only interested in autonomous amplitude dynamics. Under that change, ˙z=eiωt( ˙u+i ω u) yields an equation for uwith linear part (µ+i ω)u−i ω u = µ u. In other words, ˙u=µ u −(g+i β)|u|2u, so ωdrops out. Of course, if one cares about the absolute phase or wants to study phase interactions with external forcing, it may be convenient to keep ω. H.2.2 Nonlinear part: gand β The cubic coefficient in the SL normal form appears as the complex constant (g+i β). In principle, one can also nondimensionalize time and rescale zto remove one additional H STUART–LANDAU OSCILLATOR: PARAMETERS AND GEOMETRY 97 parameter. For example, define new units t′=g t, u(t′) = rg µ0 z(t), where µ0is some reference scale (e.g. µ0= 1). In these units, the equation becomes du dt′=˙z g=1 gh(µ+i ω)z−(g+i β)|z|2zi. Writing µ= ˜µ g and ω= ˜ω g, and using z=pµ0/g u, one finds du dt′=˜µ+i˜ωu−1 + i˜ β|u|2u, where ˜ β=β/g. In this scaled form: •The linear growth parameter becomes ˜µ=µ/g •The frequency becomes ˜ω=ω/g •The nonlinear amplitude coefficient is now unity. •The phase-amplitude coupling remains as the single ratio ˜ β=β/g. Thus, by an appropriate choice of time and amplitude scaling, one can reduce the SL normal form to ˙u= (˜µ+i˜ω)u−(1 + i˜ β)|u|2u, with only three essential real parameters ˜µ, ˜ω, ˜ β. In many analyses, using the transformation described above, one chooses the rotating frame to set ˜ω= 0, leaving ˙u= ˜µ u −(1 + i˜ β)|u|2u, with just two parameters ˜µand ˜ β. In that minimal form, ˜µcontrols the distance from bifurcation (radial growth), while ˜ βcontrols the amount of nonlinear frequency correction (phase–amplitude coupling). H.2.3 Summary: Which Parameters Are “Necessary”? •At the level of the full, original SL equation, one typically lists four real parameters µ, ω, g, β. •By scaling amplitude (i.e. set g= 1 in new units), the nonlinear damping coefficient is removed, leaving only the ratio β/g as the relevant nonlinear parameter. •By moving to a rotating frame, the linear frequency ωcan be subtracted off, so the form no longer explicitly contains ω. •Consequently, the minimal normal form that still captures amplitude growth and phase–amplitude coupling has only ˜µand ˜ β. •If one wishes to retain physical units (time in seconds, amplitude in volts, etc.), then ωmay remain as the observable oscillation frequency, and gremains the precise cubic coefficient. However, any analysis of scaling laws or bifurcation structure can be conducted in the reduced form with fewer parameters. H STUART–LANDAU OSCILLATOR: PARAMETERS AND GEOMETRY 98 Therefore, while the generic derivation of the SL normal form yields four parameters, two of them can be eliminated by choosing convenient time and amplitude scales (and, if desired, a rotating frame). The remaining parameters for the unfolded Hopf–SL dynamics are the real bifurcation parameter (distance ˜µ) and the dimensionless nonlinear frequencyshift ratio (˜ β). H.3 Alternative Form Emphasizing the Limit-Cycle Radius Here we cast the Stuart-Landau equation as a damped harmonic oscillator with dynamical damping and frequency. The equation ˙z= (µ+i ω)z−(g+i β)|z|2z can be rewritten to make explicit how |z|2is driven toward its steady-state value µ/g. Group the real and imaginary parts of the nonlinear coefficient: ˙z=µ−g|z|2z+iω−β|z|2z. Equivalently, ˙z=µ−g|z|2 | {z } real radial growth/dampingz+iω−β|z|2 | {z } instantaneous frequencyz or ˙z=µ−g r2 | {z } real radial growth/dampingz+iω−β r2 | {z } instantaneous frequencyz or ˙z=α(r) + iω(r)z reminiscent of the damped harmonic oscillator, (Equation 2.25) but with dynamical damping and frequency. The first term acts as feedback controller of the amplitude, a dynamical damping term keeping it close to the limit cycle radius, where a(r∞) = µ−g r2 ∞= 0. It can be seen as an oscillatory homeostatic mechanism. In this form: •The real part µ−g|z|2multiplies z. Writing z=r eiθ yields the radial equation ˙r=µ−g r2r, so that |z|=ris driven toward r∞=rµ g(µ > 0). Thus one sees directly that |z|2→µ/g as t→ ∞. •The imaginary part ω−β|z|2multiplies i z. In polar form this gives ˙ θ=ω−β r2, so the instantaneous oscillator frequency is shifted by β|z|2. On the limit cycle r=r∞=pµ/g, the asymptotic frequency is ω−β(µ/g). H STUART–LANDAU OSCILLATOR: PARAMETERS AND GEOMETRY 99 Because the term µ−g|z|2vanishes exactly when |z|2=µ/g, it is immediately clear from this form that |z|2must approach µ/g in order for the radial growth ˙rto vanish. Hence the limit-cycle amplitude emerges naturally as |z| → pµ/g. H.4 Stuart–Landau in Push/Pull Form Starting from ˙z= (µ+i ω)z−(g+i β)|z|2z, z =x+i y, r2=x2+y2, the Cartesian form is    ˙x=µ x −ω y −g r2x+β r2y, ˙y=µ y +ω x −g r2y−β r2x. Define the amplitude-dependent coefficient α, a dynamic damping term, and Ω, a dynamic instantaneous frequency, α(r) = −µ+g r2,Ω(r) = ω−β r2. Introduce the nonlinear differential operator L[f] := d dt +α(r)f=˙ f−(µ−g r2)f, where r2=x2+y2. Then the Stuart–Landau equations can be compactly written as L[x] = −Ω(r)y, L[y] = + Ω(r)x In this form, all growth (µ), saturation (g r2), and nonlinear phase effects (β r2) are absorbed into the operator Las a dynamic damping constant, while the right-hand side retains a “pure” rotation at instantaneous frequency Ω(r). When the limit cycle is reached, a= 0, we are back at the undamped harmonic oscillator. H.5 DC-Shifted Formulation of the Oscillator In modeling oscillatory neural dynamics, one often needs to account for a tonic bias or constant drive, or more generally constant or very slowly changing (compared to the natural timescale of the system) “forcing”. In fact, the signals generated by the model oftentimes need to be positive quantities (e.g., firing rates), or at least not centered at zero (membrane potential perturbations). What is the natural way to do this in the HO and SL cases? Harmonic Oscillator case. Consider the damped, unforced oscillator in complex form: ˙z= (a+iω)z, a ∈R, ω > 0. This equation give solutions centered at zero, which is undesirable if we want to relate them with firing rate models or membrane potentials. E.g., a constant electric field H STUART–LANDAU OSCILLATOR: PARAMETERS AND GEOMETRY 100 produces a DC shift in membrane potential or firing rate in more realistic models. How can this simple model accommodate it? A simple additive term provides the desired solution behavior (a DC shift). With a bias c, the equation becomes ˙z= (a+iω)z+c. (H.1) Introduce a constant shift C∈Cand define the change of variables (a translation) u=z+C. Then, the equation for uis ˙u= ˙z= (a+iω) (u−C) + c= (a+iω)u−(a+iω)C+c. The last term is just a complex constant and we can set it to zero with C=c a+iω so that ˙u= (a+iω)u— we recover the harmonic oscillator equation in the new coordinates. Hence, the generalized, shifted harmonic oscillator in Equation H.1 is a translation of the harmonic oscillator with a center at −C=−c/(a+iω) — see Figure H.1 (top left) for examples. The explicit solution is z(t) = z0+c a+iωe(a+iω)t−c a+iω. Hence, for a= 0, the motion is a uniform circular orbit around the displaced center −c/(iω); for a < 0, trajectories form logarithmic spirals converging to the point −c/(a+ iω). Stuart-Landau. Similarly to the HO case, a naive way to do so is to add a constant term c∈Cto the Stuart–Landau (SL) equation: ˙z= (µ+iω)z−(g+iβ)|z|2z+c. This DC–shifted SL now has •a limit cycle displaced away from the origin, •broken U(1) symmetry, •and modified amplitude/frequency balance. To recover a zero-mean description, one sets z=v+z0,where z0solves (µ+iω)z0−(g+iβ)|z0|2z0+c= 0. Substitution gives ˙v= (µ+iω)v−(g+iβ)v+z02(v+z0)+(µ+iω)z0+c−(g+iβ)|z0|2z0 | {z } =0 , so that the constant offset disappears. However, the new equation for vcontains extra quadratic and linear terms arising from the expansion of |v+z0|2(v+z0); it is no longer in the simple SL form. H STUART–LANDAU OSCILLATOR: PARAMETERS AND GEOMETRY 101 Despite this complication, normal form theory guarantees that any smooth perturbation of a Hopf normal form can be brought back—through a further smooth, near-identity change of variables—into the canonical Stuart–Landau structure, up to higher–order corrections. Le´on and Nakao (2023)245 provide expressions for the frequency and amplitude shifts up to second order in DC shift. More generally, v=u+H(u, ¯u), for a suitable cubic function H, eliminates all non-resonant quadratic and shifted cubic terms, leaving ˙u= (˜µ+i˜ω)u−(˜g+i˜ β)|u|2u to leading order. The new parameters ˜µ, ˜ω, ˜g, ˜ βabsorb the effects of the original bias c and shift z0. Thus, even with a DC offset, the local oscillatory dynamics near the Hopf bifurcation remain of Stuart–Landau type: a self-saturating limit cycle with phase–neutral drift or fixed point (depending on parameters). However, introducing a constant DC bias to the Stuart–Landau oscillator modifies its equilibrium position, breaks the symmetry, and shifts the critical Hopf bifurcation threshold, thus altering both amplitude and frequency of oscillations (see Figure H.1 for an example). Small biases result in second-order reductions in oscillation amplitude and frequency shifts, whereas sufficiently large DC offsets can completely eliminate the limit cycle. Analytically, this can be described as an imperfect Hopf bifurcation: the effective bifurcation parameter (µ) is renormalized by the DC term, requiring a higher original parameter value to sustain oscillations. Hence, persistent oscillations occur only when the system’s gain overcomes the bias-induced suppression; otherwise, the oscillator settles into a stable equilibrium without oscillations. H STUART–LANDAU OSCILLATOR: PARAMETERS AND GEOMETRY 102 Figure H.1: Top: Harmonic Oscillator (HO) and Stuart–Landau (SL) with a bias: phase portraits. Simulations display the effects of a small or larger DC bias offset, which induce frequency and amplitude changes, or destroy the limit cycle in SL while producing simple displacements in HO. We integrate ˙z= (µ+iω)z−(g+iβ)|z|2z+c with µ= 1, g= 1, ω= 2, β= 0, using RK4 (∆t= 0.01). Cases: (i) centered SL (c= 0), (ii) small DC shift (c= 0.5 + 0i) yielding a displaced limit cycle, (iii) large DC shift (c= 2 + 0i) leading to convergence to the fixed point. Initial condition for all runs: z0= 2R0with R0=pµ/g = 1. Trajectories shown after transients to emphasize the displaced cycle and the decay. Bottom: SL time-series overlay. The three cases are shown, highlighting changes in frequency and amplitude as well as extinction. I OSCILLATIONS, TOPOLOGY AND SIMPLICITY 103 I Oscillations, Topology and Simplicity The term oscillation is used across mathematics, engineering and neurophysiology, yet each community formalizes it (see Table I.2). In this section we show how all routes—dynamical systems, spectral tests, and Koopman eigenanalysis share the same core intuition: an oscillation is present when the data are better explained (or more concisely described) starting from the baseline of circular motion (S1topology) than by any simpler alternative. We will formalize this through an algorithmic definition. An intuitive classical definition is the following. Definition 1. A dynamical variable is said to oscillate when it exhibits sustained, approximately periodic departures around a reference value such that the system returns to a similar state after a characteristic interval T(its period), or, equivalently, at a dominant frequency f= 1/T. The repetition can be exact (strictly periodic) or approximate (quasi-periodic, weakly modulated, or noise-jittered); what matters is the recognisable cycle. This phrasing captures the intuition of “a signal that roughly repeats” while making explicit (i) the presence of a characteristic timescale and (ii) tolerance for imperfect cycles. This practical notion—a pattern that roughly repeats—is compatible both with modern spectrum-based detectors and with the dynamical-systems idea of an attracting limit cycle. But it can be generalized to include damped behavior: Definition 2. An oscillation is a self-sustained or externally driven fluctuation that revisits comparable states at quasi-regular intervals, identifiable either as a closed trajectory in state space or as a narrow-band spectral peak above the aperiodic background. This sentence bridges physics, signal processing, and neuroscience, preparing the ground for the algorithmic-information view developed in the following sections. We revise first the definitions in different sub-fields. Oscillations and the Koopman Operator. Alimit cycle is a closed orbit γof an autonomous ODE to which at least one neighbouring trajectory spirals (stable Floquet multipliers).249 The system’s intrinsic period Tis exact, and, as we explain, next, its Koopman operator possesses an eigenpair (ψ, λ =iω) whose argument θ= arg ψadvances uniformly, embedding the dynamics on the circle S1.250 Thus, periodic behavior can be precisely described from the Koopman operator perspective.250 Relaxation oscillators (e.g. FitzHugh–Nagumo) fit the same definition but spend long intervals on slow manifolds, producing nonsinusoidal waveforms with rich harmonic content. This approach fundamentally reframes the analysis. Instead of studying the nonlinear evolution of the system’s state vector x(t)∈Rn(governed by ˙ x=f(x)), the Koopman operator Ktdescribes the linear evolution of “observables” g(x), which are simply functions of the state. The operator is defined by how it advances any such function along the system’s trajectories: (Ktg)(x0) = g(x(t)), where x(t) is the flow starting from x0. The crucial insight is that Ktis a linear operator acting on the (infinite-dimensional) space of observables, even when the underlying dynamics f(x) are nonlinear. I OSCILLATIONS, TOPOLOGY AND SIMPLICITY 104 Table I.1: How major fields articulate the notion of oscillation. Field Typical wording Key nuance Classical physics & engineering “Repetitive or periodic variation, typically in time, of a quantity about a central value (often an equilibrium) or between two states.”246 Emphasizes small deviations from equilibrium and strictly periodic motion (e.g. undamped spring, AC current). Non-linear dynamics / mathematics Aperiodic orbit (or limit cycle) satisfying x(t+T) = x(t); at least one nearby trajectory spirals into it.1 Formal, coordinate-free; includes self-sustained oscillators (Van der Pol, Hodgkin–Huxley) and supports stability analysis. Neuroscience “Rhythmic or repetitive neural activity observed at all levels of the CNS.”247 Focus on multi-scale biological generators; amplitude mainly indexes population synchrony, not a single source. Signal processing / spectral view Oscillation = “narrow-band peak that rises above the aperiodic 1/fbackground in the power spectrum.”248 Detects oscillations without explicit time-domain periodicity; robust for noisy or burst-like data. Because Ktis linear, we can use spectral methods. For a limit cycle, the operator’s infinitesimal generator L(where Kt=eLt) possesses an eigenpair (ψ, λ =iω) corresponding to the system’s fundamental frequency ω= 2π/T. This special observable ψ, the Koopman eigenfunction, evolves simply in time: ψ(x(t)) = (Ktψ)(x0) = eλtψ(x0) = eiωtψ(x0). Consequently, its argument θ= arg ψadvances uniformly (θ(t) = θ0+ωt), embedding the complex, multi-dimensional dynamics onto a simple rotation on the circle S1. Relaxation oscillators (e.g. FitzHugh–Nagumo) fit the same definition, but their corresponding eigenfunctions ψare more complex (capturing all the harmonics), resulting in nonsinusoidal waveforms. It is important to clarify what the eigenfunction ψrepresents. It is a special observable, meaning it is a function of the state vector, ψ(x), that maps the n-dimensional state space (where the dynamics are nonlinear) to the complex plane C(where the dynamics are linear). Its special property is that when evaluated along a trajectory x(t), its value evolves with perfect simplicity according to ψ(x(t)) = eiωtψ(x0). Crucially, one must distinguish the function ψ(x)—which is a static, complex-valued map on the state space—from its evolution in time,ψ(x(t)). The function ψ(x) itself is generally not periodic. In essence, ψacts as a ”magic” coordinate transformation. While the state x(t) traces a complex orbit, the scalar observable ψ(x(t)) simply rotates in the complex plane at a constant frequency ω. The level sets of its phase, θ= arg ψ(x), are the system’s isochrons: surfaces of points in the state space that all share the same asymptotic phase on the limit cycle. I OSCILLATIONS, TOPOLOGY AND SIMPLICITY 105 Koopman perspective as compression. Finding an eigenfunction ψwith generator eigenvalue λ=iω is itself a compression. Once this function ψis known, the full n-dimensional trajectory x(t) can be encoded (losslessly near the attractor) by a single complex phase variable ψ(x(t)), or separated into its phase θ(t) = arg ψ(x(t)) and amplitude r(t) = |ψ(x(t))|.250 In effect, the Koopman transform finds a “magic” coordinate system ψin which the nonlinear dynamics become simple linear rotation. It replaces a complex waveform with uniform rotation on S1, turning geometry into a one-line program: output r(t)·cos(θ(t)). Connection to Diffeomorphisms and Lie Groups. The flow of the dynamical system, Φt:M→M, maps an initial state x0to its position at time t,x(t)=Φt(x0). For a smooth vector field f, this flow Φtis a diffeomorphism (a smooth, invertible map with a smooth inverse). The set of all such flows {Φt}t∈Rforms a one-parameter group under composition: Φt◦Φs= Φt+s. This group is a subgroup of the full infinite-dimensional Lie group of all diffeomorphisms of the state space M, denoted Diff(M). The Koopman operator Ktis the pullback (or composition) operator induced by this flow. It is defined as (Ktg) = g◦Φt, or (Ktg)(x) = g(Φt(x)). The Koopman operator family {Kt}also forms a one-parameter group, (Kt◦ Ks)g=g◦Φs◦Φt=g◦Φs+t=Kt+sg. Crucially, this provides a linear representation of the (nonlinear) flow group {Φt}on the (infinite-dimensional, linear) space of observable functions. This relationship extends to their generators. The generator of the flow group {Φt}is the vector field fitself, which is an element of the Lie algebra Vect(M) (the space of vector fields on M). The generator of the Koopman group {Kt}is the operator L. This generator Lis precisely the Lie derivative with respect to the vector field f,L=Lf, which acts on observables gas Lg = (f·∇)g. Thus, the Koopman framework ”lifts” the nonlinear dynamics from the state space Mto a linear representation on a function space, where the generator is the Lie derivative. Signal-processing criteria . In experimental neurophysiology, oscillations are usually detected rather than proven. Two widely used operational rules are (i) the BOSC power + duration test251 and (ii) spectral parameterisation (“FOOOF”) that labels any narrow-band bump lying above the aperiodic 1/fbackground as oscillatory.248 Both are statistical surrogates for asking whether a periodic template explains the data substantially better than a broadband model. Spectral peaks are suggestive of reduced entropy or increased compressibility. Table I.2 aligns the main modeling traditions—from classical limit–cycle theory to modern information–theoretic views—while the text links them through the common intuition that circular motion in an abstract coordinate revealed by compression. Next, we discuss the formalization of this intuition and generalization of the above definitions using the language of algorithmic information theory and compression (AIT).252 I.1 Algorithmic-information definition Algorithmic Information Theory (AIT) takes a computational perspective and quantifies the information content of individual objects via computation rather than probability. For a binary string x, the (prefix) Kolmogorov complexity KU(x) is defined as the length J LINEAR OPERATORS, GREEN–LAPLACE TOOLS, AND E-I OSCILLATIONS: A PEDAGOGICAL VIEW 112 with α=a 2mand ωd=√4mb−a2 2m. For fixed a, b > 0 and m→ ∞,α→0, ωd∼pb/m, and arctan(ωd/α)→π/2, hence tpeak ∼π 2rm b→ ∞. Inertia therefore increases the causal delay. Importantly, the overdamped log formula must not be used once the poles become complex; the underdamped expression governs the peak. J.3 E–I motifs, Barkhausen conditions, and where the phase lag comes from A pedagogical route to oscillation is to rewrite the undamped oscillator as a pair of coupled first-order filters. Let z=x+iy and consider ˙z= (a+iω)z, a ≥0, ω > 0. Separating real and imaginary parts gives ˙x=ax −ωy, ˙y=ay +ωx. Each equation is a leaky integrator driven by the other in 90◦phase. When a= 0 the loop produces sustained oscillations; when a > 0 the envelope decays as e−at. This “push–pull” view is the simplest template to keep in mind as we turn to E-I populations. Barkhausen as a phase-gain budget. For a feedback loop with transfer L(s), necessary conditions for linear self-oscillation at ω0are |L(jω0)|= 1 and arg L(jω0)=0◦ (mod 360◦)vonWangenheim10. The phase condition pins ω0by balancing element lags; the magnitude condition pins the product of gains, including the slope κof the nonlinearity at the bias point. Nonlinear saturation then stabilizes amplitude ˚ Astr¨om&Murray08. Wilson–Cowan with first-order synapses. Linearize about a fixed point with slope κand unit time constants: ˙x=−(1 −weeκ)x−weiκ y +Ix,˙y=−y+wieκ x −wiiκ y. The E→I→E loop has TEI (s) = weiwieκ2 (s+ 1)2, which can supply up to −180◦of phase—not enough alone to close the loop at a finite ω. Excitatory self-coupling adds TEE(s) = weeκ s+ 1, so the characteristic equation 1 −TEE(s)−TEI(s) = 0 can satisfy Barkhausen at some ω0. A bias Ixensuring κ > 0 is essential Wilson&Cowan72. This matches the intuition that two first-order elements need additional phase (or an explicit transmission delay) to reach 360◦. J LINEAR OPERATORS, GREEN–LAPLACE TOOLS, AND E-I OSCILLATIONS: A PEDAGOGICAL VIEW 113 Jansen–Rit with second-order synapses. For the cortical column with second-order synapses, ¨x+ 2a˙x+a2x=Aa κ(−w y +Ix),¨y+ 2b˙y+b2y=Bb κ(w x), the linearized synapses are band-pass filters Hexc(s) = Aa κ (s+a)2, Hinh(s) = Bb κ (s+b)2. The E→I→E loop transfer T(s) = −w2Hexc(s)Hinh(s) can contribute −360◦of phase on its own, so self-excitation is not required to meet the phase condition. A bias Ixmaintaining κ > 0 sets |T(jω0)|= 1 at the selected frequency Jansen&Rit95. Empirically, ω0follows the loop’s effective delay, which is controlled by synaptic poles a, b and any axonal conduction delay. Information flow and effective loop delay. For narrowband loops, each element’s phase behaves as φk(ω)≈ −ω τknear resonance, so Pkφk(ω0) = 0◦implies ω0≈ 2πn/τloop, where τloop ≈ −Pkφk(ω0)/ω0. Second-order synapses provide larger group delay around their passband than single-pole synapses. This is one reason gamma-range E-I oscillations arise robustly once the synaptic dynamics are at least second order Buzs´aki&Wang12. First-order worked example (frequency response to a sinusoid). To fix ideas, solve ˙x+ax =ejωt with a > 0. In Laplace, (s+a)X(s) = 1 s−jω ⇒X(s) = 1 (s+a)(s−jω). Partial fractions and inversion give x(t) = e−at −(a+jω)+ejωt a+jω. The transient term dies as t→ ∞; the steady state is x(t) = H(jω)ejωt with H(jω) = 1/(a+jω). Hence |H(jω)|= 1/√a2+ω2and arg H(jω) = −arctan(ω/a), and τg(ω) = a/(a2+ω2). This calculation makes explicit how a pole at −asets both decay and phase lag. paragraphTime-domain worked example (convolution with an alpha kernel). Consider h(t) = Aate−atH(t) and an input σ(t). Then y(t) = (h∗σ)(t) and L{h}=Aa/(s+a)2. Multiplication in s-space becomes Y(s) = Aa (s+a)2Σ(s)⇐⇒ (∂t+a)2y(t) = Aa σ(t), realizing the second-order operator directly from the kernel. This is the synaptic analog of driving a critically damped mass-spring-damper. J LINEAR OPERATORS, GREEN–LAPLACE TOOLS, AND E-I OSCILLATIONS: A PEDAGOGICAL VIEW 114 Pointers for further reading A rigorous, operator-theoretic treatment of L[h] = δand the jump conditions appears in Stakgold&Holst11. For feedback, phase, and the Barkhausen criterion in context see ˚ Astr¨om&Murray08 and the clarification in vonWangenheim10. For neural mass modeling with firstand second-order synapses see Wilson&Cowan72 and Jansen&Rit95; for conductance-based synaptic kinetics and canonical PSP shapes see Destexhe et al.94; and for a broader review of E–I mechanisms of cortical rhythms see Buzs´aki&Wang12. References Wilson, H. R., & Cowan, J. D. (1972). Biophysical Journal, 12, 1–24. 10.1016/S00063495(72)86068-5. Jansen, B. H., & Rit, V. G. (1995). Biological Cybernetics, 73, 357–366. 10.1007/BF00199471. Ermentrout, G. B., & Terman, D. H. (2010). Mathematical Foundations of Neuroscience. Springer. 10.1007/978-0-387-87708-2. Stakgold, I., & Holst, M. (2011). Green’s Functions and Boundary Value Problems (3rd ed.). Wiley. 10.1002/9780470906538. von Wangenheim, L. (2010). On the Barkhausen and Nyquist stability criteria. Analog Integrated Circuits and Signal Processing.10.1007/s10470-010-9506-4. ˚ Astr¨om, K. J., & Murray, R. M. (2008). Feedback Systems. Princeton Univ. Press. Book page. Destexhe, A., Mainen, Z. F., & Sejnowski, T. J. (1994). Journal of Computational Neuroscience, 1, 195–230. 10.1007/BF00961441. Buzs´aki, G., & Wang, X. -J. (2012). Neuron, 72, 203–229. 10.1016/j.neuron.2012.09.041.