Before and beyond the Wilson-Cowan equations

TL;DR

This paper explores the origins, biological basis, and extensions of the Wilson-Cowan equations, a fundamental model in computational neuroscience, providing historical context and discussing how to go beyond this classical framework.

Contribution

It offers a detailed historical and theoretical analysis of the Wilson-Cowan equations and proposes directions for extending their applicability in neural modeling.

Findings

Clarifies the biological and theoretical foundations of the Wilson-Cowan equations

Discusses limitations of the classical equations and potential extensions

Provides context for the equations' role in neural dynamics modeling

Abstract

The Wilson-Cowan equations represent a landmark in the history of computational neuroscience. Among the insights Wilson and Cowan offered for neuroscience, they crystallized an approach to modeling neural dynamics and brain function. Although their iconic equations are used in various guises today, the ideas that led to their formulation and the relationship to other approaches are not well known. Here, we give a little context to some of the biological and theoretical concepts that lead to the Wilson-Cowan equations and discuss how to extend beyond them.

Equations

$$\tau \dot{a_i}$$

$$\tau \dot{a}(x)$$

$$\tau \dot{s_i}$$

$$\dot{s_i}$$

$$=-s_i+\sum _j{w}_{ij}{f}_j\left(s_j\right)+I_i$$

$$a_i=f_i(s_i)$$

$$f_1(Q)=\sum _{n,p,h}P(Q;n,p,h)P_{nph}[I_0(Q)]I_0(Q)={\int }_{\Omega }O(Q,X)f_0(X)dX$$

$$F(t+\tau )=R{P}_{(q-1)}Fk$$

$$P_{(q-1)}={\mathrm{e}}^{-\overline{N}}\frac{{\overline{N}}^{(q-1)}}{(q-1)!}$$

$$\overline{N}(x,t)={\int }_{-\infty }^0{\int }_{-\infty }^{\infty }w(x-x^{\prime })F(x^{\prime },t^{\prime })\alpha (t-t^{\prime })\mathrm{d}{x}^{\prime }\mathrm{d}{t}^{\prime }$$

$$H\psi =\kappa F$$

$$F(x,t)$$

$$F\approx f(\psi )+A_1\frac{\partial \psi }{\partial t}+A_2\frac{{\partial }^2\psi }{\partial t^2},$$

$$H=a+b\frac{\partial }{\partial t}+c\frac{{\partial }^2}{\partial t^2}+{\nabla }^2.$$

$$\psi (\mathbf{r},t)=\int F\left(\mathbf{R},t-p/v\right)w(|\mathbf{r}-\mathbf{R}|)d\mathbf{R},$$

$${\nabla }^2\psi =$$

$$+\left(\frac{1}{v^2}-4\pi A_2\right)\frac{{\partial }^2\psi }{\partial t^2}-4\pi f(\psi )$$

$${\nabla }^2\psi =\mu \psi +\dot{\psi }-f(\psi )$$

$$\dot{\psi }=-\mu \psi +f(\psi )$$

$$\left[1-{\int }_{t-r}^t{A}_i\left(t^{\prime }\right)d{t}^{\prime }\right]f_i(z_i)\delta t,$$

$$z_i\equiv {\int }_{\infty }^t\alpha \left(t-t^{\prime }\right)\left[W_{iE}{A}_E\left(t^{\prime }\right)-W_{iI}{A}_I\left(t^{\prime }\right)+Q_i\left(t^{\prime }\right)\right]d{t}^{\prime },$$

$$f(x)={\int }_0^{x(t)}p(\theta )\mathrm{d}\theta .$$

$$A_E(t+\tau )=$$

$$\displaystyle\times f_{e}\Bigg{(}\int_{-\infty}^{t}\alpha(t-t^{\prime})\Big{[}W_{EE}A_{E}(t^{\prime})-W_{EI}A_{I}(t^{\prime})+Q_{E}(t^{\prime})\Big{]}dt^{\prime}\Bigg{)}$$

$$A_I(t+\tau )=$$

$$\displaystyle\times f_{I}\Bigg{(}\int_{-\infty}^{t}\alpha(t-t^{\prime})\Big{[}W_{IE}A_{E}(t^{\prime})-W_{II}A_{I}(t^{\prime})+Q_{I}(t^{\prime})\Big{]}dt^{\prime}\Bigg{)}$$

$$a_i(t)=\frac{1}{s}{\int }_{t-s}^t{A}_i\left(t^{\prime }\right)\mathrm{d}{t}^{\prime }$$

$${\int }_{t-r}^{t}A\left(t^{\prime }\right)\mathrm{d}{t}^{\prime }\to ra(t)$$

$${\int }_{-\infty }^t\alpha \left(t-t^{\prime }\right)A\left(t^{\prime }\right)\mathrm{d}{t}^{\prime }\to ka(t)$$

$${\tau }_E\frac{\mathrm{d}{a}_E}{\mathrm{d}t}$$

$${\tau }_I\frac{\mathrm{d}{a}_I}{\mathrm{d}t}$$

$$\tau \frac{\partial A_E(\mathbf{x},t)}{\partial t}=$$

$$\times f_E[{\int }_{\Omega }d{\mathbf{x}}^{\prime }{\rho }_E({\mathbf{x}}^{\prime })W_{EE}\left(\mathbf{x}-{\mathbf{x}}^{\prime }\right)A_E\left({\mathbf{x}}^{\prime },t\right)$$

$$-{\int }_{\Omega }d{\mathbf{x}}^{\prime }{\rho }_I({\mathbf{x}}^{\prime })d{\mathbf{x}}^{\prime }{W}_{EI}\left(\mathbf{x}-{\mathbf{x}}^{\prime }\right)A_I\left({\mathbf{x}}^{\prime },t\right)+I_E(\mathbf{x},t)]$$

Full text

1 Introduction

Wilson and Cowan published the first of two classic papers in 1972 [1], wherein they obtained a set of coupled mean field equations providing a coarse-grained description of the dynamics of a network of excitatory and inhibitory neuron populations. The second paper in 1973 included spatial dependence [2]. Their work provided a demarcation in the practice of theoretical neuroscience. Prior to Wilson and Cowan, the nascent field of computational neuroscience was still grappling with how best to model neural dynamics. Thereafter, it began to coalesce onto a common set of neural network equations and studying their properties became the focus.

By the mid twentieth century, it became an established idea that at least some information processing in the brain is performed via global (population) activity (as opposed to the single neuron level) [3]. A number of studies looked into the dynamics of neuronal populations in order to find a theoretical framework for studying the collective behavior of neurons [4, 5, 6]. Given the complex nature of the brain and our lack of access to details of the microscopic processes, building a dynamical theory naturally called for the use of statistical approaches, including coarse-graining and mean field theory.

Course-graining considers the system at a lower resolution, thereby reducing the degrees of freedom of the system. However, some information will be lost. The strategy is to ensure that the lost information is not crucial to understanding the phenomenon of interest. The classic example is the dynamics of gas molecules in a closed room. If we are only interested in whether the room will be comfortable, then the coarse-grained measures of temperature, pressure, and relative densities of the constituent molecules (e.g. water, nitrogen, oxygen, carbon dioxide, etc.) is sufficient. The actual microscopic dynamics of the over $10^{23}$ molecules is irrelevant. However, in systems that are highly coupled, such as the brain, it is not clear what information is important. Mean field theory is a coarse-graining approach that captures the relevant mean dynamics while discounting (hopefully irrelevant) fluctuations around the mean. Mean field theory was originally developed for magnetism [7, 8], which involves the alignment of atomic spins under the competing effects of coupling between the spin magnetic moments of each atom and disorder quantified by the temperature. Reminiscent of neurons in the brain, the dynamics are highly complex since a perturbation of one spin may lead to unending reverberating loops where each spin influences other spins which in turn affect themselves. Mean field theory cuts through this complexity by imposing self-consistency between the role of a spin as the influencer of other spins and the one being influenced by other spins. More concretely, the mean field equations are derived by computing the mean state of one spin given fixed yet unknown values of all the coupled spins then equating the values of the other spins to that of the first spin. If the number of coupled spins is large enough and the spins are uncorrelated enough then the variance of the spin state will go to zero (by the central limit theorem) and the mean will be the only relevant quantity. The neglect of correlations and higher order statistics is the hallmark of mean field theory.

The original Wilson-Cowan equations were a coarse grained mean field system for a continuous activity variable representing the proportion of a local population of neurons that is firing or active at any moment of time. However, this interpretation is not rigidly adhered to and often the activity variable is deemed a physical quantity like the mean firing rate or an electrochemical potential of a single neuron or small group of neurons. The Wilson-Cowan and related equations became ubiquitous because they were able to model many elements of observed neural phenomena. Here we present a brief (and grossly incomplete) review of the history of some key ideas that led to their inception. We discuss when the assumptions of the equations are valid, when they break down, and how to go beyond them to take into account of fluctuations and correlations.

2 What are the Wilson-Cowan equations?

The use of neural network equations is pervasive. They have been referred to as Wilson-Cowan equations [1, 9, 10], Amari equations [11], Cohen-Grossberg equations [12], rate equations [13, 14, 15], graded-response equations [16, 17], continuous short term memory (STM) equations [18, 19, 20], and neural field equations [21, 22]. Here, we will adopt the term activity equations. Most of these equations come in one of three forms. The actual activity equations in Wilson and Cowan’s 1972 paper are

[TABLE]

where $a_{i}$ is neural activity of a local population indexed by $i$, $f_{i}$ is a nonlinear response or gain function, $w_{ij}$ is a coupling weight matrix, $I_{i}$ is an external input, $\tau$ is a decay time constant, and $r$ is a refractory time constant. Wilson and Cowan explicitly separated the effects of excitatory and an inhibitory neurons. Their 1973 paper introduced a continuum version:

[TABLE]

where $x$ is a spatial variable on a domain $\Omega$. The other two common forms are

[TABLE]

which Grossberg [20] calls the additive and shunting short term memory (STM) equations respectively. Grossberg and colleagues use (2) and (3) with excitatory and inhibitory effects separated. Amari and subsequent neural field papers use the continuum version of (2), often with excitation and inhibition lumped together so neurons no longer obey Dale’s law but simplifies the analysis.

The term Wilson-Cowan equations often refer to either (1) with $r=0$ or (2). For constant external input, these two equations are actually identical [20] as can be seen by setting $s_{i}=I_{i}+\sum_{j}w_{ij}a_{j}$, and taking the derivative:

[TABLE]

From this perspective, if $a$ is neural activity then $s$ is best viewed as a synaptic input or drive [23]. We will use the term Wilson-Cowan (WC) equations to exclusively refer to (1). For (2) we will use the term graded-response model. For a quasi-stationary state where $s_{i}\approx\sum_{j}w_{ij}f_{j}\left(s_{j}\right)+I_{i}$, $a_{i}$ can be interpreted as a neuronal mean firing rate

[TABLE]

and $f_{i}(s)$ is a single-neuron gain function (FI curve) specifying rate as a function of input. Equation (1) would then be considered a rate-based model equivalent to the graded-response model (2).

3 Before the Wilson-Cowan equations

There are three essential ingredients that activity equations include: 1) a continuous time dependent activity variable, 2) a linear weighted sum over these variables in the input, and 3) a nonlinear gain or activation function linking input to output. Underlying these ingredients are coarse-graining and mean field theory. Here we explore some of the history behind these ideas. A list of milestones is in the box.

It was well known throughout the last century that neurons fired action potentials if given sufficient stimulus. Thus it was reasonable to model them as binary state machines that activate when the inputs exceed a threshold. In 1943, McCulloch and Pitts [24] showed that a network of such binary threshold neurons with both excitatory and inhibitory connections is capable of producing any logical operation and thus performing universal computations. This launched the study of neural networks in computer science, artificial intelligence, and neuroscience.

Shortly thereafter a dissenting opinion emerged that sprouted a branch of research that would eventually lead to the Wilson-Cowan equations and dynamical systems modeling of the brain. Shimbel and Rapoport [25] argued that the McCulloch and Pitts “neural net” approach fell short of building plausible theories for neural systems. They argued that countless numbers of equivalent networks could produce the same desired output and that neural nets were not robust to the failure of a few neurons. Additionally, it was unlikely that genes predetermined the details of the microscopic structure of biological systems but rather imposed statistical traces on the macroscopic properties. They proposed that the goal should be to seek the statistical properties that govern the evolution and function of a biological system, rather than proposing specific network connectivity to perform a particular task.

They derived a dynamic equation governing the firing probability of neurons located in a given region of a brain network in terms of a recursive map. Like McCulloch and Pitts they assumed binary neurons and supposed that a particular synapse receives input from $n$ axon terminals, which they called “bulbs” partitioned into $p$ “bulb groups”. The probability that a neuron will fire depends on an integer number of bulbs exceeding a threshold $h$. Neuron firing dynamics are decomposed into a set of probabilities governing the connectivity and synaptic properties. The fraction of neurons at location $Q$ firing at time $t=1$ was presumed to be

[TABLE]

where $f_{0}(X)$ is the firing probability of neurons at $X$ on a domain $\Omega$ at time $t=0$, $O(Q,X)$ is the probability that a bulb group at $Q$ originates from $X$, $P(Q;n,p,h)$ is the probability that a neuron at $Q$ has a firing threshold $h$ and synapses of type $(n,p)$, and $P_{nph}[I_{0}(Q)]$ denotes the conditional probability that an active neuron at $Q$ is of $(n,p,h)$ type. A crucial idea for these equations is that the probability of firing at a given time only depends on the probability of firing in the previous time and not on any higher order statistic. This is an implicit mean field theory assumption. Although $P_{nph}[I_{0}(Q)]$ are derivable in principle, it would be a daunting task because the probability of activation depended on the microscopic details of synaptic configurations.

Beurle [4] introduced coarse-graining, which resolved the synaptic structure problem. He was inspired by experimental findings on the distribution of dendritic fibers as a function of distance from the cell body [26, 27]. Anatomical studies had started to focus on connectivity patterns of neuronal aggregates (as opposed to the connections of individual neurons), providing probability distributions instead of specific circuits [28]. Berule’s coarse-graining scheme represented neural activity as a continuous field in space and time and was the first neural field theory (see [22] for a review). Beurle assumed a network of purely excitatory cells with uniform thresholds and refractory period. A neuron becomes active if it receives input from a threshold number of active neurons. If $q$ is the threshold then the rate at which the cells become active at a time $t+\tau$ is

[TABLE]

where $F(x,t)$ is the proportion of cells active in any region per unit time, $k$ is a scale constant, $\tau$ is the cell time constant, $R$ is the proportion of cells which are sensitive to input, and

[TABLE]

is the proportion of cells which have an integrated excitation of $(q-1)$ and are ready to spike with the arrival of one further impulse. The mean number of activated neurons given some input $\bar{N}$ obeyed

[TABLE]

where $\alpha(t)$ is the temporal response function and $w(x)=be^{-|x|/a}$ is a spatially dependent connection function inspired by experimental studies [26, 27]. Again, this was a mean field theory formulation as activity depends only on past activity and not higher moments (cumulants) of the activity.

Beurle analyzed spontaneous random activity and waves. He was interested in waves of activity because of then-recent electroencephalograph studies [29, 30, 31, 32]. He showed that his model with minimal structure was capable of supporting plane, spherical and circular waves as well as vortices. In addition, he was driven by the idea that switching of neural waves may be related to shifting of attention (or other perceptions). However, his model also showed that the stationary fixed point is unstable; a slight deviation from the critical point either leads the activity to cease or saturate (due to refractoriness). Therefore, he erroneously concluded that there is no stable continuous random activity. Despite the achievements of the theory, the lack of inhibition as well as some technical flaws drove others to look for alternative formulations. Nevertheless, Beurle established most of the concepts that eventually led to the Wilson-Cowan equations.

The next milestone was Griffith [5, 6]. After establishing the role of inhibition in stabilizing brain activity [33], Griffith took a different approach and derived a reaction-diffusion equation for neural activity. He desired a classical field theory with a field representing the synaptic excitation or input to the cells and a source that quantifies the activity of the cells (e.g. spiking frequency). He built a phenomenological model by enforcing a set of constraints on the equations.

Griffith adopted a linear inhomogeneous differential equation for $\psi$ of the form:

[TABLE]

where $H$ is a linear operator and $\kappa$ is a constant. He defined $F$ as

[TABLE]

where $f$ is a nonlinear activation function and the temporal response function $\chi(\tau)$ obeys causality. Assuming that $\psi$ is differentiable and quasi-stationary ($\partial^{2}f/\partial\psi^{2}\ll 1$), Griffith obtained the following general form for $F$ by Taylor expanding to second order:

[TABLE]

He then assumed that the field operator $H$ obeyed translational and rotational invariance, which implied

[TABLE]

Substituting a field solution ansatz with constant velocity $v$,

[TABLE]

where $w$ is a connectivity function, into equation (14), he obtained

[TABLE]

for $w(p)=Ap^{-1}e^{-\frac{1}{2}\beta p}$, which Griffith argued was compatible with Sholl’s experimental findings [26]. Rescaling and assuming a shallow activation function where $\partial f/\partial\psi$ varies slowly, in the limit of $v\to\infty$ equation (16) can be reduced to

[TABLE]

where $\mu>0$ is a constant and $f$ is also redefined as a general nonlinear function. From (17) it is apparent that for spatial homogeneity one finds a mean field equation

[TABLE]

which is the graded-response model (2) for uniform connectivity, which as shown above is equivalent to the Wilson-Cowan equation.

Griffith showed that Beurle’s formalism is equivalent to his if the activation function $f(\psi)$ is adjusted appropriately. He originally assumed a linear response function, namely $f(\psi)=f(\psi_{c})+b(\psi-\psi_{c})$, but argued that a sigmoidal $f$ could have two stable critical points and an unstable fixed point in between that corresponded to Beurle’s unstable steady state. Griffith was the first to show the possibility of stable spontaneous activity depending on the gain function. However, he did not account for excitatory and inhibitory synaptic connections explicitly.

4 The Wilson and Cowan equations

By the mid 1960’s most of the concepts that would be incorporated into the Wilson-Cowan equation had been proposed although no one had put them together into a cohesive whole. Wilson and Cowan developed a coarse-grained description of neuronal activity where the distinction between excitatory and inhibitory cells was taken into account explicitly. They were motivated by physiological evidence from Hubel and Wiesel [34, 35], Mountcastle [36], Szentagothai and Lissak [37] and Colonnier [38], which suggested the existence of certain populations of neurons with similar responses to external stimuli.

In line with Sholl [3] and Beurle [4], Wilson and Cowan argued that a microscopic (single-neuron) description of neural activity is probably not well suited for understanding higher level functions that entail more complex processes such as sensory processing, memory and learning. Using dynamical systems analysis, they showed that their equations exhibit multi-stability and hysteresis, which could serve as a substrate for memory [39, 40], and limit cycles, where the frequency of the oscillation is a monotonic function of stimulus intensity. The combination of mathematical tractability and dynamical richness is the reason for the lasting legacy of their equations.

Wilson and Cowan derived their activity equations from first principles using a probabilistic framework on an aggregate of heterogeneous threshold neurons coupled with excitatory and inhibitory synapses with random connectivity by which spatial interactions could be neglected. They extended the model to the case of spatial inhomogeneity the following year [2].

They assumed that that neurons had a sensitive phase where input exceeding a threshold would cause them to fire and a refractory period in which they would not fire regardless of input. Defining $A_{i}$ as the proportion of cells of type $i\in\{E,I\}$ (excitatory or inhibitory) active at time $t$, the fraction that are refractory will be $\int_{t-r}^{t}A_{i}(t^{\prime})dt$, with refractory period $r$ and the fraction that are sensitive is then $1-\int_{t-r}^{t}A_{i}(t^{\prime})dt$.

The activation function $f(z)$ gives the expected proportion of cells that would respond to a given level of input $z$ if nonrefractory. If this and the fraction of sensitive neurons are independent, then the updated fraction of active cells would be

[TABLE]

where

[TABLE]

$Q_{i}$ is an external input, $\alpha(t-t^{\prime})$ is a response function governing the time evolution of the spike arrivals, and $W_{jk}$ is the average number of synapses from type $k$ to $j$. Wilson and Cowan noted that the input to a cell and the sensitive fraction could be correlated, which would violate mean field theory since both involve $A_{i}$. They argued that this correlation could be negligible for highly interconnected networks due to the presence of spatial and temporal fluctuations in the average excitation within the population, and also due to the variability of thresholds supported by experiments [41, 42, 43].

To estimate the shape of the activation function, they assumed that even if all the neurons receive the same input, any variability in other parameters such as the firing threshold or the number of afferent connections could lead to a probability distribution for the number of spiking neurons. Assuming that the variability mainly stems from the threshold, the expected proportion of neurons that receive at least a threshold excitation (per unit time) would be 111A similar argument holds for a uniform threshold and a distribution of neuronal afferent synapses, for which $f(x)=\int_{\theta/x(t)}^{\infty}c(w)\mathrm{d}w$ where $c(w)$ is the probability density for synaptic connections, and the lower bound on the integral is because all neurons with at least $\theta/x(t)$ connections would cross the threshold.

[TABLE]

where $p(\theta)$ is the threshold probability density. If $p(\theta)$ is unimodal then $f(x)$ would be a sigmoid function. Later, Feldman and Cowan [44] expanded on that interpretation showing that the activation function in the Wilson-Cowan equations can be regarded as the average of single-neuron activation functions.

Putting this all together gives

[TABLE]

After obtaining the general equations (21) and (22), Wilson and Cowan derived a set of differential equations that carried the biologically relevant aspects of their general equations. Defining time coarse-grained variables

[TABLE]

they argued that if $\alpha(t)\approx 1$ for $0\leq t\leq r$ and decays rapidly for $t>r$, then

[TABLE]

where $k$ and $r$ are constants. Inserting into (21) and (22), Taylor expanding, and rescaling gives

[TABLE]

After deriving the original mean field equations, Wilson and Cowan considered spatial inhomogeneity. Inspired by Beurle [4], they extended their original equations to:

[TABLE]

[TABLE]

where the activity variables are

[TABLE]

where $\alpha\left(t-t^{\prime}\right)=\alpha_{0}e^{-\left(t-t^{\prime}\right)/\tau}$, $W_{ij}\left(\mathbf{x}-\mathbf{x}^{\prime}\right)$ is the spatially dependent connection weight, $\rho_{i}(\mathbf{x})$ defines the density of neurons in a small volume around $\mathbf{x}$, and $n_{E}$ and $n_{I}$ represent the proportions of excitatory and inhibitory neurons activated per unit time. Although Wilson and Cowan developed the model to incorporate spatial dependence of connections, the formulation could go beyond spatial dependency as $\mathbf{x}$ could in general be any abstract quantity.

5 Graded Potential models

The artificial intelligence branch of neural networks following McCulloch and Pitts focused initially on the computational capabilities of the binary neuron, e.g. [45, 46, 47]. However, Hartline and Ratliff [48] considered the continuous firing rate of single cells in a neural network to successfully model the liminus retina, which could be directly measured. This was followed by Grossberg, whose initial goal was to understand how the behavior of an individual can adapt stably in real-time to complex and changing environmental conditions [20]. The core of his approach led to a class of continuous neural networks defined by the nonlinear coupling between activity and synaptic (adaptive) weights. He proved that the computational units of these networks are not individual activities or connections, but are the pattern of these variables (Grossberg,1968b, 1969a,1969b,1970a) [49, 50, 18, 18, 20].

Grossberg first considered model (2), which he termed the additive STM equation since the equation exhibits bistability as shown by Griffith. He then added dynamics for the synaptic weights to model long term memory. Inspired by the structure of the Hodgkin-Huxley model, Grossberg next derived an equation for neural networks that more closely modeled the shunting dynamics of individual neurons resulting in (3) [20]. The shunting STM model is approximated by the additive STM model when the activities $s_{i}$ in (3) are far from saturation. These networks are capable of performing a rich repertoire of behaviors including content addressable memory and oscillations [19]. Cohen and Grossberg (1983) proved that under general conditions including symmetric weights, general activity equations possess a Lyapunov function, indicating that all orbits will flow to local minima and can serve as a basis for models of associative memory.

In 1982, Hopfield using the analogy to the statistical mechanics of spin glasses, showed that a discrete time binary neuron model with symmetric connections has a non-increasing energy function (Lyapunov function) and thus stable attractors that can act as memory states [51]. The foundations of this work had been set in the early 1960’s, when several studies demonstrated that (artificial) adaptive networks could perform recognition tasks [52, 53]. These advances were followed by a number of studies on associative memory and pattern recognition, including those by Grossberg [49, 50, 18, 18] and Amari [54, 55, 11]. Little and Shaw [56, 57] pointed out the analogy between neural networks and spin glasses. Hopfield’s contributions [17] attracted a lot of interest from the physics community [58, 59], but it is rarely acknowledged that his discovery was a special case of the general Cohen-Grossberg theorem and that there was much work in this area that preceded his.

6 Spike response formalism

The Wilson-Cowan equations were highly impactful in modeling a number of neural phenomena such as pattern formation, waves, and slow oscillations. However, Feldman and Cowan showed that the Wilson-Cowan equations (1) are only valid for quasi-stationary activity [28]. Gerstner and van Hemmen [60] also showed that for stationary activity, any single-neuron model can be reduced to a network of graded-response neurons, but this is not true for coherent oscillations [60].

To address this deficiency, Gerstner sought a general formulation of global dynamics whereby one can systematically estimate the accuracy of rate-based models. To that end, he developed the spike-response formalism, which is a systematic derivation of effective population (macroscopic) dynamics, given that the single neuron (microscopic) dynamics are known [61]. His approach echoes the ideas underlying the Wilson-Cowan equations, with a focus on a realistic model of single-neuron dynamics.

In the spike-response formalism, the membrane potential of a neuron is modeled by the combination of the time dependent refractory response of the neuron to its own activity and the summed responses to the incoming spikes [62, 60]. As a result the synaptic potential of a single neuron can be described by the integro-differential equation

[TABLE]

where $h_{i}$ is the membrane potential, $\eta^{\mathrm{refr}}$ is the refractory function, $t_{i}^{f}$ denotes the spike times for neuron $i$, and $\kappa\left(s,s^{\prime}\right)$ denotes the response kernel in which $s$ is the time that has passed since the last postsynaptic spike represented by $S_{j}^{(F)}$ [62]. Different models of single-neuron dynamics can be reduced to a spike-response model with appropriate kernel functions [60, 63, 61].

In order to find the connection to rate-based models, Gerstner applied mean field theory. He considered a uniform population of neurons for which the activity $A_{p}(t)$ of a single pool $p$ is defined as the proportion of spiking neurons per unit time within the pool. Neurons in a pool are equivalent; connection weights and the response kernel $\kappa(s,s^{\prime})$ only depend the pool identity. However, the addition of noise, which turned out to be essential in his formulation, does cause variability among spike trains within a single pool. To formulate his pool dynamics, Gerstner assumed that the spike trains could be approximated as renewal processes if the synapses are weak $|W_{pq}|\ll 1$. This led to a closed set of mean field equations describing the synaptic input $h(p,t)$, pool activity $A_{p}(t)$ and the firing noise probability imposed by the dynamics. His formalism is suited to model neural network dynamics at arbitrarily short time scales.

In line with previous studies [44, 60], Gerstner showed that differential equation activity models, while excellent for modeling asynchronous firing, break down for fast transients and coherent oscillations. He derived a correction to the quasi-stationarity assumption, estimating the error of rate-based models. He showed that in order for the rate-based description to be valid one should have $\langle s\rangle\dot{A}(t)/A(t)\ll 1\,$ where $\langle s\rangle$ denotes the inter-spike interval (inverse of mean firing rate). In addition, he also re-derived the Wilson-Cowan equations by adjusting the appropriate kernels and provided a more rigorous derivation wherein the noise-induced neuron firing probability plays the same role as the sigmoidal activation function in the original Wilson-Cowan derivation.

7 Mean field theory for known microscopic dynamics

The Wilson-Cowan equations and variants have successfully modeled and predicted a variety of neural phenomena. However, the question remains as to how quantitatively accurate they are in modeling a neural system with known microscopic neuron and synaptic dynamics, the spike response formalism notwithstanding. Here, we explicitly derive a mean field theory for the neural activity of a deterministic network of coupled spiking neurons [64, 65].

Consider a network of $N$ conductance-based neurons:

[TABLE]

where $V_{i}$ is the membrane potential of neuron $i$, $m_{i}$ represents a single or set of auxiliary variables, $s_{i}$ represents the synaptic input or drive, $h$’s are continuous functions specifying the dynamics of the respective variabls, $\tau$’s are time constants, $w_{ij}$ is a matrix of synaptic weights between neurons, and $A$ is a function that activates whenever $V_{j}$ exceeds a threshold indicating an action potential. We do not explicitly distinguish between excitatory or inhibitory neurons but this is reflected in the parameter values, which can vary from neuron to neuron, and the synaptic weight matrix (e.g. obey Dale’s law).

If the individual neurons have limit cycle dynamics, as expected for spiking neurons, and the coupling between individual neurons is not excessively strong (i.e. a single spike does not strongly perturb the spiking dynamics, although many can and will), then the neuron dynamics can be reduced to a phase variable around or near the limit cycle [66]. The system takes the simpler form of

[TABLE]

where $\theta_{i}$ is the phase of neuron $i$, $F_{i}$ is the phase velocity, and $t_{s}^{j}$ are the spiking times of neuron $j$, which we set to $\theta_{j}=\pi$.

Our goal is to generate a coarse-grained mean field description of the network specified by (32) and (33). We quantify the neuron activity in terms of a density of the neuron phase:

[TABLE]

where $\delta(\cdot)$ is the Dirac delta function. This allows us to write

[TABLE]

which is also the firing rate of the neuron, i.e. the neuron velocity or flux at $\pi$. Inserting (35) into (33) gives

[TABLE]

We obtain the dynamics of $\eta_{i}$ by imposing local neuron conservation. Although $\eta_{i}$ is not differentiable, we can still formally write

[TABLE]

which is called the Klimontovich equation in the kinetic theory of plasmas [67, 68]. Equations (37) and (36) fully specify the spiking neural network dynamics but are no simpler to solve than the original equations. They need a regularization scheme to make them useful. One approach is to average over a selected population. This was the strategy employed in sections 3, 4, and 6 and in previous mean field treatments of networks [69, 70] where a population average is taken. In the limit of $N\rightarrow\infty$, the phases are assumed to be sufficiently asynchronous so that $\eta$ is differentiable and (37) becomes a well behaved partial differential equation. With the addition of Gaussian white noise, it takes the form of a Fokker-Planck equation.

An alternative is to average over an ensemble of networks, each prepared with a different initial condition drawn from a specified distribution [71, 72, 73, 64, 74, 75, 76]. Taking the ensemble average $\langle\cdot\rangle$ over (36) and (37) gives

[TABLE]

We see that the dynamics of the mean of $\eta_{i}$ and $s_{i}$ depend on the higher order moment $\langle F_{i}(\theta_{i},s_{i})\eta_{i}\rangle$. We can construct an equation for this by differentiating $F_{i}(\theta_{i},s_{i})\eta_{i}$, inserting (36) and (37) then taking the average again but this will include even higher order moments. Repeating will result in a moment hierarchy that is more complicated than the original equations.

However, if all higher order moments factorize into products of the means (i.e. cumulants are zero) then we obtain the mean field theory of the spiking network (32) and (33):

[TABLE]

which does not match any of the mean field activity equations in Sec 2. However, since $F_{j}(\pi,t)\langle\eta_{j}(\pi,t)\rangle$ is the ensemble mean firing rate of neuron $j$, then we see that (39) has the form of (2) with time dependent firing rate dynamics given by (38). If $\langle\eta_{i}\rangle$ were to go to steady state then (38) with $\partial_{t}\langle\eta_{i}(\theta,t)\rangle$ set to zero can be solved to yield

[TABLE]

where $C_{j}$ is determined by the normalization condition $\int_{\theta}\eta_{j}(\theta)d\theta=1$. We can then obtain (2) with gain function

[TABLE]

However, (38) is dissipation free and thus will not relax to steady state from all initial conditions without a noise source or finite size effects [72, 73, 64].

The formalism can be applied to any type of neuron and synaptic dynamics. The neuron model (32) and (36) was chosen for simplicity and also so that the mean field equations would be similar to those of Sec 2. The inclusion of higher order time derivatives or nonlinear terms could result in very different mean field equations. A similar argument to arrive at mean field theory could be applied with the population average, although this will be further complicated by the neuron dependent weight $w_{ij}$, which we will address below. No matter what average is used, we still do not know if or when mean field theory is valid. To answer this question we need to compute the corrections to mean field theory and see when they are small.

8 Beyond mean field theory

Consider again the Klimontovich formulation of the spiking network equations (37) and (36)

[TABLE]

where we have expressed the initial conditions as forcing terms and $\beta=1/\tau_{s}$.

We quantify the ensemble average by defining an empirical probability density functional for $\eta_{i}(\theta,t)$ and $s_{i}(t)$ as functions constrained to (40) and (41), to wit

[TABLE]

${\cal P}[\eta,s|\eta^{0},s^{0}]$ is a Dirac delta functional in function space. This formal expression can be rendered useful with the functional Fourier transform expression of the Dirac delta functional [77]:

[TABLE]

where $\tilde{\eta}_{i}(\theta,t)$ and $\tilde{s}_{i}(t)$ are transform variables, called response fields, and the functional integral is over all possible paths or histories of these functions. Applying the functional Fourier transform gives the functional or path integral expression

[TABLE]

Assuming fixed $s_{i}^{0}$ and a distribution for initial density ${\cal P}[\eta^{0}]$, we can integrate or marginalize over the initial condition densities to obtain

[TABLE]

If we set $\eta_{i}^{0}(\theta)=\delta(\theta-\theta_{i}(t=0))$, the distribution over initial densities is given by the distribution over the initial phase $\rho_{i}^{0}(\theta)$. Thus we can write $\int{\cal D}\eta^{0}P[\eta^{0}]=\int\prod_{i}d\theta\rho_{i}^{0}(\theta)$. The initial condition contribution is given by the integral

[TABLE]

Hence, the system given by (40) and (41) can be mapped to the distribution functional

[TABLE]

with action $S=S_{\eta}+S_{s}$ given by

[TABLE]

The exponential in the initial data contribution to the action (which is a generating function for a Poisson distribution) can be simplified via the Doi-Peliti-Janssen transformation [71, 72, 73, 64, 74, 75, 76]:

[TABLE]

resulting in

[TABLE]

where we have not included noncontributing terms that arise after integration by parts.

The probability density functional (44) with the action (46) and (47) describes the ensemble statistics for the entire spiking and synaptic drive history of the network. By marginalizing over initial conditions, the functional density is no longer sharply peaked in function space. The expression cannot be evaluated explicitly but can be computed perturbatively by employing Laplace’s method and expanding around the extremum of the integrand. This extremum is given by the point in function space with minimal variation, $\delta S=0$ (i.e. the least action principle). The resulting equations of motion, if they exist, are the mean field equations.

The accuracy of mean field theory can be assessed by computing the variance around the mean field solution using perturbation theory. The inverse network size $1/N$ is a candidate expansion parameter. The action $S$ is comprised of the transformed phase density $\psi_{i}$, which is not differentiable, but the sum over neurons provides an opportunity to regularize it. Consider first a homogeneous network, $F_{i}=F$, and uniform global coupling, $w_{ij}=w$. Then (40) and (41) become

[TABLE]

The dynamics of all neurons are identical. We can drop the neuron index and rewrite the action as

[TABLE]

Applying the least action principle, which amounts to setting the functional derivatives to zero [74, 65] returns the mean field theory equations

[TABLE]

where

[TABLE]

The factor of $N$ in the action indicates that deviations away from mean field will be suppressed as $N\rightarrow\infty$. A perturbative expansion in $1/N$ can be constructed using Laplace’s method to compute all higher order cumulants, which go to zero as as $N\rightarrow\infty$ [64, 65].

The same argument can be applied to a heterogeneous network with spatially dependent coupling. If the coupling function and neuron specific parameters are continuous functions of the neuron index, then mean field theory applies for the density and synaptic drive averaged over a local population of neurons at each location in the network in the limit of $N\rightarrow\infty$ [65].

9 Final words

We have provided a brief and biased view of the ideas behind the Wilson-Cowan equations and their descendants. We show that one can arrive at a similar set of activity equations with first order time derivatives, nonlinear activation (gain) function, and a linear sum over inputs, from various starting points. Despite the fact that the equations are not generally quantitatively accurate, they seem to capture many facets of neural activity and function. They are a set of effective equations that have stood the test of time. Mean field theory is implicitly assumed in all of these equations. Correlations are deliberately ignored to make the equations tractable. However, they are omnipresent in the brain and there is much ongoing effort to account for them. Nonetheless, we predict that forty seven years from now, the Wilson-Cowan equations will still have great utility.

10 Acknowledgments

This work was supported by the Intramural Research Program of the NIH, NIDDK.