Neural mass and field models generate brain rhythms using the notion of population firing rates, and in some sense can be thought of as sitting above more detailed networks of interacting conductance-based spiking neuron models. The latter are hard to analyze in the raw, given that they are both high-dimensional and nonlinear, and are typically studied with tools from computational neuroscience, as exemplified by the work of the Human Brain Project [1]. In contrast neural mass and field models are much more amenable to mathematical analysis, as reviewed in Coombes et al. [2] and Cook et al. [3], and although these low-dimensional models are typically not derived from any underlying microscopic spiking dynamics they can be motivated by a number of phenomenological arguments for the evolution of coarse-grained neuronal variables. Nonetheless, they are only expected to provide appropriate levels of description for many thousands of near identical interconnected neurons with a preference to operate coherently. As such, they are not ideally suited to studying phenomenon that are known to be associated with changes of synchrony, such as the post stimulus response ubiquitously seen in human neuroimaging studies, see e.g., Mullinger et al. [4]. However, a new type of neural mass model has recently been developed that is able to capture the phenomenon of event-related synchronization/desynchronization (ERS/ERD) that is believed to underlie the changes in power seen in brain spectrograms, or more specifically for that seen in magneto-encephalography (MEG) data showing post-movement beta rebound [5]. Importantly, this new mean-field model is an exact description of a globally coupled network of heterogeneous θ-neuron models in the thermodynamic limit. The θ-neuron, or Ermentrout–Kopell canonical model, is the normal form for the saddle-node on a limit cycle bifurcation, and for constant stimulation can generate low firing rates typical of those seen in real cortical neurons [6]. Interestingly, the resulting mean-field model has a population firing rate that depends on the degree of population synchrony, and has a richer dynamical repertoire than standard neural mass models. Moreover, it can incorporate realistic models of synaptic currents (both chemical and electrical) and the mass model is easily incorporated into network studies (both discrete and continuous). In this perspective article I review some of the main features of this new neural population model, its application to neuroimaging studies, extensions of the model to incorporate further biologically important components, and some mathematical challenges for developing similar mean-field reductions for other spiking models.

The θ-neuron describes a single neuron using a phase θ ∈ [0, 2π) such that a spike is generated whenever θ passes through π from below. For a stimulus I it evolves according to the ordinary differential equation (1)θ.=1-cosθ+(1+cosθ)I. The θ-neuron is formally equivalent to a quadratic integrate-and-fire (QIF) model for voltage dynamics [7] under the transformation v = tan(θ/2), and here we note that this can also be written as a Möbius transformation v = iM(e^iθ; −1, 1, 1, 1), where (2)M(z;a,b,c,d)=az+bcz+d,  a,b,c,d∈ℂ, ad-bc≠0. A network of θ-neurons can be described with the introduction of an index n = 1, …, N and the replacement I → μ_n + I_n, where I_n describes the synaptic input current to neuron n and μ_n is a constant drive. For a globally coupled network, the current arising from a chemical synapse can be written in the form I_n = g(t)(v_syn − v_n) for some global conductance g, synaptic reversal potential v_syn, and local voltage v_n. In the voltage representation this leads to the network equations (3)v.n=μn+vn2+g(vsyn-vn),  n=1,…,N, Where the voltage is reset to v_r → −∞ whenever v_th → ∞. Here, μ_n is a random variable drawn from a Lorentzian distribution with center μ₀ and width at half maximum Δ. For a treatment of the model with a Gaussian and q-Gaussian distribution of drives, see Klinshov et al. [8] and Pyragas and Pyragas [9], respectivley. We will work in the thermodynamic limit (N → ∞), and choose a model of conductance change that is driven by delta–Dirac spikes in the form: (4)Qg=kR,  R(t)=limN→∞1N∑n=1N∑m∈ℤδ(t-Tnm), Where k is a strength of coupling, R is the population firing rate, and Tnm denotes the mth time that neuron n spikes. Here Q is a linear differential operator whose Green's function is chosen to capture the post-synaptic response to the arrival of a spike. A common choice is the α-function η(t) = α²te^−αtH(t), where H is a Heaviside step function. The Q for this choice is the second order operator Q = (1 + α⁻¹d/dt)². For a schematic of the spiking network described below, see Figure 1.

To date two independent approaches have been developed for obtaining mean-field equations for the spiking network model described above (with N → ∞), both finding exact solutions to the continuity equation describing the distribution of states. The first approach, due to Luke et al. [10], considers the θ-neuron representation and makes use of the Ott–Antonsen (OA) ansatz [11]. Here a Fourier series representation for the distribution of phases with all Fourier coefficients restricted to be powers of a function a(μ, t) ∈ ℂ is shown to be an exact solution provided that a(μ, t) satisfies a certain differential equation. Moreover, the average of a(μ, t) over the distribution of drives, denoted 〈a(μ, t)〉_μ, has a physical interpretation as the Kuramoto order parameter for synchrony Z(t)=limN→∞N-1∑n=1Neiθn(t). Recent work of Cestnik and Pikovsky generalizes this approach to describe relaxation to the stable OA manifold in an exact manner [12]. The second approach due to Montbrió et al. [13], considers the QIF representation and shows that a Lorentzian choice for the distribution of voltages, with center y(μ, t) and width at half maximum x(μ, t), is also a solution provided that the pair (x, y) obey a set of coupled differential equations. Moreover, after averaging over the distribution of drives, this pair can also be related to physical variables. Introducing the average voltage V(t)=limN→∞N-1∑n=1Nvn(t) then (R(t),V(t))=(π-1〈x(μ,t)〉μ,〈y(μ,t)〉μ). Thus, in both approaches an evolution equation for physically meaningful macroscopic variables is generated. Introducing the complex number W(t) = πR(t) + iV(t) then this evolution equation from the approach of Montbrió et al. can be written succinctly as the complex Ricatti equation (5)W.=Δ-gW+i[μ0+gvsyn-W2]. It is noteworthy that the distribution of phases arising from the OA approach can be written as a Lorentzian distribution in e^iθ. Thus, exploiting a result of McCullagh [14] that if the random variable X has a Lorentzian distribution with complex parameter ζ, then the random variable Y = M(X; a, b, c, d) has a Lorentzian distribution with parameter M(ζ; a, b, c, d), and remembering that the phase and voltage representations are related by a Möbius transformation, it can be shown that Z(t) and W(t) are related by the conformal map Z = M⁻¹(W^*; −1, 1, 1, 1), or more specifically: (6)Z(t)=1-W*(t)1+W*(t). Thus, the dynamical equation for Z(t) can either be obtained directly following the approach of Luke et al. or by applying (Equations 6–5). Either way, this gives the complex Ricatti equation for the evolution of synchrony as (7)Z.=-i(Z-1)22+(Z+1)22[-Δ+iμ0+ivsyng]-(Z2-1)2g. A method for analyzing periodic solutions to complex Ricatti equations of the form given by Equations (5), (7) (in the self-consistent case that g(t) is periodic) has recently been proposed in Omel'chenko [15], exploiting the fact that every periodic solution corresponds to a fixed point of a Möbius transformation. For a numerical bifurcation analysis of Equation (7), see Coombes and Byrne [16] where the basic mechanism for generating oscillations at the single population level is shown to be a Hopf bifurcation (though note that oscillations can also emerge out of the blue along a branch of isolas in models with two interacting populations). The mean-field model can also exhibit canards and bursting solutions under slow periodic forcing [17], as well as chaos for moderately fast periodic forcing [18].

In contrast to more traditional phenomenological neural mass models, as exemplified by that of Wilson and Cowan [19], the population firing rate is a derived quantity that takes the form R(t) = π⁻¹Re W(t) = π⁻¹Re ((1 − Z^*)/(1 + Z^*)). Thus, the firing rate is determined by the degree of spike synchrony within a population (measured by the Kuramoto order parameter), rather than just keeping track of the fraction of active neurons in a given interval of time (as in the original heuristic derivation of the Wilson–Cowan equations). In cases where spike synchrony is not important, then the Wilson-Cowan equations would still be expected do a good job of qualitatively describing network dynamics. The main difference between the heuristic Wilson–Cowan model and next generation neural mass model is that the latter explicitly includes a dynamic component for describing spike synchrony within a population. Moreover, the next generation neural mass equations are closed in slightly different way in that R only depends indirectly on g through the dynamical evolution of W or Z [20]. Thus, we may select from one of the two equivalent perspectives depending on context. For example, when modeling MEG studies that highlight ERS/ERD it is convenient to use the Z representation for the Kuramoto order parameter. In contrast, for studies involving voltage sensitive dyes it is more natural to use the W representation and extract the average voltage from V = Im W. For a detailed derivation of the mean-field equations described above, see Coombes and Wedgwood [21, Chapter 8].

With the inclusion of an external time-dependent input, the mean-field equations described above have previously been shown to generate behavior consistent with neuroimaging experiments exhibiting so-called beta rebound [5]. Modulations in the beta power of human brain rhythms are known to occur during and after movement. A sharp decrease in neural oscillatory power is observed during movement followed by an increase above baseline on movement cessation. Both are captured in the response of the mean field model to simple square wave forcing, as illustrated in Figure 2.

The extension of the mean-field equations to describe interacting excitatory and inhibitory sub-populations, with both reciprocal and self connections, is straightforward [16, 22–26], as is the further extension to whole brain models utilizing (connectivity and delay) data from the Human Connectome Project [27], and the corresponding limit of spatially continuous neural fields [28, 29]. Moreover, for certain forms of noise in the microscopic equations the mean-field reduction can also be generalized [30–36]. A variation on the Lorentzian ansatz (in the voltage representation) combined with a moment closure approach has also been shown to be useful for approximating networks of Izhikevich neurons [37, 38], namely QIF neurons with a form of adaptation [39].

The mean-field reduction remains viable with the inclusion of electrical synapses. These are electrically conductive links between two adjacent nerve cells at a so-called gap junction (modeled as a passive ohmic resistance). The synaptic delay for a chemical synapse is typically in the range 1 − 100 ms, while the synaptic delay for an electrical synapse may be only about 0.2 ms. Little is known about the functional aspects of gap junctions, but they are thought to be involved in the synchronization of neurons, see e.g., Bennet and Zukin [40]. For their inclusion in the next generation framework, see Laing [41], Pietras et al. [42], Montbrió and Pazó [43], Byrne et al. [44], and Clusella et al. [45]. Recent work has shown that they can have a strong effect on network excitability and that gap junction coupling strength can play a key role in generating empirically-observed patterns of large scale spatio-temporal brain activity [46].

In recent years there has been a growing use of the next generation population modeling framework in a variety of different neuroscience contexts, including abnormal beta rebound in schizophrenia [47], beta bursts seen in single trial data [48], the generation of gamma oscillations [49] and theta-nested gamma oscillations [50], syllable segregation in speech comprehension [51], the generation of functional connectivity in whole brain networks [48, 52], seizure propagation [53], cross-frequnecy coupling [54] and communication through coherence [55], working memory (with the inclusion of short-term synaptic plasticity) [56], non-invasive transcranial brain stimulation [57], and neurodegeneration [58].

The mean-field reduction discussed here leads to low dimensional equations that are in the spirit of previous neural mass and field equations, albeit with a direct link to an underlying microscopic spiking dynamics [59]. They have proven to have very rich dynamics, at the node [13, 16, 60], network [46, 48, 52], and continuum level [28, 29, 41, 44]. Nonetheless, their pattern forming properties remain relatively unexplored, providing an open challenge to the applied mathematics community. However, it is worth bearing in mind the caveats of heterogeneity, global coupling, and the thermodynamic limit in the reduction. For identical neurons one could instead turn to the Watanabe–Strogatz ansatz (instead of OA) as considered by Laing [61]. Moving away from global coupling is difficult without some form of approximation, and here the techniques developed by Thibeault et al. [62] might be useful for modular networks, and those of Restrepo and Ott for assortative networks [63] as considered in Chandra et al. [64] and Laing and Bläsche [65]. Moving away from the thermodynamic limit, to address finite size effects, may be possible with a perturbative approach (around the asynchronous state), such as the Bogolyubov–Born–Green–Kirkwood–Yvon hierarchy with an appropriate moment-closure approximation [66], or a path-integral formalism to determine the evolution of system covariances [67]. Interestingly, it has recently been shown that finite size effects can be captured with the inclusion of stochastic forcing by shot noise into the mean-field equations [68].

In essence the mean-field reduction discussed here is well suited to describing systems which dynamically evolve between an incoherent state and a partially synchronized state, which is often the case in phase oscillator networks with interactions that are prescribed by sinusoidal functions (and the large global network generates a unimodal distribution of phases). In systems with higher harmonics in their interactions, it is possible for cluster states to emerge, and the OA ansatz breaks down (as it cannot describe distributions with more than one peak). Thus, there is an interest in either generalizing the OA approach or developing an alternative. As a generalization it might be worth considering whether there are other models like the QIF that can be described by a multi-peak circular distribution [69] of network phases ρ(θ) (after integrating over the distribution of drives) of the form (8)ρ(θ)∝∏m=1M1-|Zm|2|eiθ-Zm|2, Where Zm(t)=limN→∞N-1∑n=1Neimθn(t) are the Kuramoto–Daido order parameters. This choice effectively recovers the OA ansatz when M = 1. However, this would again only cover a special case of a particular choice of idealized single neuron model. Xiao et al. [70] have recently developed a data-informed mean-field approach that couples to the average network voltage (and note that the next generation approach also has a coupling to the average voltage), and it may be that this is a more promising line of attack for describing networks of biophysically realistic neurons.

Finally, it is important from a modeling perspective to realize that the brain consists of more than just the cortex and that for many sub-cortical structures intrinsic nonlinear ionic currents can dominate the firing rate response. It remains an open challenge to develop mean-field models that can incorporate such important features, though for slow intrinsic currents some progress has been possible, see e.g., Modhara et al. [71] for a mean-field model of the thalamus (sensory gateway to the cortex). Similarly, it is important to recognize that there is more to brain tissue than neurons, and that for many clinical applications, some form of coupling to the extracellular space is beneficial; a case in point being epilepsy, and the important role that extracellular K⁺ has in the genesis and propagation of seizure states [72].

The original contributions presented in the study are included in the article/supplementary material, further inquiries can be directed to the corresponding author.

SC is the single author of this article and contributed to the article and approved the submitted version.