## 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

*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.

*t*, and

*n*is the number of active neurons at location

_{i}*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

*i*th component is

*n*± 1, respectively. The rate at which a neuron at location

_{i}*i*becomes active is given by the firing rate or gain function , an implicit function of the weight function

*w*and external inputs

_{ij}*I*. One of the crucial elements of the ensuing calculation is making the connection between the gain function , which appears in equation 2.2, and

_{i}*f*(

*I*(

_{i}*t*) + ∑

_{j}

*w*), which appears in equation 2.1.

_{ij}n_{j}*n*(

_{i}*t*)

*n*(

_{j}*t*′)

*n*(

_{k}*t*′

^{}′) ⋅ ⋅ ⋅ 〉 where the expectation value is over all statistical realizations of the Markov process. The first moment,

is a measure of the mean activity in the network. We obtain an equation for *a _{i}*(

*t*) by multiplying equation 2.2 by

*n*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

_{i}*n*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.

_{i}n_{j}*a*(

_{i}*t*). This is the naive mean field assumption where all cumulants are zero. For example, the second cumulant (i.e., variance) 〈

*n*(

_{i}*t*)

*n*(

_{j}*t*′)〉 −

*a*(

_{i}*t*)

*a*(

_{j}*t*′) is set to zero. For our master equation, 2.2, this assumption yields (see the computation in the next section) 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.

*a*(

_{i}*t*) can be written as 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

*F*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.

_{i}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

*a*(

_{i}*t*), we multiply the master equation 2.2 by

*n*and sum over all configurations : The first two terms on the right-hand side simplify to The first equality results from reindexing the summation over

_{i}*n*from (0, ∞) to (1, ∞). We leave the summation indicated as over all configurations because the factor of

_{k}*n*prevents the 0 term from contributing. We have also separated out the terms where

_{k}*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

*a*(

_{i}*t*). The next two are given by and The normal-ordered cumulants can be computed using a recursive algorithm. The algorithm involves replacing all moments

*N*

_{ijk⋅⋅⋅}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

*N*−

_{ij}*a*. To compute the normal-ordered version, we replace the moments appearing in the expression with the normal-ordered-corrected forms, that is, set

_{i}a_{j}*N*→

_{ij}*N*−

_{ij}*a*δ

_{i}_{ij}. The term subtracted results from the contraction of the

*i*,

*j*indices:

*N*→

_{ij}*N*δ

_{i}*ij*=

*a*δ

_{i}_{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 To obtain the normal-ordered cumulant, we first make the replacement

*N*→

_{ijk}*N*−

_{ijk}*N*δ

_{ij}_{jk}−

*N*δ

_{ik}_{ij}−

*N*δ

_{jk}_{ik}−

*a*δ

_{i}_{ij}δ

_{jk}. We then must correct for all appearances of the second moment

*N*, 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

_{lm}*C*

_{ijk⋅⋅⋅}= 0 except the first,

*a*. An alternative rationale for normal ordering will appear in section 4.

_{i}*N*

_{ijk⋅⋅⋅}, as given by the Taylor expansion of , We have implicitly defined the notation that

*F*

^{jk⋅⋅⋅}

_{i}is the derivative of

*F*with respect to

_{i}*n*,

_{j}*n*, …. 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

_{k}*N*

_{ijk⋅⋅⋅}will contribute a term proportional to

*a*⋅ ⋅⋅ to the normal-ordered cumulant expansion. This means that we can write the expansion in the form 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

_{i}a_{j}a_{k}*a*. These corrections form the series where by

_{i}*F*

^{jm}we mean the

*m*th derivative of with respect to

*n*. 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

_{j}*a*,

_{i}a_{j}*a*, ….

_{i}a_{j}a_{k}*f*(

*s*), where

_{i}*s*= ∑

_{i}_{i}

*w*+

_{ij}a_{j}*I*(

_{i}*t*), we assume that this resummation produces 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.

*C*(

_{ij}*t*). At each order

*m*in the series (i.e., the order that contains the

*m*th moment), there are

*m*(

*m*− 1)/2 terms that have one factor of

*C*(

_{ij}*t*). These terms sum to give which can be rewritten as Other terms are at least second order in

*C*, 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

_{ij}*C*(

_{ij}*t*) and obtain Equations 3.21 and 3.22 constitute a closed set of equations for the mean and variance of a Wilson-Cowan network of neurons.

*a*(

_{i}*t*) represents the Poisson rate, and

*C*(

_{ij}*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 *C _{ij}*(

*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

*C*(

_{ij}*t*) will receive input only when the activity is precisely near threshold. Also notice that the strength of the input to

*C*(

_{ij}*t*) is proportional to the weight

*w*between the neurons in question, as well as the mean activity. An initial check on the equations is that

_{ij}*C*decouples from

_{ij}*a*in the case where

_{i}*f*(

*s*) is linear or constant.

*a*(

*x*,

*t*) and correlation

*C*(

*x*,

*y*,

*t*): and 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 *C _{ij}*(

*t*) instead of allowing higher moments to interact with the mean activity.

*a*

^{0}

_{i}to be some steady-state solution to this equation and linearize equation 2.4 around this solution. The perturbations δ

*a*(

_{i}*t*), from this steady-state solution obey the equation We rewrite this equation as where the matrix Γ

_{ij}is defined by If all of the eigenvalues of Γ

_{ij}are positive, then the solution

*a*

^{0}

_{i}is stable. Likewise, negative eigenvalues indicate instability. Criticality is the condition of marginal stability, in which one or more of the eigenvalues are 0.

*C*(

_{ij}*t*), we see We assume that the mean field solution

*a*

^{0}

_{i}is stable. In addition, we assume, per the truncation hypothesis, that the steady-state value of

*C*(

_{ij}*t*) does not appreciably alter either

*a*or, therefore, the matrix Γ

_{i}_{ij}. In this case, Γ

_{ij}has all positive eigenvalues and is diagonalizable. Define to be the diagonalization of Γ

_{ij}. We also define the shorthand for the driving terms in equation 3.22. The steady-state solution is given by Notice that each term contributing to the magnitude of

*C*

^{0}

_{ij}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.

*C*are the sums of the eigenvalues of the mean field equation, λ

_{ij}_{i}+ λ

_{j}. In the case that

*a*

^{0}

_{i}is stable, not only will

*C*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 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

_{ij}*f*(

*s*).

*i*has

*N*presynaptic neurons (i.e., neurons for which

_{i}*w*≠ 0). Further, let the average connectivity weight over all inputs be

_{ij}*w*

_{0}. In this case, we can approximate

*w*≈

_{ij}*w*

_{0}/

*N*. The steady-state value of

_{i}*C*(

_{ij}*t*), which is determined by the driving term

*A*(

_{ij}*t*), is seen to be inversely proportional to the number of presynaptic neurons due to the linear dependence on

*w*of

_{ij}*A*. In the most extreme case, the number of presynaptic neurons is the entire network, so

_{ij}*N*=

_{i}*N*, and the correlation function

*C*(

_{ij}*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

*N*= min

_{m}*N*and scale all weights so that they can be written as where

_{i}*w*is the maximum total input to any given neuron. This allows us to treat

_{M}*N*as an effective system size. Larger

_{m}*N*reduces the effects of fluctuations at a given distance to a bifurcation.

_{m}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 *C*^{0}_{ij}, is given by the product of the smallest eigenvalue of Γ_{ij} and the number of presynaptic neurons *N _{i}*. 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 *a _{i}*(

*t*). The higher equations in the hierarchy are affected by this input only through its effect on the mean activity

*a*(

_{i}*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

*a*(

_{i}*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.

*P*are obtained by taking derivatives of the generating function with respect to

*J*. For example, . The cumulants can be obtained by taking the derivatives of

_{i}*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 where is called the action and we use the notation

*u*·

*v*= ∫

*d*(

^{d}xdt 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

*a*

_{0}(

*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 and 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

*t*⩽

*t*′.

*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 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

*n*

_{0}(

*x*). The delta functional is defined by the generalized Fourier transform, and is integrated along the imaginary axis. We can now integrate over the stochastic variable ξ to obtain a noise-generating functional defined by We choose ξ such that which is consistent with a Poisson activation at rate

*F*and a Poisson decay at rate α. We next transform to the new variables: 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

*t*→

*t*′. The properties of the response field provide and so that

*F*becomes Hence, restoring the Ito convention requires the replacement

*n*

^{2}→

*n*

^{2}+

*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.

*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 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

*J*

_{1}=

*J*and . We can “raise” an index μ via multiplication with the matrix 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 For functions of more than one spatial variable, this inner product notation will generalize to the trace of the matrix product: We will also write the action as

*S*[Φ

_{μ}] = ∫

*d*[Φ

^{d}x dt L_{μ}], where

*L*is the integrand of the action in equation 4.3.

*a*

_{μ}= 〈Φ

_{μ}〉. The “effective action” Γ[

*a*

_{μ}] for

*a*

_{μ}is derived by a Legendre transformation, where the conditions are enforced. In analogy with classical mechanics, the extrema of the effective action, give the equations of motion or the activity equations for

*a*

_{μ}(

*x*,

*t*).

*W*[

*J*

^{μ}] using the Legendre transformation, equation 4.17, gives where we have suppressed the

*x*and

*t*arguments. Defining a new variable Ψ

_{μ}= Φ

_{μ}−

*a*

_{μ}and using equation 4.19 gives where we define Γ

^{μ}[

*a*

_{μ}] ≡ δΓ/δ

*a*

_{μ}. We now expand

*S*[Ψ

_{μ}+

*a*

_{μ}] in a functional Taylor series to obtain 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, We introduce a scaling parameter for the action,

*h*(which we will set equal to one), via 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}+

*h*

^{2}Γ

_{2}+ ⋅ ⋅ ⋅ If we set Γ

_{0}=

*S*, we obtain Computing the corrections involves taking expectation values of the operator

*e*

^{−V/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).

*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.

*a*(

*x*,

*t*) and

*C*(

*x*,

*y*,

*t*), we define the composite cumulant generating functional: We now perform a double Legendre transform to obtain the effective composite action, with conditions and 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 and for

*C*

_{μν}, we get together with the conditions 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 (

*C*

_{11}) and the linear response (

*C*

_{1,−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

*C*

_{1,−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

*w*=

_{ij}*w*

_{0}/

*N*for some

*w*

_{0}. The generalized Wilson-Cowan equations become We can simplify this further by taking the initial conditions

*a*(0) =

_{i}*a*

_{0}and

*C*(0) = 0 and assuming homogeneous external input

_{ij}*I*. This corresponds to initial conditions in the network determined by a Poisson distribution with mean

*a*

_{0}.

*C*(0) = −

_{ij}*a*(0)δ

_{i}_{ij}would indicate starting with precisely

*a*(0) active neurons at

_{i}*i*at time

*t*= 0. Then symmetry reduces the equations to where we have defined 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 The steady-state value of

*C*(

*t*) is given by The relative size of the fluctuations in steady state is determined by the product

*N*Γ[

*a*

_{0}]. Large networks or networks distant from a bifurcation experience reduced correlations.

*f*(

*s*) to be 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 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.

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 *a _{i}*(

*t*).

*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 We plot

*a*(

_{i}*t*) and

*C*(

_{ii}*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:

*a*(0) = 2 and

_{i}*C*(0) = 0. The simulations of the Markov process are averaged over 10

_{ij}^{5}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.

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.

*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 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

*C*(

_{ij}*t*) in order to separate the effects of sourcing

*a*(

_{i}*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.

## 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 *a _{i}*(

*t*) is the probability of neuron

*i*being active at time

*t*. The unboundedness of the variable

*a*(

_{i}*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

*h*= 1 for bookkeeping purposes. The generalized effective action is given by with

*J*

^{μ}and

*K*

^{μν}given by The path integral representation is thus which we transform to a new variable Ψ

_{μ}(

*x*,

*t*) = Φ

_{μ}(

*x*,

*t*) −

*a*

_{μ}(

*x*,

*t*), set and expand to obtain where we use We consider the expansion Γ = Γ

_{0}+

*h*Γ

_{1}+

*h*

^{2}Γ

_{2}, where Γ

_{2}contains all terms of order

*h*

^{2}and higher. Setting Γ

_{0}=

*S*gives We now fix

*K*

^{μν}according to which gives where 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 Substituting this into expression A.10 gives us We can extract the normalization of the functional integral using and the identity to obtain 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

*h*

^{2}. Thus, we have

*C*

^{μν}are which we can invert to get In particular, if we ignore loop corrections (consider only first order in

*h*, recalling that Γ

_{2}is

*O*(

*h*

^{2})), we get 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*^{−1}_{m} where *N _{m}* 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.

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 *C*_{1,−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).

## 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.