## Abstract

Population rate or activity equations are the foundation of a common approach to modeling for neural networks. These equations provide mean field dynamics for the firing rate or activity of neurons within a network given some connectivity. The shortcoming of these equations is that they take into account only the average firing rate, while leaving out higher-order statistics like correlations between firing. A stochastic theory of neural networks that includes statistics at all orders was recently formulated. We describe how this theory yields a systematic extension to population rate equations by introducing equations for correlations and appropriate coupling terms. Each level of the approximation yields closed equations; they depend only on the mean and specific correlations of interest, without an ad hoc criterion for doing so. We show in an example of an all-to-all connected network how our system of generalized activity equations captures phenomena missed by the mean field rate equations alone.

## 1.  Introduction

Modeling the brain is confounded by the fact that there are a very large number of neurons, and the neurons are heterogeneous and individually complex. Given current analytical and computational capabilities, we can either study neuronal dynamics in some biophysical detail for a small or medium set of neurons or consider a large population of abstract simplified neural units. We then can extrapolate only to the desired regime of large numbers of biophysical neurons. In particular, there is a dichotomy between network models that incorporate Hodgkin-Huxley or integrate-and-fire spiking dynamics and models that include only the rate or activity of neural units. While the rate description has yielded valuable insights into many neural phenomena, it cannot describe physiological phenomena thought to be important for neural processing such as synchronization, spike-time-dependent plasticity, or any correlated activity at the spike level. Likewise, it is difficult to analyze or simulate a large network of spiking neurons. Our goal is to derive an intermediate description of neural activity that is complex enough to account for spike-level correlations yet simple enough to be amenable to analysis and numerical computation for large networks.

Rate or activity equations have been a standard tool of computational and theoretical neuroscience, early important examples being the work of Wilson and Cowan, Cohen and Grossberg, Amari and Hopfield (Wilson & Cowan, 1972, 1973; Amari, 1975, 1977; Hopfield, 1984; Cohen & Grossberg, 1983). Models of this type have been used to investigate pattern formation, visual hallucinations, content-addressable memory, and many other topics (Ermentrout & Cowan, 1979; Hopfield, 1984; Ermentrout, 1998; Bressloff, Cowan, Golubitsky, Thomas, & Wiener, 2001; Coombes, 2005). Naturally these equations are so-called because they describe the evolution of a neural activity variable often ascribed to the firing rate or synaptic drive of a population of interacting neurons (Ermentrout, 1998; Gerstner, 2000). These equations are considered to represent the neural dynamics averaged over time or population of a more complicated underlying process. In general, these activity equations make an implicit assumption that correlated firing is unimportant. They are a mean field theory that captures the dynamics of the mean firing rate or activity that is independent of the influence of correlations, which in some cases may alter the dynamics considerably. As an example, the effects of synchrony, which have been proposed to be important for neural processing (Gray & Singer, 1989; Beshel, Kopell, & Kay, 2007), are not included. Here, we give a systematic prescription to extend rate models to account for these effects.

An analogy for our problem and approach can be made to the field of equilibrium statistical mechanics. The statistics of such systems (e.g., the Ising model) in thermal equilibrium are described by a partition function, an integration over all configurations available to the system. For the Ising model, this refers to all possible configurations of the individual spins. The partition function is akin to the generating function for a statistical distribution from which the moments or cumulants can be obtained. For the Ising model, the first moment corresponds to the mean magnetization, and the second moment describes the mean correlation between the spins. The linear response of the system is the magnetic susceptibility, which describes the reaction of the system to an external input. In general, the partition function cannot be summed or integrated explicitly. However, these moments can be obtained perturbatively by using the method of steepest descents to approximate the partition function. This then yields a systematic expansion, and the lowest order is called mean field theory, since all higher cumulants are zero. By computing the expansion to higher order, the effects of correlations and fluctuations can be included.

This procedure requires full knowledge of the underlying microscopic theory that is to be averaged over. In neuroscience, the underlying model is not completely known; it would require full knowledge of the different types of neurons, their membrane and synaptic kinetics, and their synaptic connectivity. However, given a particular mean field theory, one can ask about the minimal constraints this theory places on the microscopic theory and its asymptotic expansion. Thus, although the full microscopic theory cannot be reconstructed, by constraining the expansion, the mean field theory can dictate the minimal structure of any extension of a set of rate equations. In this letter, we consider a well-known neural rate equation and deduce the minimal structure we expect for a consistent extension that includes correlations.

Buice and Cowan (2007) previously adapted a path integral formalism used in nonequilibrium statistical mechanics (Doi, 1976a, 1976b; Peliti, 1985) to analyze the dynamics of a Markov model for neural firing. They derived a generating functional (expressed as an infinite-dimensional path integral), which is specified by an “action” for the complete dynamical distribution of the model, and showed that the mean field theory for that system corresponded to a Wilson-Cowan type of rate equation. They then analyzed the scaling properties for the correlations near criticality. They showed how mean field theory could be corrected by using steepest descents to generate a systematic expansion that describes the effect of correlations. Here, we show that a moment hierarchy can be constructed by taking explicit averages of the Markov model. Each equation in the moment hierarchy is coupled to higher moments in the hierarchy. The hierarchy can be made useful as a calculational tool for the statistics of the dynamics if it can be truncated. We show that the moment hierarchy and the generating functional are equivalent and that the equations of the hierarchy are the equations of motion of the action in the generating functional. The truncation condition for the perturbation series of the path integral is also a truncation condition for the hierarchy. This provides for both a compact description of network statistics and a natural truncation or closure condition for a moment hierarchy. We can also show using the path integral formalism that the Markov model is a natural minimal extension of the Wilson-Cowan rate equation.

Approaches to neural network modeling using statistical mechanics are not new (Cowan, 1991; Hopfield, 1982, 1984; Ohira & Cowan, 1993; Peretto, 1984; Amit, Gutfreund, & Sompolinsky, 1985). Those works were largely concerned with models adhering to detailed balance, whereas we make the explicit assumption that neural dynamics admits an absorbing state that violates detailed balance. In the absence of internal activity and external stimulation, there will be no activity in the network. Other studies using a stochastic description of neural dynamics have considered the neurons in a background of Poisson activity with disorder in the connectivity (Amit & Brunel, 1997a, 1997b), or considered neural activity as a renewal process (Gerstner, 1995, 2000; Gerstner & Kistler, 2002). Van Vreeswijk and Sompolinsky (1996, 1998) demonstrated that disorder in network activity can arise purely as a result of disorder in the connectivity, without stochastic input. Kinetic theory and density approaches are investigated in Nykamp and Tranchina (2000), Cai, Tao, Shelley, and McLaughlin (2004), Ly and Tranchina (2007), and mean field density approaches to the asynchronous state appear in Abbott and Van Vreeswijk (1993) and Treves (1993). Golomb and Hansel (2000) study synchrony in sparse networks via a reduction to a phase model. Fokker-Planck approaches for networks appear in Fusi and Mattia (1999), Brunel and Hakim (1999), Brunel (2000), and Brunel and Hansel (2006). Responses of single neurons driven by noise appear in Plesser and Gerstner (2000), Salinas and Sejnowski (2002), Fourcaud and Brunel (2002), and Soula, Beslon, and Mazet (2006). Approaches to correlated neural activity including finite size effects appear in Ginzburg and Sompolinsky (1994), Mattia and Del Giudice (2002), Soula and Chow (2007), and El Boustani and Destexhe (2009). El Boustani and Destexhe (2009) develop a moment hierarchy for a Markov model of asynchronous irregular states of neural networks that is truncated through a combination of finite size and a scaling condition. Our work extends the results of Ginzburg and Sompolinsky (1994) by providing the systematic higher-order expansion without explicitly requiring the consideration of the rest of the hierarchy. We also provide conditions for the truncation of the expansion and consider the network response to correlated input. Our expansion is not a finite size expansion, although it can reduce to a finite size expansion under certain conditions (such as normalized all-to-all connectivity in the network).

In section 2, we revisit the original Wilson-Cowan framework and propose a Markov model that has the minimal stochastic dynamics to produce the Wilson-Cowan equations. This will be more rigorously justified in section 4. Section 3 presents the derivation of a moment hierarchy for this Markov model. After truncating, we provide a posteriori justification for the truncation. It will be seen in section 4 that the validity of this truncation was in fact natural and did not require ad hoc assumptions. The truncation conditions turn out to be related to the proximity to a bifurcation point as well as the extent of connectivity in the network. We also make more precise the sense in which our Markov model is “minimal” by introducing the path integral formulation. The field theory formalism that appears in this letter arose in the context of reaction-diffusion problems. (See Janssen & Tauber (2005) and Tauber, Howard, & Vollmayr-Lee (2005) for reviews of this formalism applied to reaction-diffusion and percolation processes.) We demonstrate a simple example all-to-all system in section 5 and show some simulation results.

## 2.  Rate Equations Reconsidered

We consider a population rate equation of the form
2.1
where a(x,t) is a measure of local neural activity at location , α is a decay constant (often equated with the membrane time constant or a synaptic time constant), f(s) is the firing rate or gain function describing how input affects the activity, I(x,t) is a time-dependent external input to location x, and w(x,y) is a weight function describing how a neuron at location y affects neurons at location x. Equation 2.1 is the standard form of rate equation seen in the Wilson-Cowan equations (Wilson & Cowan, 1972, 1973). While we use this form in our letter, our results can be adapted to any other type of rate or activity equation. The exact nature of the activity a(x,t) is open to interpretation (Gerstner, 1995; Ermentrout, 1998). It could be envisioned as the time average, population average, or ensemble average of neural firing or synaptic activity. In any case, the picture is that of some underlying process whose degrees of freedom have been marginalized to generate an effective theory with simpler variables.

We imagine that the typical rate equation is produced by some marginalization process over both disorder and extra degrees of freedom. Hence, it may be possible to derive a generating function for the statistics of the marginalized process. The lowest order in the steepest-descent expansion of the generating function describes mean field theory, which gives the rate equation. Since the operation of marginalization is dissipative, we cannot recreate the underlying microscopic process exactly with only the rate equation alone. However, the mean field theory places constraints on the structure of the dynamics, enabling us to investigate the structure of higher-order statistics implied by the structure of mean field theory. In the original derivation by Wilson and Cowan (1972, 1973), the activity variable was presumed to describe the fraction of neurons firing per unit time within some region of the brain. Two main features of this interpretation bear emphasizing. First, the rate equations were originally understood to be equations providing the dynamics of the probability that a neuron at x will fire at time t. There is therefore an implied underlying probabilistic model. Second, the probability a(x,t) applies to all neurons within some region of the brain, not just a single neuron. Thus, a spatial averaging component is implicit in the equations. The original Wilson-Cowan rate equations thus described the dynamics of the probability for a neural aggregate in the brain. Another feature implicit in the Wilson-Cowan equations is that these probabilities are independent for each neuron. This implies that the Wilson-Cowan picture is one in which neurons fire with Poisson statistics with firing rate determined by a(x,t), a picture supported by neural recordings (Softky & Koch, 1993).

Given this perspective, one might consider what processes may underlie rate equations. One route is to treat the fundamental, small-scale dynamics as a probabilistic process, for example, a Markov process. In this case, the basic description for neural activity will be provided by a master equation governing the evolution of probabilities for different neural configurations. This route obscures the source of uncertainty in neural activity in favor of directly modeling the probabilistic activity. This tactic has been used to model the so-called asynchronous irregular states seen in some neural models (Van Vreeswijk & Sompolinsky, 1996; El Boustani & Destexhe, 2009). Another route would be to employ the strategy of kinetic theory (Nicholson, 1992; Ichimaru, 1973) and define a continuity (i.e., Klimontovich) equation for the probability density of a network of deterministic neurons. (For an example of this approach applied to coupled oscillators, see Hildebrand, Buice, & Chow, 2006; Buice & Chow, 2007.) In that work, the probabilistic aspects of the model arise from the distribution of driving frequencies and initial conditions. Ultimately the difference in the two approaches is the origin of stochasticity—whether it is implicit in the dynamics of the neurons or an emergent property of the interaction of deterministic neurons (e.g., chaos). In either case, the final product is an effective stochastic dynamical system. In this letter, we will follow the approach of assuming an underlying probabilistic model given by a master equation, so that any emergent chaos has already been absorbed into the dynamics. We then seek a minimal stochastic model that will produce the Wilson-Cowan rate equation at the mean field level. We can then formulate equations governing the fluctuations of this model. In this section, we motivate such a minimal model qualitatively, leaving a more rigorous approach for section 4.

Our primary interest is in tracking the statistics of active neurons. A simple master equation whose mean field statistics for neural activity is represented by the Wilson-Cowan rate equations is given by
2.2
where is the probability of the network having the configuration described by at time t, and ni is the number of active neurons at location i. Neurons relax back to the inactive or quiescent state with rate α, which appears as a decaying transition in the master equation. Configurations and denote the configuration , where the ith component is ni ± 1, respectively. The rate at which a neuron at location i becomes active is given by the firing rate or gain function , an implicit function of the weight function wij and external inputs Ii. One of the crucial elements of the ensuing calculation is making the connection between the gain function , which appears in equation 2.2, and f(Ii(t) + ∑jwijnj), which appears in equation 2.1.
In general, we cannot solve for in equation 2.2 explicitly. One strategy is to derive an expansion of in terms of its moments 〈ni(t)nj(t′)nk(t′) ⋅ ⋅ ⋅ 〉 where the expectation value is over all statistical realizations of the Markov process. The first moment,
2.3

is a measure of the mean activity in the network. We obtain an equation for ai(t) by multiplying equation 2.2 by ni and taking the sum over all configurations . However, this equation is not closed (i.e., it depends on the second and possibly higher moments). An equation for the second moment can be similarly constructed by multiplying equation 2.2 by ninj and summing over all configurations. The resulting equation will depend on the third and higher moments. Continuing this process will result in a moment hierarchy with as many equations as there are locations, which could be infinite. In general, no finite subset of this hierarchy is closed. This means that if we wish to have a closed set of equations, we need to make some approximation that allows us to truncate the hierarchy.

The simplest way to close or truncate the moment hierarchy is to assume that all the higher-order moments factorize into products of ai(t). This is the naive mean field assumption where all cumulants are zero. For example, the second cumulant (i.e., variance) 〈ni(t)nj(t′)〉 − ai(t)aj(t′) is set to zero. For our master equation, 2.2, this assumption yields (see the computation in the next section)
2.4
which is similar in form to the rate equation 2.1 for a discrete domain and with the firing rate function given by . However, the fact that statistics in the brain are observed to be near Poisson is devastating for this naive mean field assumption because every cumulant is comparable to the mean, implying that any such truncation of the resulting hierarchy is not justifiable.
Here, we describe an alternative means of truncating the hierarchy consistent with near-Poisson firing statistics. In this case, we observe that the equation for ai(t) can be written as
2.5
where g is some functional dependent on the higher moments in the hierarchy. In order to choose a reasonable approximation and get a finite system of closed equations, we must identify a finite set of higher moments in terms of which we can express the remaining moments. We are guided by the indication that neuron firing statistics are near Poisson. Not coincidentally, the solution to the master equation in the case where the gain function Fi is constant or linear is exactly Poisson with mean rate (i.e., stochastic intensity) determined by the Wilson-Cowan rate equation. In order to truncate the hierarchy, we perform a change of variables to measure the deviations of each cumulant from the value under Poisson statistics. This new hierarchy is truncatable, as will be demonstrated a posteriori. From the perspective of solving the master equation, this new hierarchy is the natural one because the underlying statistical model is a point process.

The moment hierarchy approach does not make any approximation. It is a change of variables from the distribution to moments of that distribution. The approximation arises when we truncate this hierarchy in order to render the equations tractable. The simplest truncation is mean field theory. The first-order corrections to mean field theory are given by truncating at the next order. Truncation of the moment hierarchy requires some justification. We will demonstrate below that this justification in the neural case may be provided by the large spatial extent of neural connectivity and the distance of the system from a bifurcation.

## 3.  Truncation of the Moment Hierarchy

We will derive a moment hierarchy from our master equation, 2.2, and then show how it can be truncated. To get an equation for the first moment ai(t), we multiply the master equation 2.2 by ni and sum over all configurations :
3.1
The first two terms on the right-hand side simplify to
3.2
The first equality results from reindexing the summation over nk from (0, ∞) to (1, ∞). We leave the summation indicated as over all configurations because the factor of nk prevents the 0 term from contributing. We have also separated out the terms where i = k; the only term that survives is one of these. Note that we have made no approximations thus far. The terms involving the function take the form
3.3
Unlike the first term, we cannot directly represent this term as a function only of ai(t) due to its nonlinear nature. The equation for the mean is therefore
3.4
where will include higher-order moments.
We continue by constructing an equation for the second moment,
3.5
and third moment,
3.6
to obtain
3.7
3.8
Since we expect solutions to be near Poisson, we transform the hierarchy to describe the departure of moments from Poisson statistics. For a Poisson distribution, the cumulants are all equal to the mean of the distribution. Hence, we introduce what are called normal-ordered cumulants, which measure the “deviations” of the cumulants from Poisson values. The first normal-ordered cumulant is the same as the first moment ai(t). The next two are given by
3.9
and
3.10
The normal-ordered cumulants can be computed using a recursive algorithm. The algorithm involves replacing all moments Nijk⋅⋅⋅ with “normal-order-corrected” moments recursively by subtracting terms with coincident or “contracted” indices, which reduce the order of the moment. For example, the ordinary second cumulant is simply Nijaiaj. To compute the normal-ordered version, we replace the moments appearing in the expression with the normal-ordered-corrected forms, that is, set NijNijaiδij. The term subtracted results from the contraction of the i, j indices: NijNiδij = aiδij. The third cumulant is more complicated but follows the same strategy. The important thing to note here is that the higher moments must be independently corrected for each group of contracted indices. The ordinary third cumulant is given by
3.11
To obtain the normal-ordered cumulant, we first make the replacement NijkNijkNijδjkNikδijNjkδikaiδijδjk. We then must correct for all appearances of the second moment Nlm, resulting in equation 3.10. This algorithm systematically removes the underlying Poisson contributions (at all tensor ranks) and leaves us with the normal-ordered cumulants. For a Poisson distribution, all Cijk⋅⋅⋅ = 0 except the first, ai. An alternative rationale for normal ordering will appear in section 4.
Transforming the first three equations of the hierarchy, 3.4, 3.7, and 3.8 yields
3.12
3.13
3.14
where we must reexpress the expectation values involving Fi in terms of the normal-ordered cumulants.
Since is defined in terms of the vector of active neuron numbers, , its expectation value will be naturally expressed in terms of the moments, Nijk⋅⋅⋅, as given by the Taylor expansion of ,
3.15
We have implicitly defined the notation that Fjk⋅⋅⋅i is the derivative of Fi with respect to nj, nk, …. Note that this expansion applies only to the expectation value of . We need to reexpress this series as an expansion in terms of the normal-ordered cumulants. This transformation of variables will rearrange the terms and result in a new series with new coefficients that sums to the same result as the original series. For example, every term proportional to Nijk⋅⋅⋅ will contribute a term proportional to aiajak ⋅ ⋅⋅ to the normal-ordered cumulant expansion. This means that we can write the expansion in the form
3.16
as we would expect. However, there are also contributions from the normal ordering corrections. The simplest are those that arise due to every index being coincident or contracted at each order. This produces a term linear in ai. These corrections form the series
3.17
where by Fjm we mean the mth derivative of with respect to nj. Were a polynomial, this procedure would truncate at the highest order of the function. However, for an arbitrary general function, these corrections quickly become unwieldy as one proceeds through the orders of the expansion to include corrections to the terms that go as aiaj, aiajak, ….
Our perspective has been to interpret the Wilson-Cowan equation as describing the mean field of some Poisson process with activating and decaying transitions. Hence, the Wilson-Cowan equation should be the mean field solution of the normal-ordered cumulant hierarchy, equation 3.12. However, the gain function in the Markov equation is not the same as the gain function in the Wilson-Cowan equation. The Wilson-Cowan gain function is the normal-ordered version of the Markov gain function. Thus, in order for the Wilson-Cowan gain function to have the form of f(si), where si = ∑iwijaj + Ii(t), we assume that this resummation produces
3.18
where the higher-order terms (h.o.t.) are dependent on the higher normal-ordered cumulants according to the Taylor series expansion of f(s). It will be seen in section 4 that there always exists a master equation gain function F such that f is expressible in this form and the resummation works to produce the same f(s) (and derivatives thereof) at every order in this expansion.
We now return to the series 3.17 to consider terms with precisely one factor of Cij(t). At each order m in the series (i.e., the order that contains the mth moment), there are m(m − 1)/2 terms that have one factor of Cij(t). These terms sum to give
3.19
which can be rewritten as
3.20
Other terms are at least second order in Cij, higher normal-ordered cumulants, or corrections from normal ordering. Our first generalized activity equation, excluding terms dependent on third and higher-normal ordered cumulants, is therefore
3.21
We can take this same approach in order to compute an equation for Cij(t) and obtain
3.22
Equations 3.21 and 3.22 constitute a closed set of equations for the mean and variance of a Wilson-Cowan network of neurons. ai(t) represents the Poisson rate, and Cij(t) measures the deviation of the variance from Poisson statistics. These equations are the minimal consistent extension of the Wilson-Cowan rate equations to include the effects of higher-order statistics. Higher-order corrections can be incorporated by adding terms involving higher-order cumulants into equations 3.21 and 3.22 and including equations for these higher-order cumulants. We also note that the gain function need not be analytic everywhere for this expansion to work. It can contain a countable number of noncontinuous or nondifferentiable points. The equation would be corrected with the inclusion of impulse function terms at these singular points.

An immediate noteworthy consequence is that Cij(t) will have substantial input only when the activity is such that f′(s) is large. As an example, suppose f(s) is a simple sigmoid function. In this case, f′(s) is peaked at threshold (where we define threshold to be the central point of half maximum) and zero far away from threshold. Reasonably we have the result that correlated activity will increase only when the input to a neuron is near threshold. If the slope of the sigmoid is such that f(s) is a step function or near a step function, then Cij(t) will receive input only when the activity is precisely near threshold. Also notice that the strength of the input to Cij(t) is proportional to the weight wij between the neurons in question, as well as the mean activity. An initial check on the equations is that Cij decouples from ai in the case where f(s) is linear or constant.

To consider the dynamics of large-scale neural activity, we can take the continuum limit of these equations to get equations for the mean activity a(x,t) and correlation C(x, y, t):
3.23
and
3.24
These are the generalized activity equations. Had we wished to include even higher moments, we could have continued through the hierarchy. For simplicity of illustration, we truncate the hierarchy at this level.

### 3.1.  Criticality and Truncation of the Hierarchy.

Although we have derived equations for the mean activity and equal-time correlation, there are some outstanding issues. The primary concern that must be addressed is that we require some justification for the truncation of the hierarchy at the level of the two-point correlation function Cij(t) instead of allowing higher moments to interact with the mean activity.

Consider the mean field equation without the correction due to correlated activity 2.4. Define a0i to be some steady-state solution to this equation and linearize equation 2.4 around this solution. The perturbations δai(t), from this steady-state solution obey the equation
3.25
We rewrite this equation as
3.26
where the matrix Γij is defined by
3.27
If all of the eigenvalues of Γij are positive, then the solution a0i is stable. Likewise, negative eigenvalues indicate instability. Criticality is the condition of marginal stability, in which one or more of the eigenvalues are 0.
Returning to the equation for Cij(t), we see
3.28
We assume that the mean field solution a0i is stable. In addition, we assume, per the truncation hypothesis, that the steady-state value of Cij(t) does not appreciably alter either ai or, therefore, the matrix Γij. In this case, Γij has all positive eigenvalues and is diagonalizable. Define
3.29
to be the diagonalization of Γij. We also define the shorthand
3.30
for the driving terms in equation 3.22. The steady-state solution is given by
3.31
Notice that each term contributing to the magnitude of C0ij is attenuated by a sum of eigenvalues. The magnitude of the eigenvalues λi determines the distance of the system from a bifurcation or criticality; it is the distance from the onset of an instability. Thus, the further the system is from criticality, the more attenuated the fluctuations and the more justified we are in truncating the hierarchy. Conversely, the closer the system is to criticality, the more the approximation breaks down. At criticality, this solution, 3.31, becomes singular. This is an indication that criticality is a fluctuation-dominated, as opposed to mean field–dominated, regime. A similar argument will extend to any equation in the hierarchy, and we are left with an intuitively satisfying result: stability smoothes out fluctuations. This argument is what allows us to disregard the effects of still higher moments on the mean activity and truncate by considering only the two-point correlation function's effect on the mean.
It is also worth noting that the eigenvalues relevant for the dynamics of the two-point correlation Cij are the sums of the eigenvalues of the mean field equation, λi + λj. In the case that a0i is stable, not only will Cij be stable, but it will relax to equilibrium faster in general than the mean field solution. In kinetic theory, this is akin to the Bogoliubov approximation, in which the collision term is computed by solving for the steady state of the two-point correlation on the assumption that the correlation function reaches steady state on a timescale shorter than the mean field. In our case, this approximation leads to
3.32
Note immediately that the input from correlated activity decouples in the case that the firing rate function is in a linear or constant region, since the coupling is proportional to the second derivative of f(s).
The extent of neural interconnections also has an effect on the size of correlations and the ability to truncate the hierarchy. Consider that neuron i has Ni presynaptic neurons (i.e., neurons for which wij ≠ 0). Further, let the average connectivity weight over all inputs be w0. In this case, we can approximate wijw0/Ni. The steady-state value of Cij(t), which is determined by the driving term Aij(t), is seen to be inversely proportional to the number of presynaptic neurons due to the linear dependence on wij of Aij. In the most extreme case, the number of presynaptic neurons is the entire network, so Ni = N, and the correlation function Cij(t) becomes simply a finite size effect, going as 1/N. Smaller system sizes in general will have larger correlations, which is intuitively sound. More generally, as long as we can bound the total input to any given neuron, we can define Nm = minNi and scale all weights so that they can be written as
3.33
where wM is the maximum total input to any given neuron. This allows us to treat Nm as an effective system size. Larger Nm reduces the effects of fluctuations at a given distance to a bifurcation.

We have two competing effects. On one hand, we have the system size governing the magnitude of correlations. On the other, the distance of the system to a bifurcation likewise affects the size of fluctuations. The relative trade-off of the two, from the definition of C0ij, is given by the product of the smallest eigenvalue of Γij and the number of presynaptic neurons Ni. We will demonstrate this relationship more precisely when considering the all-to-all network in section 5.

Finally, it is worth pointing out the effect of input upon the hierarchy. If the input is another Poisson process, then the only equation affected is the one for ai(t). The higher equations in the hierarchy are affected by this input only through its effect on the mean activity ai(t) and the firing rate function f(s). In general, this suggests that large external inputs will actually reduce fluctuations, depending on the form of f(s), in the sense of driving the system toward Poisson-like behavior, a reasonable result. In particular, if f(s) is a saturating function, then the correlations will decouple from the equation for ai(t) and the source terms for higher correlations will be driven to zero, leaving the system described completely by the rate equations. The analogous situation for a ferromagnet is driving the system with a large external magnetic field.

## 4.  Path Integral Solution of the Master Equation

We have thus far demonstrated how a minimal Markov model consistent with the Wilson-Cowan rate equation can be used to derive generalized equations in a hierarchy of moments. Although we truncated this hierarchy at second order, one can in principle truncate at any desired cumulant, although the calculations become successively more cumbersome. Here we show that the moment hierarchy is equivalent to a path integral or field-theoretic approach, which systematizes the perturbation theory for the statistics of the network by providing rules for the construction and evaluation of the cumulants. Another major benefit is that it provides a systematic means for obtaining moment truncations or closures. The path integral representation of the master equation, 2.2, was derived by Buice and Cowan (2007) by modifying a method originally developed for reaction diffusion systems (Doi, 1976a, 1976b; Peliti, 1985). We quickly review the representation and then detail how the generalized equations can be derived from this representation.

The moment-generating function for the probability density is given by
4.1
where the sum is over all configurations of . Moments of P are obtained by taking derivatives of the generating function with respect to Ji. For example, . The cumulants can be obtained by taking the derivatives of W[J] = ln Z. Field theory generalizes the generating function over a set of discrete variables to a generating functional over functions or fields. The result is a functional or path integral over all possible paths of time evolution for the system, weighted by the probability of that particular evolution. While it is sometimes possible and desirable to take the spatial continuum limit of the neural system, this is not necessary for the path integral approach. Here we use continuum spatial notation, although the results carry through in the case where x indexes a discrete variable. Buice and Cowan (2007) showed that the generating functional for the master equation, 2.2, is given by
4.2
where
4.3
is called the action and we use the notation u · v = ∫ddxdt u(x, t)v(x, t). The fields φ and , which are obtained from the configuration variables , are defined below (Janssen & Tauber, 2005; Tauber et al., 2005). The asterisk denotes convolution of the weight function with the inputs, and the term is the cumulant generating functional of the initial distribution and takes into account arbitrary distributions in the initial condition. For example, if the initial state is described by Poisson statistics, we have , where a0(x) is the mean of the Poisson distribution at x. Analogous to the generating function for discrete variables, functional derivatives with respect to yield the normal-ordered cumulants, such as
4.4
and
4.5
Within this formalism (Doi, 1976a, 1976b; Peliti, 1985; Buice & Cowan, 2007), expectation values of products of φ are the normal-ordered cumulants found in the moment hierarchy. The normal-ordered cumulant C(x, y, t) from equation 3.9 results from setting t = t′ in C(x, t; x′, t′). The field is a “response” field, and expectation over functions of it yields Green's functions or response functions for the dynamics (Martin, Siggia, & Rose, 1973). The Ito convention is taken for the Langevin equation, 4.6, so that moments that involve combinations of and φ that are evaluated at the same time are zero. More specifically, the convention taken is such that if tt′.
We can heuristically derive the action (see equation 4.3), which was derived explicitly in Buice and Cowan (2007), and show that it represents a minimal model of the Wilson-Cowan equation where the activity is to be interpreted as a stochastic intensity or rate of a Poisson process. Consider an effective Wilson-Cowan Langevin equation,
4.6
where n(x,t) is the neural activity at location x and time t, ξ(x, t) is an effective stochastic forcing with probability density functional P[ξ], and the firing rate function F is not necessarily that which appears in equation 2.1, but rather is that in the master equation, 2.2. We will show that Poisson noise is necessary to match the Buice and Cowan action, equation 4.3. The probability density functional for n(x,t) can be written formally as
4.7
where δ[ · ] is the functional generalization of the Dirac delta function. The probability density 4.7 constrains the field n(x,t) to obey the Langevin equation, 4.6, with initial condition n0(x). The delta functional is defined by the generalized Fourier transform,
4.8
and is integrated along the imaginary axis. We can now integrate over the stochastic variable ξ to obtain a noise-generating functional defined by
4.9
We choose ξ such that
4.10
which is consistent with a Poisson activation at rate F and a Poisson decay at rate α. We next transform to the new variables:
4.11
The transformation 4.11 to the new fields serves to simplify the noise-generating functional and results in an action that has the form of equation 4.3, but with a different gain function. This new action is reconciled with equation 4.3 by the normal-ordering operation. The reason this is necessary is that the Ito convention used to interpret the action would not hold uniformly for the transformed fields since in performing the transformation, equation 4.11, the and φ fields inside the gain function are evaluated at the same time and moments between these particular instances of the fields would not necessarily be zero. This inconsistency can be corrected by redefining (i.e., normal ordering) the terms in the gain function. As an example, consider the firing rate function to be F(n) = (w · n)2. After transforming, it becomes , which mixes response and configuration operators at the same time point. To restore Ito convention, we normal-order so that response and configuration variables are no longer mixed. We do this by considering the n(x,t) operators at separate times t, t′ and computing how the operators approach each other as tt′. The properties of the response field provide and so that F becomes
4.12
Hence, restoring the Ito convention requires the replacement n2n2 + n (albeit in the variables and similarly for higher powers of n). This normal ordering will adjust the gain function to be , leaving us with the action, equation 4.3. From the perspective of the original master equation, 2.2, the transformation to is equivalent to expanding solutions to the master equation around a Poisson solution. Lefèvre and Biroli (2007) provide a lengthier discussion of the relationship between path integrals formulated in terms of the density and those formulated in terms of the deviations from Poisson behavior.

### 4.1.  Closed-Activity Equations from the Path Integral.

We derive the generalized activity equations for the normal-ordered cumulants directly from the generating functional by Legendre transforming to an effective action and then calculating the extrema of the effective action (Zinn-Justin, 2002; Cornwall, Jackiw, & Tomboulis, 1974; Buice & Cowan, 2007). We first perform the computation for the mean activity a(x,t) and then show how to generalize to arbitrary numbers of cumulants. The generating functional, equation 4.2, can be written more compactly as
4.13
where we define Φμ(x, t), where μ ∈ {1, − 1} such that Φ1(x, t) = φ(x, t) and . We use Einstein summation convention (i.e., when the same index appears twice—one upper, one lower—in equations, summation will be implied). Similarly we define Jμ(x, t) via J1 = J and . We can “raise” an index μ via multiplication with the matrix
4.14
so that and J−1 = J. Thus, we have . In order to streamline our notation, we will also define the dot product as the integral
4.15
For functions of more than one spatial variable, this inner product notation will generalize to the trace of the matrix product:
4.16
We will also write the action as Sμ] = ∫ddx dt Lμ], where L is the integrand of the action in equation 4.3.
Define aμ = 〈Φμ〉. The “effective action” Γ[aμ] for aμ is derived by a Legendre transformation,
4.17
where the conditions
4.18
4.19
are enforced. In analogy with classical mechanics, the extrema of the effective action,
4.20
give the equations of motion or the activity equations for aμ(x, t).
In general, we will not be able to compute the equation of motion exactly since the path integral in equation 4.13 cannot be computed exactly. However, we can perform a steepest-descent asymptotic expansion of equation 4.13 and compute the activity equation perturbatively. In field theory, this is known as the loop expansion, because the terms in the expansion can be represented by Feynman diagrams whose order is given by the number of loops that the diagram possesses. Substituting for W[Jμ] using the Legendre transformation, equation 4.17, gives
4.21
where we have suppressed the x and t arguments. Defining a new variable Ψμ = Φμaμ and using equation 4.19 gives
4.22
where we define Γμ[aμ] ≡ δΓ/δaμ. We now expand Sμ + aμ] in a functional Taylor series to obtain
4.23
where Lμ represents the functional derivative of L[aμ] with respect to aμ and similarly for Lμν and higher derivatives. Vμ, aμ] represents the remaining terms in the Taylor expansion. By definition, these are at least cubic in Ψμ. Hence,
4.24
We introduce a scaling parameter for the action, h (which we will set equal to one), via
4.25
We will show that an expansion in powers of h is consistent with the truncation used in deriving the moment hierarchy in section 3. The reason is that it organizes the terms in the expansion so that the true small parameters in the system, namely, inverse distance from criticality and inverse number of inputs, are manifested. We thus consider an asymptotic expansion Γ = Γ0 + hΓ1 + h2Γ2 + ⋅ ⋅ ⋅ If we set Γ0 = S, we obtain
4.26
Computing the corrections involves taking expectation values of the operator eV/h as well as other operators with respect to the gaussian functional with covariance (L−1)μν, which can be expanded as an infinite series of gaussian moments. Fortunately, we can describe the terms in this series graphically using Feynman diagrams. A result of this analysis is that we can arrange the corrections to the effective action according to the number of loops in the Feynman diagrams, the order in h being given by the number of loops (Zinn-Justin, 2002). S[aμ] provides the no loop- or “tree”-level computation. The lowest-order correction that these terms can produce is O(1), which would be an O(h) correction to the effective action (i.e., one loop correction). This is because the corrections will be given by moments of operators, which go as 1/h under a gaussian functional distribution whose covariance goes as h. The terms Γ1 and higher produce still higher-order corrections. We discuss in appendix A the connection between the h expansion, an artificial parameter, and the effective small parameters in the system (i.e., the inverse of the distance to criticality and the inverse of the extent of connectivity within the network, as addressed in section 3.1).
To lowest order, we obtain Γ[aμ] = S[aμ], which implies from equation 4.3 that the equations of motion to lowest order are given by
from which we obtain (because there is no initial condition “source” term for ) and the mean field Wilson-Cowan equation, 2.1. We can go to higher order by performing a loop expansion on the path integral in equation 4.26, and this next correction was computed in Buice and Cowan (2007).

The importance of this approach is that as we consider successive orders in the loop expansion, the effective action closes the system automatically. If we could calculate Γ[aμ] for our model of interest, then we would have the exact equation of motion for the true mean of the theory. In essence, we are trading a closure problem for an approximation problem. The advantage gained is that we do not have to worry about the contributions of higher moments explicitly and can consider explicitly the criteria allowing us to truncate the expansion. If there is an explicit loop expansion parameter, this truncation is straightforward. If not, as in our case, we must explicitly assess the regimes in which any truncation is valid. Even in the case where a truncation fails, the loop expansion can provide guidance in terms of identifying classes of diagrams (i.e., terms in the expansion) that are relevant in appropriate limits, which could be summed.

We can now generalize this procedure for equations of motion for an arbitrary number of cumulants by considering a generating functional for an arbitrary number of “composite operators” (Cornwall et al., 1974). In the case of the first and second cumulants, a(x,t) and C(x, y, t), we define the composite cumulant generating functional:
4.27
We now perform a double Legendre transform to obtain the effective composite action,
4.28
with conditions
4.29
4.30
and
4.31
4.32
The equations of motion are obtained by setting Jμ = 0 and Kμν = 0 in equations 4.31 and 4.32. We calculate the equations of motion to lowest order for this system in appendix B. The results are
4.33
and for Cμν, we get
4.34
4.35
together with the conditions
4.36
4.37
The two-point correlation designated as Cμν(x, t; x′, t′) generalizes the cumulant C(x, y, t) in equation 3.9 to include both the unequal time two-point correlation (C11) and the linear response (C1,−1). The equation for C(x, y, t) is obtained by taking the equation for , which results in equation 3.24. Note that we have also produced an equation for the Green's function of the theory C1,−1 as well as its time-reversed counterpart. Time reversal involves swapping the fields φ and . In the time-reversed theory, plays the role of activity. Time reversal does not give an equivalent theory since our Markov process is not time reversal invariant in general. In field theory language, this system of equations is known as the 2PI equations, and we adopt this moniker for convenience.1 We can then analogously derive nPI equations for any number of normal-ordered cumulants.

With the moment hierarchy approach, in order to produce better approximations, we are required to truncate further in the hierarchy. This can quickly produce unwieldy equations. The loop expansion provides an alternative in that corrections to the generalized equations can be produced with a diagrammatic expansion, namely, the one that calculates Γ2[aμ, Cμν], from which the corrections to the equations can be calculated.

## 5.  All-to-All Networks, Finite Size Effects, and Simulations

We consider the example of an all-to-all system, wherein each neuron connects to the entire network. Mean field theory should work well in this case because the number of postsynaptic neurons reduces the coupling of the fluctuations. In this case, the fluctuations reasonably reduce to corrections due to the finite size of the network, as we would expect. We take the weight function to be a constant, normalized by the number of neurons in the system wij = w0/N for some w0. The generalized Wilson-Cowan equations become
5.1
5.2
We can simplify this further by taking the initial conditions ai(0) = a0 and Cij(0) = 0 and assuming homogeneous external input I. This corresponds to initial conditions in the network determined by a Poisson distribution with mean a0. Cij(0) = −ai(0)δij would indicate starting with precisely ai(0) active neurons at i at time t = 0. Then symmetry reduces the equations to
5.3
5.4
where we have defined
5.5
5.6
Note that as N → ∞ the source term for C(t) vanishes, which implies that C(t) decouples from the equation for a(t), which then reduces to the standard Wilson-Cowan equation. The matrix Γij in this case is the function
5.7
The steady-state value of C(t) is given by
5.8
The relative size of the fluctuations in steady state is determined by the product NΓ[a0]. Large networks or networks distant from a bifurcation experience reduced correlations.
We now examine the phase plane of this simplified system. For concreteness, consider the firing rate function f(s) to be
5.9
where Θ is the Heaviside step function. At mean field level (i.e., consider C(t) to be zero) equilibria are determined by solutions of the equation
5.10
Figure 1 graphically displays the solutions of this equation. From the figure, we see that the equation exhibits a bifurcation as the value of α (the slope of the straight line in the figure) is decreased. The critical value for this bifurcation is α = 1.0. We will refer to the nonzero stable equilibrium as the “activated” state. The generalized activity equations, 5.3 and 5.5, also exhibit a bifurcation. The phase plane is shown in Figure 2 for α = 0.5, 0.9 and 1.0 with N = 100 and in Figure 3 for N = 10 with the same values of α. As expected, the steady-state value of C is larger for N = 10. Note that the nullclines for C(t) display divergences associated with Γ approaching zero. In addition, note that the location of the bifurcation is different. For N = 100, the bifurcation happens near α = 0.9, and with N = 10, the bifurcation happens for 0.9 < α < 1.0. Because the generalized equations are a coupled system, it is possible that a more interesting bifurcation structure may be manifested. In addition to the fixed points that exist in mean field theory, there is a new fixed point. Whereas the bifurcation at mean field level is a pitchfork bifurcation, that of the generalized equations is a saddle node bifurcation, with an unstable fixed point emerging as the a and C nullclines cross. There is always a fixed point at a = 0 and C = 0 because it is an absorbing state.
Figure 1:

Graphical depiction of solutions to equation 5.11 for various values of α.

Figure 1:

Graphical depiction of solutions to equation 5.11 for various values of α.

Figure 2:

Phase planes for the all-to-all generalized equations with (a) α = 0.5, (b) α = 0.9, and (c) α = 1.0. N = 100. Solid (black) lines are a nullclines; dotted lines are C nullclines.

Figure 2:

Phase planes for the all-to-all generalized equations with (a) α = 0.5, (b) α = 0.9, and (c) α = 1.0. N = 100. Solid (black) lines are a nullclines; dotted lines are C nullclines.

Figure 3:

Phase planes for the all-to-all generalized equations with α = 0.5 on the left, α = 0.9 in the center, and α = 1.0 on the right. N = 10. Solid (black) lines are a nullclines; dotted lines are C nullclines.

Figure 3:

Phase planes for the all-to-all generalized equations with α = 0.5 on the left, α = 0.9 in the center, and α = 1.0 on the right. N = 10. Solid (black) lines are a nullclines; dotted lines are C nullclines.

Importantly, we see that we can alter the bifurcation structure by adding a forcing or source term to the correlation function C(t) equation, linearly shifting the C nullcline. This removes the stable fixed point for high a (the one associated with the activated state in mean field theory). Because of this, we see that we can disrupt the activated state by stimulating the system with an input such that the correlation is sourced more strongly than the mean. We can use this to “turn off” the activated state by synchronizing the network. These correlations drive the system to the absorbing state of the full model. To reverse this deactivation, we simply drive the system with Poisson noise (i.e., force the equation for a(t) but not C(t)). Compare this to the effect demonstrated in Laing and Chow (2001) in which synchronized activity associated with fast synapses led to the destabilization of activity that the Wilson-Cowan equation predicts to be stable. For a saturating firing rate function (more generally, a function such that f′(s) < 0 in the appropriate region), increased correlations inhibit the mean activity ai(t).

We now demonstrate the utility of the generalized activity equations, 3.21 and 3.22, for describing the full Markov system 2.2 away from a bifurcation point. In order to simulate the Markov model, we use a form of the Gillespie algorithm and take expectation values over many time evolutions of the system. We use the firing rate function, equation 5.10. It is important to remember when comparing results with simulations that we use the whose normal ordered form is f(s). In this all-to-all case, we need only consider the correction arising from the quadratic portion of the firing rate function since the corrections will go as powers of the weight function, which in this all-to-all case means they carry factors of 1/N. In particular, to lowest order we have
5.11
We plot ai(t) and Cii(t) versus t for various values of α and N in Figures 4 through 7. (Note that we are numerically evaluating the generalized equations, not the simplified equations 5.3 and 5.5, and comparing them with Monte Carlo simulations of the Markov system; the plots overlay data for each of the N neurons.) We initialize the network with Poisson initial conditions: ai(0) = 2 and Cij(0) = 0. The simulations of the Markov process are averaged over 105 instances. One can see that the equations match the simulations quite well away from the critical point. As one approaches the bifurcation, however, the simulations begin to deviate. At α = 0.5, mean field and the generalized equations each match the simulated Markov process. As one approaches the mean field bifurcation point of α = 0.9, the mean field equations no longer match well, but the generalized equations account for the deviation. From Figure 5, one can see that as we approach α = 1.0, the estimate of the correlation from the generalized equations becomes poorer.
Figure 4:

a(t) versus t for (a) α = 0.5, (b) α = 0.9, and (c) α = 1.0. N = 100. Dotted lines are solutions of mean field theory. Dashed lines are solutions of the generalized equations. Solid lines are expectations values of data from simulations of the Markov process.

Figure 4:

a(t) versus t for (a) α = 0.5, (b) α = 0.9, and (c) α = 1.0. N = 100. Dotted lines are solutions of mean field theory. Dashed lines are solutions of the generalized equations. Solid lines are expectations values of data from simulations of the Markov process.

Figure 5:

C(t) vs. t for (a) α = 0.5, (b) α = 0.9, and (c) α = 1.0. N = 100. Dashed lines are solutions of the generalized equations. Points (black) are expectations values of data from simulations of the Markov process.

Figure 5:

C(t) vs. t for (a) α = 0.5, (b) α = 0.9, and (c) α = 1.0. N = 100. Dashed lines are solutions of the generalized equations. Points (black) are expectations values of data from simulations of the Markov process.

The plots for N = 10 in Figures 6 and 7 demonstrate the breakdown of the generalized equations. There is already a significant deviation of both mean field and the generalized equations at α = 0.5. Naturally the discrepancy is accounted for by the poor estimation of the correlation at this level. As we near α = 1.0, the estimate of the correlations begins to grow, whereas the simulated correlation is dropping to zero.

Figure 6:

a(t) versus t for (a) α = 0.5, (b) α = 0.9, and (c) α = 1.0. N = 10. Dotted lines are solutions of mean field theory. Dashed lines are solutions of the generalized equations. Solid lines are expectations values of data from simulations of the Markov process.

Figure 6:

a(t) versus t for (a) α = 0.5, (b) α = 0.9, and (c) α = 1.0. N = 10. Dotted lines are solutions of mean field theory. Dashed lines are solutions of the generalized equations. Solid lines are expectations values of data from simulations of the Markov process.

Figure 7:

C(t) versus t for (a) α = 0.5, (b) α = 0.9, and (c) α = 1.0. N = 10. Dashed lines are solutions of the generalized equations. Solid lines are expectations values of data from simulations of the Markov process.

Figure 7:

C(t) versus t for (a) α = 0.5, (b) α = 0.9, and (c) α = 1.0. N = 10. Dashed lines are solutions of the generalized equations. Solid lines are expectations values of data from simulations of the Markov process.

One can see that although the theory begins to deviate from the simulations near criticality, we still capture the loss of stability of the active state, even for N = 10. This is due to the growth of correlations that negatively feed back on the mean due to the negative second derivative of f(s). We can use this feature to observe the effect of correlated input directly by adding a term to the Markov process that provides multiplicative gaussian noise. In particular, we add a transition at rate
5.12
where η is a gaussian noise source. The purpose of the step function is to prevent an individual neuron from getting a kick that will drive the activity negative. Note that we are not using an “input” term as we have defined it. Because the firing rate function is strictly positive, we cannot use a stimulus such that the mean transition rate is strictly zero. However, we can use a stimulus such that the source to the correlation function is much stronger than the mean. For our example, to maximize the effect, we do not use an input to the firing rate function but instead add another transition to the Markov process separate from the firing rate function and the decay transition. This allows us to source only the correlation function. Although this is artificial, this transition adds terms to the coupled equations that would approximate those from some input with zero mean and large correlations. In the following simulations, we used α = 0.5, σ = 100, and N = 100. The correlated input was given between t = 20 and t = 30. In the absence of external correlated input, these parameters result in the active state being stable, as shown in Figure 4. As can be seen in Figure 8, the use of the global noise source results in the “switch-off” behavior predicted by the generalized equations. If we instead use a Poisson process to provide this external stimulation, one can also see in Figure 4 that the network responds and then reverts back to the active state. The reason for the explicit shut-off is that the system has an absorbing state to which it is driven. The more important point is that the correlated input is acting as a source of inhibition, whereas the Poisson input serves as an excitatory input. A linear system or one in which the firing rate function f(s) is such that f′(s)>0 will not exhibit this behavior. We chose this particular example for sourcing Cij(t) in order to separate the effects of sourcing ai(t) as well. Given a more complicated noise source, one would need to examine the phase plane or solve the equations after determining the effects of the noise source on the normal-ordered cumulants.
Figure 8:

Response of all-to-all network to correlated input. α = 0.5, N = 100. (a) The response to correlated input with σ = 100. (b) The response to a Poisson process with rate λ = 10. Note the change in scale between the two plots.

Figure 8:

Response of all-to-all network to correlated input. α = 0.5, N = 100. (a) The response to correlated input with σ = 100. (b) The response to a Poisson process with rate λ = 10. Note the change in scale between the two plots.

## 6.  Discussion

We have demonstrated a formalism for constructing a minimal extension of a rate equation to include correlated activity. We have explicitly computed the lowest-order results of this extension, which provide coupled equations for the mean activity, two-point correlations, and linear responses. These results indicate that correlations can have an important impact on the dynamics of a rate equation, affecting both stability and the structure of bifurcations. Our argument relied on inferring a “minimal” Markov process. Our use of the Doi-Peliti path integral formalism guides our assertion that our inferred Markov process is the simplest one compatible with both the rate equations and their interpretation as measuring some stochastic counting process. Thus, a general extension for any type of rate equation should share the same basic structure that we have described here. We performed this construction on a Markov process consistent with the Wilson-Cowan equation, but our prescription would work equally well with any Markov process.

In keeping with this idea, our results have something in common with other approaches to understanding correlations in neural networks. El Boustani and Destexhe (2009) attempt to derive a Markov model for the asynchronous irregular states of an underlying neural system and explore the moment hierarchy of that Markov model. We take the opposite approach, beginning with a presumed set of rate equations and asking what possible restrictions can be placed on the correlation functions knowing only the dynamics of the rate model. Their hierarchy is truncated via scaling and finite size, whereas our hierarchy's truncation (and the truncation of the loop expansion) arises through the distance to a bifurcation in the rate equations. Ginzburg and Sompolinsky (1994) propose a simple Markov model and study its moment hierarchy. For the correlations, they achieve results similar to the tree level of our loop expansion. An important point of departure is that we consider the recurrent effects the correlations have on the mean field, which we demonstrate can be sufficiently significant to alter the structure of the bifurcation.

As we predict, our theory breaks down sufficiently close to a bifurcation. Examining the dynamics near the critical point requires a different form of analysis such as a renormalization argument. An example was presented in Buice and Cowan (2007), where it was argued that a large class of neural models would exhibit scaling laws near a bifurcation like those of the directed percolation model (Janssen & Tauber, 2005). The predictions of this scaling coincide with measurements made in cortical slices of “avalanches” (Beggs & Plenz, 2003). If criticality is important for neural function (Beggs, 2008), then these sorts of approaches will be more important for future work, and our rate model extension will be less applicable.

In contrast, supporting the potential usefulness of our rate model extension is the fact that large neural connectivity suppresses correlations and aids the truncation of the hierarchy, an analogous result to the Ginsburg criterion in equilibrium statistical mechanics. In addition, we demonstrated that Poisson-like input in general pushes the system into configurations in which the correlations are suppressed relative to the mean. All of this suggests that our extension will be at least as applicable as the rate models themselves.

Regarding that applicability, both the Markov process and the rate equations assume a large degree of underlying asynchrony in the network. The expansion we describe should be appropriate for networks in which a relatively small amount of synchrony at the level of individual neurons is developing. The coupled correlation function captures this effect. If the population is being dominated by neuron-level synchrony, then the Markov process should no longer hold as a description of the system. Population-level synchrony as captured by the original rate model, however, should have no effect on our analysis. In other words, there will be correlation-effects acting on oscillating states, for example, such as presumably correlation-induced modulation of the frequency of the oscillation. We will demonstrate this in future work.

An important outstanding point is that we have posited this Markov process based on the original interpretation of the Wilson-Cowan equations as dynamical equations for the probabilistic activity of neurons. Although our Markov process is the most “natural” given the transitions in the Wilson-Cowan equations, there is no a priori reason to suppose that this Markov process reflects the probabilistic dynamics of a physiologically based neural model or of real neurons precisely because there is nothing that mandates this interpretation. Per the renormalization analysis of Buice and Cowan (2007), measuring scaling laws in cortex will provide a means of identifying classes of models (by identifying the relevant universality class), but this will in no way distinguish among models within the same class. Distinguishing models within the same class will require the measurement of nonuniversal quantities. This would likely mean some relatively precision measurements of response functions in cortical activity.

Our model is one in which activity is Poisson distributed in the absence of connectivity. Likewise, we have shown that connectivity, away from a bifurcation, only slightly modifies this. Although we introduced the model in terms of populations, its original description arose in the context of the low firing rate reduction of the Cowan two and three state neural Markov processes (Cowan, 1991), where ai(t) is the probability of neuron i being active at time t. The unboundedness of the variable ai(t) arises from the low firing rate assumption. In Buice and Cowan (2007), the model is described as counting the number of recent spikes from a single neuron at i that still produce an effect upon the activity of the network; for a given spike, this effect decays with rate constant α. Either of these descriptions provides justification for the Poisson nature of solutions in terms of single neuron activity. If we have large numbers of equivalent neurons, then an averaging process reduces the effect of correlations, consistent with what we have shown here. Our all to all example produces a network in which the deviations from Poisson are governed by finite size effects. A population average over such equivalent neurons will reduce the fluctuations according to the size of the population. In the limit of an infinite number of identical neurons, we are left with mean field theory.

We feel our approach is a useful starting point for understanding effects beyond mean field. We have demonstrated a correlation-induced loss of stability in an all-to-all network. This effect should carry over to nonhomogeneous solutions such as bump solutions or traveling waves. Likewise, correlations will modify important aspects of mean field solutions such as dispersion relations. Our approach enables this dispersion relation to be calculated. In addition to stability, our equations are a useful starting point toward understanding the wandering of bump solutions. They also provide a natural means of studying beyond mean field effects of correlation-based learning. A model of visual hallucinations in cortex based on the Turing mechanism has explained many hallucinatory effects (such as the various Kluver form constants). Since the Turing mechanism is based on bifurcations, it is an interesting question to what extent the coupling with correlations affects the results of the hallucination analysis. Our approach may provide this model with a means of explaining further hallucinations not covered by the model in Bressloff et al. (2001).

It remains an important question how to connect our Markov and generalized rate model approach with models of deterministic neurons. While the formalism admits almost any gain function, there remains the question of connecting this gain function to, for example, the transfer function for some neural model of which the Markov process is some approximation. This is, of course, not a question of the analysis of Markov models but of the applicability of rate models as high-level descriptions of more detailed neural models. Answering this question will likely involve a kinetic theory formulation of the neural systems, such as the one pursued in Hildebrand et al. (2006) and Buice and Chow (2007). Having derived the generalized equations, it is also now important to explore their further consequences for phenomena such as pattern formation. There are also many avenues to extend this model and this approach. The Markov process can be enlarged to account for synaptic adaptation by adding a synaptic time variable to the neural configuration. Likewise, noise in the transitions themselves, whether spatial or temporal, is easily incorporated into the action. A reduction of the resulting theory would no longer satisfy the Markov property, although there may be certain natural assumption (such as slow dynamics for the auxiliary field) that could allow one to regain Markovicity with an approximate model. This would allow us to construct extended Wilson-Cowan equations that incorporate these and other aspects of neural dynamics. These questions will be explored in future work.

## Appendix A:  Composite Operator Effective Action and the 2PI Equations

Here we derive the 2PI equations. We begin with the generating functional:
A.1
where we have introduced a parameter h = 1 for bookkeeping purposes. The generalized effective action is given by
A.2
with Jμ and Kμν given by
A.3
A.4
The path integral representation is thus
A.5
which we transform to a new variable Ψμ(x, t) = Φμ(x, t) − aμ(x, t), set
A.6
and expand to obtain
where we use
We consider the expansion Γ = Γ0 + hΓ1 + h2Γ2, where Γ2 contains all terms of order h2 and higher. Setting Γ0 = S gives
We now fix Kμν according to
A.7
which gives
A.8
where
A.9
We need to extract the order h contributions from the functional integral. We expect the effect of Γ0,μν is to replace Lμν with the full inverse two-point function (C−1)μν. This in turn will affect the normalization of the integral. Because of this, we expect the order h contribution to the effective action to be
A.10
Substituting this into expression A.10 gives us
A.11
We can extract the normalization of the functional integral using and the identity to obtain
A.12
The factor of the determinant serves as a normalization for the functional integral, which is now in a form that will contribute to the effective action only at order h2. Thus, we have
A.13
Now we can calculate the equations of motion to a given loop order from equations 4.31 and 4.32. Using equation 4.31, the equations of motion for the mean field are
A.14
The equations for Cμν are
A.15
which we can invert to get
A.16
In particular, if we ignore loop corrections (consider only first order in h, recalling that Γ2 is O(h2)), we get
A.17
where L−1μν[aμ] is the inverse of Lμν[aμ]. In the absence of interactions, L−1μν is the two-point function, as expected.

We can use the loop expansion to draw some conclusions about the applicability of perturbation theory. Since Γ2[aμ, Cμν] is second order and is the sum of vacuum two-particle irreducible graphs, every graph contributing to it must be at least of two loop order. Every internal line represents a factor of Cμν, and so each graph contributing to Γ2[aμ, Cμν] must have at least two factors of Cμν, each of which will be either equal to 0 or be attenuated (in steady state) by the same exponents that attenuate the magnitude of C(x, y, t) away from a bifurcation, according to equations B.14 to 4.36. Thus, the argument that Cμν is small away from the critical point extends to every term in the expansion for the generalized equations.

The caveat here is that there is a class of diagrams that couple the lowest-order expression for a given moment to the mean field. Although these graphs are suppressed by the distance to criticality, each of these is of the same order. We are assisted by two facts. The first is that the source terms for each of these moments at lowest order will be proportional to derivatives of the firing rate function. If f(s) is sufficiently smooth, this will suppress higher-order contributions. In addition, each coupling will go as an additional factor of N−1m where Nm was defined in section 3.1 as the smallest number of inputs to any given neuron. Thus, the connectivity in cortex will serve to “average-out” sources to the mean from higher moments. This will be the case as long as we can bound the total input to any given neuron.

## Appendix B:  Tree Level Equations of Motion

In order to calculate the expansion for the equations of motion, we need to compute the value of both Lμν and Γ2. We compute the lowest-order correction here.

First, we find the intermediate results (which give us the classical equations of motion for a and ):
B.1
B.2
from which follows
B.3
B.4
B.5
B.6
The terms f(n)(x, t) indicate the nth derivative of f. Note that we have suppressed the argument, so that
We can now write down the equations of motion from equation A.16, minus the loop corrections. The first “diagonal” equation (for (−1, − 1)) is
B.7
The second “diagonal” equation (for equation 3.6) is
B.8
The “off-diagonal” equations are (starting with −1, 1):
B.9
and the other (1, − 1):
B.10

The “mean field” portion of the equations of motion, A.16, are obtained from equations B.1 and B.2 (by setting the left-hand side to zero). The remainder of the equations of motion are “classical” terms dependent on the correlation functions and loop corrections. The latter are given by the term in equation A.16. The term in the trace is, of course, the sum of the left-hand side of equations B.7 and B.8.

We can simplify the equations for the mean field by realizing that any term involving C−1,1 or C1,−1 can be ignored because they will appear only in the form C−1,1(x′, t; x, t), that is, at equal initial and final times. These will be zero. This can be seen as either the “initial condition” for the linear response terms or as a manifestation of the “ϵ(0)” problem in quantum field theory (Zinn-Justin, 2002).

Furthermore, we can use some results from the full theory. In particular, we have
B.11
B.12
because causality enforces that operators cannot contract with anything “in the past.”
The equation for a(x,t) is then:
B.13
Applying these same simplifications to the equations for Cμν, we get:
B.14
B.15
B.16
together with the conditions
B.17
B.18

## Acknowledgments

We thank Vipul Periwal for helpful suggestions. This work was supported by the Intramural Research Program of the NIH/NIDDK.

## Notes

1

“2PI” stands for 2 particle irreducible. The effective action Γ[aμ] is the generating functional of 1PI graphs, that means that the graphs that determine Γ[aμ] cannot be disconnected by cutting any single line of the graph. Similarly, Γ[aμ, Cμν] is 2PI in the sense that graphs contributing to it cannot be disconnected by cutting two lines.

## References

Abbott
,
L.
, &
Van Vreeswijk
,
C.
(
1993
).
Asynchronous states in networks of pulse-coupled oscillators
.
Phys. Rev. E
,
48
(
2
),
1483
1490
.
Amari
,
S.
(
1975
).
Homogeneous nets of neuron-like elements
.
Biological Cybernetics
,
17
,
211
220
.
Amari
,
S.
(
1977
).
Dynamics of pattern formation in lateral-inhibition type neural fields
.
Biological Cybernetics
,
27
,
77
87
.
Amit
,
D.
, &
Brunel
,
N.
(
1997a
).
Model of global spontaneous activity and local structured activity during delay periods in the cerebral cortex
.
Cereb. Cortex
,
7
(
3
),
237
252
.
Amit
,
D. J.
, &
Brunel
,
N.
(
1997b
).
Dynamics of a recurrent network of spiking neurons before and following learning
.
Network: Computation in Neural Systems
,
8
(
4
),
373
404
.
Amit
,
D. J.
,
Gutfreund
,
H.
, &
Sompolinsky
,
H.
(
1985
).
Spin-glass models of neural networks
.
Physical Review A
,
32
(
2
),
1007
1018
.
Beggs
,
J. M.
(
2008
).
The criticality hypothesis: How local cortical networks might optimize information processing
.
Philosophical Transactions of the Royal Society A
,
366
(
1864
),
329
343
.
Beggs
,
J. M.
,
Plenz
,
D.
(
2003
).
Neuronal avalanches in neocortical circuits
.
Journal of Neuroscience
,
23
(
35
),
11167
11177
.
Beshel
,
J.
,
Kopell
,
N.
, and
Kay
,
L.
(
2007
).
Olfactory bulb gamma oscillations are enhanced with task demands
.
Journal of Neuroscience
,
27
(
31
),
8358
8365
.
Bressloff
,
P.
,
Cowan
,
J. D.
,
Golubitsky
,
M.
,
Thomas
,
P. J.
, &
Wiener
,
M. C.
(
2001
).
Geometric visual hallucinations, Euclidean symmetry and the functional architecture of striate cortex
.
Philosophical Transactions of the Royal Society of London B
,
356
,
299
330
.
Brunel
,
N.
(
6 May 2000
).
Dynamics of sparsely connected networks of excitatory and inhibitory spiking neurons
.
Journal of Computational Neuroscience
,
8
,
183
208
.
Brunel
,
N.
, &
Hakim
,
V.
(
1999
).
Fast global oscillations in networks of integrate-and-fire neurons with low firing rates
.
Neural Computation
,
11
(
7
),
1621
1671
.
Brunel
,
N.
, &
Hansel
,
D.
(
2006
).
How noise affects the synchronization properties of recurrent networks of inhibitory neurons
.
Neural Comput.
,
18
(
5
),
1066
1110
.
Buice
,
M. A.
, &
Chow
,
C. C.
(
2007
).
Correlations, fluctuations, and stability of a finite-size network of coupled oscillators
.
Physical Review E (Statistical, Nonlinear, and Soft Matter Physics)
,
76
(
3
),
031118
.
Buice
,
M. A.
, &
Cowan
,
J. D.
(
2007
).
Field theoretic approach to fluctuation effects for neural networks
.
Physical Review E
,
75
,
051919
.
Cai
,
D.
,
Tao
,
L.
,
Shelley
,
M.
, &
McLaughlin
,
D. W.
(
2004
).
An effective kinetic representation of fluctuation-driven neuronal networks with application to simple and complex cells in visual cortex
.
Proceedings of the National Academy of Sciences of the United States of America
,
101
(
20
),
7757
7762
.
Cohen
,
M.
, &
Grossberg
,
S.
(
1983
).
Absolute stability and global pattern formation and parallel storage by competitive neural networks
.
IEEE Transactions on Systems, Man, and Cybernetics
,
13
,
815
826
.
Coombes
,
S.
(
2005
).
Waves, bumps, & patterns in neural field theories
.
Biological Cybernetics
,
93
,
91
108
.
Cornwall
,
J. M.
,
Jackiw
,
R.
, &
Tomboulis
,
E.
(
1974
).
Effective action for composite operators
.
Physical Review D
,
10
(
8
),
2428
2445
.
Cowan
,
J. D.
(
1991
).
Stochastic neurodynamics
. In
D. S. Touretzky, R. P. Lippman, & J. E. Moody
(Eds.),
Advances in neural information processing systems
,
3
.
San Francisco
:
Morgan Kaufmann
.
Doi
,
M.
(
1976a
).
Second quantization representation for classical many-particle system
.
Journal of Physics A: Mathematical and General
,
9
(
9
),
1465
1477
.
Doi
,
M.
(
1976b
).
Stochastic theory of diffusion controlled reaction
.
Journal of Physics A: Mathematical and General
,
9
(
9
),
1479
1495
.
El Boustani
,
S.
, &
Destexhe
,
A.
(
2009
).
A master equation formalism for macroscopic modeling of asynchronous irregular activity states
.
Neural Computation
,
21
,
46
100
.
Ermentrout
,
B.
(
1998
).
Neural networks as spatio-temporal pattern-forming systems
.
Reports on Progress in Physics
,
61
(
4
),
353
430
.
Ermentrout
,
G. B.
, &
Cowan
,
J. D.
(
1979
).
A mathematical theory of visual hallucination patterns
.
Biological Cybernetics
,
34
,
137
150
.
Fourcaud
,
N.
, &
Brunel
,
N.
(
2002
).
Dynamics of the firing probability of noisy integrate-and-fire neurons
.
Neural Comput.
,
14
(
9
),
2057
2110
.
Fusi
,
S.
, &
Mattia
,
M.
(
1999
).
Collective behavior of networks with linear (VLSI) integrate-and-fire neurons
.
Neural Computation
,
11
(
3
),
633
652
.
Gerstner
,
W.
(
1995
).
Time structure of the activity in neural network models
.
Phys. Rev. E
,
51
(
1
),
738
758
.
Gerstner
,
W.
(
2000
).
Population dynamics of spiking neurons: Fast transients, asynchronous states, and locking
.
Neural Computation
,
12
,
43
89
.
Gerstner
,
W.
, &
Kistler
,
W. M.
(
2002
).
Spiking neuron models
.
Cambridge
:
Cambridge University Press
.
Ginzburg
,
I.
, &
Sompolinsky
,
H.
(
1994
).
Theory of correlations in stochastic neural networks
.
Physical Review E
,
50
(
4
),
3171
3191
.
Golomb
,
D.
, &
Hansel
,
D.
(
2000
).
The number of synaptic inputs and the synchrony of large, sparse neuronal networks
.
Neural Comput.
,
12
(
5
),
1095
1139
.
Gray
,
C. M.
, &
Singer
,
W.
(
1989
).
Stimulus-specific neuronal oscillations in orientation columns of cat visual cortex
.
Proceedings of the National Academy of Sciences of the United States of America
,
86
(
5
),
1698
1702
.
Hildebrand
,
E. J.
,
Buice
,
M. A.
, &
Chow
,
C. C.
(
2006
).
Kinetic theory of coupled oscillators
.
Physical Review Letters
,
98
,
1
4
.
Hopfield
,
J.
(
1982
).
Neural networks and physical systems with emergent collective computational abilities
.
Proceedings of the National Academy of Sciences
,
79
,
2445
2558
.
Hopfield
,
J.
(
1984
).
Neurons with graded response have collective computational properties like those of two state neurons
.
Proceedings of the National Academy of Sciences
,
81
,
3088
3092
.
Ichimaru
,
S.
(
1973
).
Basic principles of plasma physics: A Statistical Approach
.
Reading, MA: W. A. Benjamin Advanced Book Program
.
Janssen
,
H.-K.
, &
Tauber
,
U. C.
(
2005
).
The field theory approach to percolation processes
.
Annals of Physics
,
315
(
1
),
147
192
.
Laing
,
C.
, &
Chow
,
C.
(
2001
).
Stationary bumps in networks of spiking neurons
.
Neural Computation
,
13
(
7
),
1473
1494
.
Lefèvre
,
A.
, &
Biroli
,
G.
(
2007
).
Dynamics of interacting particle systems: Stochastic process and field theory
.
J. Stat. Mech.
,
2007
(
7
),
P07024
P07024
.
Ly
,
C.
, &
Tranchina
,
D.
(
2007
).
Critical analysis of dimension reduction by a moment closure method in a population density approach to neural network modeling
.
Neural Computation
,
19
(
8
),
2032
2092
.
Martin
,
P. C.
,
Siggia
,
E. D.
, &
Rose
,
H. A.
(
1973
).
Statistical dynamics of classical systems
.
Physical Review A
,
8
(
1
),
423
437
.
Mattia
,
M.
,
Del Giudice
,
P.
(
2002
).
Population dynamics of interacting spiking neurons
.
Phys. Rev. E
,
66
(
5
),
051917
.
Nicholson
,
D. R.
(
1992
).
Introduction to plasma theory
.
Malabar, FL
:
Krieger
.
Nykamp
,
D. Q.
, &
Tranchina
,
D.
(
2000
).
A population density approach that facilitates large-scale modeling of neural networks: Analysis and an application to orientation tuning
.
J. Comp. Neurosci.
,
8
,
19
50
.
Ohira
,
T.
, &
Cowan
,
J.
(
1993
).
Master equation approach to stochastic neurodynamics
.
Phys. Rev. E
,
48
,
2259
2266
.
Peliti
,
L.
(
1985
).
Path integral approach to birth-death processes on a lattice
.
Journal de Physique
,
46
,
1469
1483
.
Peretto
,
P.
(
1984
).
Collective properties of neural networks: A statistical physics approach
.
Biological Cybernetics
,
50
,
51
62
.
Plesser
,
H. E.
, &
Gerstner
,
W. E.
(
2000
).
Noise in integrate-and-fire neurons: From stochastic input to escape rates
.
Neural Comput.
,
12
(
2
),
367
384
.
Salinas
,
E.
, &
Sejnowski
,
T. J.
(
2002
).
Integrate-and-fire neurons driven by correlated stochastic input
.
Neural Comput.
,
14
(
9
),
2111
2155
.
Softky
,
W. R.
, &
Koch
,
C.
(
1993
).
The highly irregular firing of cortical cells is inconsistent with temporal integration of random EPSPS
.
Journal of Neuroscience
,
13
(
1
),
334
350
.
Soula
,
H.
,
Beslon
,
G.
, &
Mazet
,
O.
(
2006
).
Spontaneous dynamics of asymmetric random recurrent spiking neural networks
.
Neural Comput.
,
18
(
1
),
60
79
.
Soula
,
H.
, &
Chow
,
C. C.
(
2007
).
Stochastic dynamics of a finite-size spiking neural network
.
Neural Comput.
,
19
(
12
),
3262
3292
.
Tauber
,
U. C.
,
Howard
,
M.
, &
Vollmayr-Lee
,
B. P.
(
2005
).
Applications of field-theoretic renormalization group methods to reaction-diffusion problems
.
Journal of Physics A: Mathematical and General
,
38
(
17
),
R79
R131
.
Treves
,
A.
(
1993
).
Mean-field analysis of neuronal spike dynamics
.
Network: Computation in Neural Systems
,
4
,
259
284
.
Van Vreeswijk
,
C.
, &
Sompolinsky
,
H.
(
1996
).
Chaos in neuronal networks with balanced excitatory and inhibitory activity
.
Science
,
274
,
1724
1726
.
Van Vreeswijk
,
C.
, &
Sompolinsky
,
H.
(
1998
).
Chaotic balanced state in a model of cortical circuits
.
Neural Computation
,
10
,
1321
1371
.
Wilson
,
H. R.
, &
Cowan
,
J. D.
(
1972
).
Excitatory and inhibitory interactions in localized populations of model neurons
.
Biophysical Journal
,
12
,
1
24
.
Wilson
,
H. R.
, &
Cowan
,
J. D.
(
1973
).
A mathematical theory of the functional dynamics of cortical and thalamic nervous tissue
.
Kybernetik
,
13
,
55
80
.
Zinn-Justin
,
J.
(
2002
).
Quantum field theory and critical phenomena
(4th ed.).
New York
:
Oxford Science Publications
.