Even with the rapid increase in computing power due to Moore's law and proposals to simulate the entire human brain notwithstanding Ailamaki et al. (2012), a realistic simulation of a functioning human brain performing non-trivial tasks is remote. While it is plausible that a network the size of the human brain could be simulated in real time Izhikevich and Edelman (2008); Eliasmith et al. (2012) there are no systematic ways to explore the parameter space. Technology to experimentally determine all the parameters in a single brain simultaneously does not exist and any attempt to infer parameters by fitting to data would require exponentially more computing power than a single simulation. We also have no idea how much detail is required. Is it sufficient to simulate a large number of single compartment neurons or do we need multiple-compartments? How much molecular detail is required? Do we even know all the important biochemical and biophysical mechanisms? There are an exponential number of ways a simulation would not work and figuring out which remains computationally intractable. Hence, an alternative means to provide appropriate prior distributions for parameter values and model detail is desirable. Current theoretical explorations of the brain utilize either abstract mean field models or small numbers of more biophysical spiking models. The regime of large but finite numbers of spiking neurons remains largely unexplored. It is not fully known what role spike time correlations play in the brain. It would thus be very useful if mean field models could be augmented with some spike correlation information.

This paper outlines a scheme to derive generalized activity equations for the mean and correlation dynamics of a fully deterministic system of coupled spiking neurons. It synthesizes methods we have developed to solve two different types of problems. The first problem was how to compute finite system size effects in a network of coupled oscillators. We adapted the methods of the kinetic theory of gases and plasmas Ichimaru (1973); Nicholson (1993) to solve this problem. The method exploits the exchange symmetry of the oscillators and characterizes the phases of all the oscillators in terms of a phase density function η(θ, t), where each oscillator is represented as a point mass in this density. We then write down a formal flux conservation equation of this density, called the Klimontovich equation, which completely characterizes the system. However, because the density is not differentiable, the Klimontovich equation only exists in the weak or distributional sense. Previously, e.g., Desai and Zwanzig (1978); Strogatz and Mirollo (1991); Abbott and van Vreeswijk (1993); Treves (1993) the equations were made usable by taking the “mean field limit” of N → ∞ and assuming that the density is differentiable in that limit, resulting in what is called the Vlasov equation. Instead of immediately taking the mean field limit, we regularize the density by averaging over initial conditions and parameters and then expand in the inverse system size N⁻¹ around the mean field limit. This results in a system of coupled moment equations known as the BBGKY moment hierarchy. In Hildebrand et al. (2007), we solved the moment equations for the Kuramoto model perturbatively to compute the pair correlation function between oscillators. However, the procedure was somewhat ad-hoc and complicated. We then subsequently showed in Buice and Chow (2007) that the BBGKY moment hierarchy could be recast in terms of a density functional of the phase density. This density functional could be written down explicitly as an integral over all possible phase histories, i.e., a Feynman-Kac path integral. The advantage of using this density functional formalism is that the moments to arbitrary order in 1/N could be computed as a steepest-descent expansion of the path integral, which can be expressed in terms of Feynman diagrams. This made the calculation more systematic and mechanical. We later applied the same formalism to synaptically coupled spiking models Buice and Chow (2013b).

Concurrently with this line of research, we also explored the question of how to generalize population activity equations, such as the Wilson-Cowan equations, to include the effects of correlations. The motivation for this question is that the Wilson-Cowan equations are mean field equations and do not capture the effects of spike-time correlations. For example, the gain in the Wilson-Cowan equations is fixed, (which is a valid approximation when the neurons fire asynchronously), but correlations in the firing times can change the gain Salinas and Sejnowski (2000). Thus, it would be useful to develop a systematic procedure to augment population activity equations to include spike correlation effects. The approach we took was to posit plausible microscopic stochastic dynamics, dubbed the spike model, that reduced to the Wilson-Cowan equations in the mean field limit and compute the self-consistent moment equations from that microscopic theory. Buice and Cowan (2009) showed that the solution of the master equation of the spike model could be expressed formally in terms of a path integral over all possible spiking histories. The random variable in the path integral is a spike count whereas in the path integral for the deterministic phase model we described above, the random variable is a phase density. To generate a system of moment equations for the microscopic stochastic system, we transformed the random spike count variable in the path integral into moment variables Buice et al. (2010). This is accomplished using the effective action approach of field theory, where the exponent of the cumulant generating functional, called the action, which is a function of the random variable is Legendre transformed into an effective action of the cumulants. The desired generalized Wilson-Cowan activity equations are then the equations of motion of the effective action. This is analogous to the transformation from Lagrangian variables of position and velocity to Hamiltonian variables of position and momentum. Here, we show how to apply the effective action approach to a deterministic system of synaptically coupled spiking neurons to derive a set of moment equations.

Consider a network of single compartment conductance-based neurons

The equations are remarkably stiff with time scales spanning orders of magnitude from milliseconds for ion channels, to seconds for adaptation, and from hours to years for changes in synaptic weights and connections. Parameter values must be assigned for 10¹¹ neurons with 10⁴ connections each. Here, we present a formalism to derive a set of reduced activity equations directly from a network of deterministic spiking neurons that capture the spike rate and spike correlation dynamics. The formalism first constructs a density functional for the firing dynamics of all the neurons in a network. It then systematically marginalizes the unwanted degrees of freedom to isolate a set of self-consistent equations for the desired quantities. For heuristic reasons, we derive an example set of generalized activity equations for the first and second cumulants of the firing dynamics of a simple spiking model but the method can be applied to any spiking model.

A convenient form to express spiking dynamics is with a phase oscillator. Consider the quadratic integrate-and-fire neuron

where I is an external current and u(t) are the synaptic currents with some weight α_i. The spike is said to occur when V goes to infinity whereupon it is reset to minus infinity. The quadratic non-linearity ensures that this transit will occur in a finite amount of time. The substitution V = tan(θ/2) yields the theta model Ermentrout and Kopell (1986):

which is the normal form of a Type I neuron near the bifurcation to firing Ermentrout (1996). The phase neuron is an adequate approximation to spiking dynamics provided the inputs are not overly strong as to disturb the limit cycle. The phase neuron also includes realistic dynamics such as not firing when the input is below threshold. Coupled phase models arise naturally in weakly coupled neural networks Ermentrout and Kopell (1991); Hoppensteadt and Izhikevich (1997); Golomb and Hansel (2000). They include the Kuramoto model Kuramoto (1984), which we have previously analyzed Buice and Chow (2007); Hildebrand et al. (2007).

Here, we consider the phase dynamics of a set of N coupled phase neurons obeying

where each neuron has a phase θ_i that is indexed by i, u is a global synaptic drive, F(θ, γ, u) is the phase and synaptic drive dependent frequency, γ_i represents all the parameters for neuron i drawn from a distribution with density g(γ), ν is the population firing rate of the network,t^l_j is the lth firing time of neuron j and a neuron fires when its phase crosses π. In the present paper, we consider all-to-all or global coupling through a synaptic drive variable u(t). However, our basic approach is not restricted to global coupling.

We can encapsulate the phase information of all the neurons into a neuron density function Buice and Chow (2007, 2011, 2013a,b); Hildebrand et al. (2007).

where δ(·) is the Dirac delta functional, and θ_i(t) is a solution to system (3–5). The neuron density gives a count of the number of neurons with phase θ and synaptic strength γ at time t. Using the fact that the Dirac delta functional in (5) can be expressed as ∑lδ(t−tjl)=θj.δ(π−θj), the population firing rate can be rewritten as

The neuron density formally obeys the conservation equation

with initial condition η(θ, γ, t₀) = η₀(θ, γ) and u(t₀) = u₀. Equation (8) is known as the Klimontovich equation Ichimaru (1973); Liboff (2003). The Klimontovich equation, the equation for the synaptic drive (4), and the firing rate expressed in terms of the neuron density (7), fully define the system. The system is still fully deterministic but is now in a form where various sets of reduced descriptions can be derived. Here, we will produce an example of a set of reduced equations or generalized activity equations that capture some aspects of the spiking dynamics. The path we take toward the end will require the introduction of some formal machinery that may obscure the intuition around the approximations. However, we feel that it is useful because it provides a systematic and controlled way of generating averaged quantities that can be easily generalized.

For finite N, (8) is only valid in the weak or distributional sense since η is not differentiable. In the N → ∞ limit, it has been argued that η will approach a smooth density ρ that evolves according to the Vlasov equation that has the same form as (8) but with η replaced by ρ Ichimaru (1973); Desai and Zwanzig (1978); Strogatz and Mirollo (1991); Nicholson (1993); Hildebrand et al. (2007). This has been proved rigorously in the case where noise is added using the theory of coupled diffusions McKean Jr (1966); Faugeras et al. (2009); Baladron et al. (2012); Touboul (2012). This N → ∞ limit is called mean field theory. In mean field theory, the original microscopic many body neuronal network is represented by a smooth macroscopic density function. In other words, the ensemble of networks prepared with different microscopic initial conditions is sharply peaked at the mean field solution. For large but finite N, there will be deviations away from mean field Buice and Chow (2007); Hildebrand et al. (2007); Buice and Chow (2013a,b). These deviations can be characterized in terms of a distribution over an ensemble of coupled networks that are all prepared with different initial conditions and parameter values. Here, we show how a perturbation theory in N⁻¹ can be developed to expand around the mean field solution. This requires the construction of the probability density functional over the ensemble of spiking neural networks. We adapt the tools of statistical field theory to perform such a construction.

The main goal of this paper was to show how to systematically derive generalized activity equations for the ensemble averaged moments of a deterministically coupled network of spiking neurons. Our method utilizes a path integral formalism that makes the process algorithmic. The resulting equations could be derived using more conventional perturbative methods although possibly with more calculational difficulty as we found before Buice et al. (2010). For example, for the case of the stochastic spike model, Buice et al. (2010) presumed that the Wilson-Cowan activity variable was the rate of a Poisson process and derived a system of generalized activity equations that corresponded to deviations around Poisson firing. Bressloff (2010), on the other hand, assumed that the Wilson-Cowan activity variable was a mean density and used a system-size expansion to derive an alternative set of generalized activity equations for the spike model. The classical derivations of these two interpretations look quite different and the differences and similarities between them are not readily apparent. However, the connections between the two types of expansions are very transparent using the path integral formalism.

Here, we derived equations for the rate and covariances (first and second cumulants) of a deterministic synaptically coupled spiking network as a system size expansion to first order. However, our method is not restricted to these choices. What is particularly advantageous about the path integral formalism is that it is straightforward to generalize to include higher order cumulants, extend to higher orders in the inverse system size, or to expand in other small parameters such as the inverse of a slow time scale. The action fully specifies the system and all questions regarding the system can be addressed with it.

To give a concrete illustration of the method, we derived the self-consistent generalized activity equations for the rates and covariances to order N⁻¹ for a simple phase model. The resulting equations consist of ordinary and partial differential equations. This is to be expected since the original system was fully deterministic and memory cannot be lost. Even mean field theory requires the solution of an advective partial differential equation. The properties of these and similar equations remain to be explored computationally and analytically. The system is possibly simpler near the asynchronous state, which is marginally stable in mean field theory like the Kuramoto model Strogatz and Mirollo (1991) and like the Kuramoto model, we conjecture that the finite size effects will stabilize the asynchronous state Buice and Chow (2007); Hildebrand et al. (2007). The addition of noise will also stabilize the asynchronous state. Near asynchrony could be exploited to generate simplified versions of the asynchronous state.

We considered a globally connected network, which allowed us to assume that networks for different parameter values and initial conditions converge toward a “typical” system in the large N limit. However, this property may not hold for more realistic networks. While the formalism describing the ensemble average will hold regardless of this assumption, the utility of the equations as descriptions of a particular network behavior may suffer. For example, heterogeneity in the connectivity (as opposed to the global connectivity we consider here) may threaten this assumption. This is the case with so called “chaotic random networks” Sompolinsky et al. (1988) in which there is a spin-glass transition owing to the variance of the connectivity crossing a critical threshold. This results in the loss of a “typical” system in the large N limit requiring an effective stochastic equation which incorporates the noise induced by the network heterogeneity. Whether the expansion we present here is useful without further consideration depends upon whether the network heterogeneity induces this sort of effect. This is an area for future work. A simpler issue arises when there are a small discrete number of “typical” systems (such as with bistable solutions to the continuity equation). In this case, there are noise induced transitions between states. While the formalism has a means of computing this transition Elgart and Kamenev (2004), we do not consider this case here.

An alternative means to incorporate heterogeneous connections is to consider a network of coupled systems. In such a network, a set of generalized activity equations, such as those derived here or simplified versions, would be derived for each local system, together with equations governing the covariances between the local systems. Correlation based learning dynamics could then be imposed on the connections between the local systems. Such a network could serve as a generalization of current rate based neural networks to include the effects of spike correlations with applications to both neuroscience and machine learning.