A stochastic-field description of finite-size spiking neural networks

Neural network dynamics are governed by the interaction of spiking neurons. Stochastic aspects of single-neuron dynamics propagate up to the network level and shape the dynamical and informational properties of the population. Mean-field models of population activity disregard the finite-size stochastic fluctuations of network dynamics and thus offer a deterministic description of the system. Here, we derive a stochastic partial differential equation (SPDE) describing the temporal evolution of the finite-size refractory density, which represents the proportion of neurons in a given refractory state at any given time. The population activity—the density of active neurons per unit time—is easily extracted from this refractory density. The SPDE includes finite-size effects through a two-dimensional Gaussian white noise that acts both in time and along the refractory dimension. For an infinite number of neurons the standard mean-field theory is recovered. A discretization of the SPDE along its characteristic curves allows direct simulations of the activity of large but finite spiking networks; this constitutes the main advantage of our approach. Linearizing the SPDE with respect to the deterministic asynchronous state allows the theoretical investigation of finite-size activity fluctuations. In particular, analytical expressions for the power spectrum and autocorrelation of activity fluctuations are obtained. Moreover, our approach can be adapted to incorporate multiple interacting populations and quasi-renewal single-neuron dynamics.


Author summary
In the brain, information about stimuli is encoded in the timing of action potentials produced by neurons. An understanding of this neural code is facilitated by the use of a wellestablished method called mean-field theory. Over the last two decades or so, mean-field theory has brought an important added value to the study of emergent properties of neural circuits. Nonetheless, in the mean-field framework, the thermodynamic limit has to be taken, that is, to postulate the number of neurons to be infinite. Doing so, small fluctuations are neglected, and the randomness so present at the cellular level disappears from the description of the circuit dynamics. The origin and functional implications of variability at the network scale are ongoing questions of interest in neuroscience. It is therefore

Introduction
Neurons communicate by sending and receiving pulses called spikes which occur in a rather stochastic fashion. A stimulus is thus translated by neurons into spike trains with a certain randomness [1]. At the microscopic scale, this variability is mostly attributed to the probabilistic nature of the opening and closing of ion channels underlying the emission of an action potential. At a mesoscopic scale, variability typically stems from the seemingly random barrage of synaptic inputs. This variability is fundamentally noise [2,3]. Many papers have been devoted to establish its origin [4,5], and the mathematical formalization to suitably describe this effect has been an intense subject of research over the past decades. Nowadays, models are capable of reproducing both the statistics of the spiking activity and the subthreshold dynamics of different cell types.
In neuroscience, it is believed that information-about external stimuli, internal states, motor outputs, etc.-is encoded in the timing of spikes produced by populations of neurons [6]. An understanding of this high-dimensional response, from an information-theoretic or dynamical point of view, is facilitated by various dimensionality-reduction methods [7]. A trivial one is to consider the population activity, i.e. the proportion of neurons firing in small time windows, which is assumed to be a meaningful coarse-grained dynamical variable [6,[8][9][10]. The importance of population activity is manifest in its extensive use in macroscopic models of neural activity, and by the constant effort put forth to derive its dynamics from singleneuron and network properties.
Most attempts to produce analytically tractable population rate models have made use (directly or indirectly) of mean-field theory [11][12][13][14]. The population activity obtained by solving mean-field models is deterministic, since this theory neglects finite-size fluctuations. Theoretically, this causes no problem. Real neural circuits, however, necessarily have a finite size. For a system made up of N independent units, the relative magnitude of fluctuations should scale as 1= ffiffiffiffi N p . The thermodynamic limit (N goes to infinity) neglects those small fluctuations, and the randomness so present at the cellular level disappears in the description of the circuit. However, these finite-size effects should be taken into account because they can drastically affect both the synchronization [15,16] and the stability [17,18] of neural systems.
Various analytical methods from statistical physics have been used to describe such activity fluctuations. A first type of approach consists in adapting the Fokker-Planck formalism to incorporate the finite-size rate fluctuations as a source term [18,19]. The spiking processes of the neurons are then assumed to be Poisson processes. One can also apply the so-called linear response theory (LRT) [20] and compute spectral quantities characterizing the neuronal spiking processes. This theory relies on the ansatz that the spiking response of a neuron embedded in a network can be written as the sum of the neuron's unperturbed spiking activity-i.e., when the neuron is isolated from the fluctuations of the network activity-and its first-order response to the fluctuations of synaptic currents. Finite-size effects do not need a special treatment in that theory, but it can only manage wide-sense stationary external inputs that maintain a comfortable distance from bifurcation points. LRT has in fact been successfully used to study finite-size fluctuations in the context of decorrelative effects of inhibitory feedback [21], the interplay between correlation structures and network topology [22], and the effect of correlations on coding [23].
A more systematic approach to study finite-size effects is to construct a master equation describing the time evolution of the probability distribution of the states of the network. In that case the formalism borrows heavily from reaction kinetics in chemistry [24,25]. The network navigates between states via activation and decay functions, which embody the stochastic transitions between quiescent and active spiking states. Using the master equation, one can then construct a system involving moments of the network's activity. A troublesome feature of this approach is that lower order moments typically depend on higher order moments, thus constituting a hierarchy of moments that must be truncated in one way or the other. The chosen truncation scheme depends on the assumptions about the underlying neural activity. One possibility is to assume that the network is in an asynchronous state-defined by a constant infinite-size mean field-so that spike-spike correlations are of order 1/N, and a direct expansion of the moments to that order becomes possible [26][27][28]. In the same spirit, one can also perform a system-size expansion-to order 1= ffiffiffiffi N p -of the master equation, and then construct the hierarchy of moments that can be subsequently truncated [29]. Another truncation based on normal-ordered cumulants assumes near-Poisson spiking [30].
The truncation of the hierarchy of moment equations breaks down near criticality, i.e., near a bifurcation. One way to circumvent this problem is to use path integral methods [31]borrowed from statistical field theory-which themselves are amenable to the use of renormalization group methods that can extract critical exponents near the bifurcation. These field theoretic methods have also been applied to describe the statistics of deterministic phase neurons [32]. All the approaches discussed above assume a network that is either far below criticality or near criticality.
As we see from this overview, analytical treatments are emerging and hold the expectancy of understanding the nontrivial effects of variability on neural circuits. To move toward the formulation of models that keep track of the intrinsic randomness, our challenge is to correct the usual mean-field equations to account for the inescapable stochastic nature of spike initiation. The aim of the present paper is then to take up the problem of finite-size fluctuations and to show that, actually, one can formulate it in the framework of stochastic partial differential equations (SPDE).
Contrary to other treatments of finite-population dynamics, our main objective is to describe the finite-size activity itself, and not its moments. To this end, we derive a SPDE (see Eq 7 below) that gives the activity of a finite population of spike-response model neurons with escape noise [33] for a fully connected inhibitory network. The equation describes the dynamics of the finite-size refractory density [34], i.e. the density of neurons in a given refractory state at a given time. The boundary condition (Eq 8) or the conservation law (Eq 9) for the refractory density is used to extract the activity. Finite-size fluctuations appear in the SPDE through a two-dimensional Gaussian white noise-in time and along the refractory dimension-whose prefactor vanishes in the infinite-size (thermodynamic) limit. Importantly, the Gaussian white noise acting in the network's description naturally emerges from the intrinsic randomness of spike initiation present at the cellular level.
The SPDE can be solved numerically via a discretization along its characteristic curves (see Eq 55 in the Methods section), and thus provides a direct mean to simulate finite-size networks, both below and above the bifurcation towards the oscillatory state. Importantly, the simulation time does not depend on the size of the population. Bifurcation analysis of the associated mean-field counterpart enables us to reveal how delayed inhibitory feedback permits the emergence of macroscopic rhythms. More insight into finite-size effects is obtained by applying a linear noise approximation, followed by a study of the spectral and autocorrelation properties of fluctuations in the asynchronous activity regime (see Eqs 18 and 20).
Perhaps the approaches closest in spirit to the one adopted in the following are those of Meyer and van Vreeswijk [27], and Deger, Schwalger and coworkers [35,36]. Meyer and van Vreeswijk treated finite-size fluctuations in homogeneous, fully-connected networks of renewal spike-response neurons using a refractory density approach, as in the current paper. However, their main goal was to derive equations for the temporal cross correlations in the stable stationary state. Thus, even though their framework is akin to ours, the objectives differ: instead of focusing on the moments, here we analytically build a stochastic dynamical equation that captures the temporal evolution of the activity itself, including finite-size fluctuations.
The paper by Deger, Schwalger et al. [35] dealt with the finite-size activity of randomly connected spike-response model neurons with adaptation. With their approximations, they ended up with quasi-renewal neurons [37] connected in an effective all-to-all fashion. In that context, they established an integral equation for the activity, and obtained the temporal autocorrelation function of the activity. However, this formulation makes it impossible to actually simulate the network's activity; this is because their correlation function includes the probability density that a spike in the past causes a spike in the future, which cannot be computed in general because the future activity-which is yet unknown-is required.
In a recent work, however, Schwalger et al. [36] solved that problem by proposing an error minimizing method that permits rapid and accurate simulations of the firing statistics of the network, for a single or multiple populations. They derived stochastic differential equations that only involve the finite-size activity; contrary to our approach, information about the refractory distribution is thus purposefully disregarded. A byproduct of this approximation is that the neuron number is not conserved anymore, unlike our approach. Moreover, using a Gaussian approximation they obtained a stochastic equation-Eq 41 in their appendix-that is similar in spirit to ours, although it differs somewhat in its details. Their stochastic equation also involves the temporal evolution of the refractory distribution, but is not thoroughly analyzed.
This paper is organized as follows. First, we present the network and neuron model that will be used throughout. Then, we obtain the SPDE and perform a linear noise approximation, which is then used to study both the stability of the deterministic (mean-field) part of the activity, and the statistical properties of the finite-size fluctuations. These include spectral and autocorrelation properties.

The neural network
We consider a fully connected (all-to-all) homogenous network of N neurons with inhibitory synapses. Neuronal dynamics are described using the spike-response model with escape noise [33]. For this model, neuronal spiking occurs according to an instantaneous firing intensity or hazard function that depends on the difference between the membrane potential and a fixed threshold. The membrane potential of neuron i at time t is given by ðk Ã y j ÞðtÞ þ ðk Ã I ext ÞðtÞ À ðV Ã y i ÞðtÞ þ h i ðtÞ; where h i (t) is the change in potential caused by inputs (input potential), À V is a refractory kernel representing the reset of the membrane potential following a spike, and κ is a filter kernel encompassing synaptic and membrane filtering of synaptic and external inputs. The spike train of any neuron is given by where the sum is over all of its spike times. Hence, the membrane potential of neuron i at time t is given by the convolution of its own spike train with the refractory kernel [ðV Ã y i ÞðtÞ], to which are added a filtered external input [(κ Ã I ext )(t)] and the total inhibitory synaptic input coming from all neurons within the network, with uniform synaptic weight −J s /N. We shall restrict ourselves to non-adaptive single-neuron dynamics, meaning that only the most recent spike of a neuron affects its potential. Thus, where is the refractory state or age of neuron i, andt i denotes its most recent spike. We are interested in the dynamics and statistics of the population activity The homogeneity of the network implies that the index i can be dropped in Eq 1: all neurons with the same refractory state have the same subsequent dynamics. Using the definition of the population activity, Eq 1 becomes uðt; rÞ ¼ ðk Ã ½I ext À J s AÞðtÞ À VðrÞ: The hazard function ρ is the probability per unit time of emitting a spike, conditional on the past history of the activity, H t fAðt 0 Þ : 0 < t 0 < tg; and of the external signal, fI ext ðt 0 Þ : 0 < t 0 < tg: For concreteness we will use an exponential hazard function, where λ 0 and δu = 1 mV are constants. This choice has no impact on the theory presented herein, other than simplifying some computations. Also, the refractory kernel is taken to be where τ is the recovery time scale. The synaptic filter is Finite-size spiking neural networks with H(t) the Heaviside step function, τ s the synaptic decay and Δ the conduction delay. The specific expression characterizing the escape rate is justified and largely used because it combines a relative mathematical simplicity with the capacity to realistically reproduce observed neuronal behaviors. Note that J s and I ext have units mVÁms and mV, respectively, because κ has units ms −1 . The synaptic kernel is thus defined because in that case the average input potential is independent of the time scale τ s . The hazard function and the synaptic kernel are depicted in Fig 1. We also provide examples of network dynamics in Fig 2. The dynamics of every neuron follow the non-adaptive version of Eq 1 (i.e., with Eq 2), together with the escape rate firing mechanism and the hazard function of Eq 5. As J s increases, the network develops an oscillatory instability (see below) and oscillations appear.
The stochastic-field equation The continuous deterministic mean-field approach to modeling neural networks fails to capture many important details. The missing detail manifests itself as small unpredictable finite size fluctuations present at the network level. Our main challenge is then to define an equation Finite-size spiking neural networks that fully entails the stochastic aspects of the network. To do so, we consider the finite-size refractory density q(t, r), such that Nq(t, r)dr gives the number of neurons with age in [r − dr, r) at time t. A time dt later, all neurons that did not fire will have aged by an amount dr = dt, whereas the number of neurons that did fire is given by a Poisson random variable P with rate Nρ(t, r)q(t, r)dtdr [25,35]. This idea encompasses the presence of fluctuations that are proportional to the mean number of firing events and therefore retains the full random character of the spiking probability. Using a Gaussian approximation of this Poisson distribution, i.e. making the approximation where N ð0; 1Þ denotes the standard normal distribution and * means "is distributed like", we show (see Methods) that q(t, r) obeys the following stochastic partial differential equation ( where η is a Gaussian (sheet) white noise with hZðt; rÞi ¼ 0; hZðt; rÞZðt 0 ; r 0 Þi ¼ dðt À t 0 Þdðr À r 0 Þ; and rðt; rÞ ¼ l 0 exp ½uðt; rÞ: The brackets denote an ensemble average over realizations of the stochastic process. The boundary condition is naturally given by the reset mechanism. Indeed, once a neuron triggers a spike, its age is reset to zero. Therefore, the boundary condition is where A(t) is the finite-size activity of the network (see for instance Fig 2). This activity A(t) is also given by the total rate at which neurons of all ages escape their trajectories in the (t, r)plane: AðtÞ ¼ Finite-size spiking neural networks The integral above is to be understood in the Itô sense. Finally, the refractory density must obey a conservation law, since each neuron possesses a given age at any given time. Note that both q(t, r) and A(t) are stochastic processes in this formulation. From the derivation above we observe that for a network with N cells, the relative magnitude of fluctuations scale as 1= ffiffiffiffi N p . The stochastic parts in Eq 7 as well as in Eq 9 disappear in the thermodynamic limit when the network is taken infinitely large. Doing so, the equation finally reduces to the classical refractory equation [33]. However, the thermodynamic limit does not allow for any characterization of fluctuations around the mean activity, because it is an entirely deterministic approach as is illustrated in Figs 3 and 4. Fig 3A shows the time evolution of the external stimulus, whereas panel B gives the spiking activity obtained from a simulation of the full network and its corresponding activity. The Finite-size spiking neural networks activity given by the standard mean field in the thermodynamic limit (N taken infinitely large, see Eq 13 below) is shown in Fig 3C. Although in this case the mean-field approximation captures the essential "shape" of the activity of the full network shown in the panel B, it completely ignores the substantial finite-size fluctuations. Finally, Fig 3D shows the activity as generated via a simulation of the stochastic-field equation (Eq 7). The stochastic-field equation effectively captures both the shape and variability of the full neural activity, which could be described as a noisy version of the mean-field activity. Similar observations can be made regarding the refractory density q(t, r) as can be seen in Fig 4. The numerical integration of Eq 7 is discussed in the Methods section.
Note that multiple populations can be modeled straightforwardly using the above formalism. To each population n would be assigned a refractory density q n (t, r) obeying the SPDE. The respective membrane potentials would be given by u n ðt; rÞ ¼ ðk n Ã ½I ext À X k J nk A k ÞðtÞ À V n ðrÞ with J nk the total synaptic strength connecting population k to population n. Different hazard functions ρ n can be chosen as well.
An analytical solution of the SPDE is exceedingly difficult, if not impossible. However, informations about the statistics of activity fluctuations can be extracted via a system-size expansion, as discussed in the next section.

Linear noise approximation
To investigate the effects of fluctuations for large but finite network size N, we perform a linear noise approximation (LNA) when I ext is constant. In our situation, the LNA states that the density function as well as the neural activity can be approximated by the sum of a deterministic and a stochastic process. The fluctuating part is scaled by a 1= ffiffiffiffi N p factor that is justified by the Van Kampen system-size expansion [24]. This system-size expansion is usually pursued only up to first order, hence the "linear" qualifier. We write with q 0 (t, r) and q ξ (t, r) the deterministic and stochastic parts, respectively. Similarly, the population activity reads again with A 0 (t) the deterministic part and A ξ (t) the stochastic part. After algebraic manipulations-including a linearization that only keeps the first-order term in 1= ffiffiffiffi N p -we find that the deterministic part obeys the usual mean-field description with boundary condition where A 0 (t) is the deterministic part of the activity. The hazard ρ 0 (t, r) is r 0 ðt; rÞ r½u 0 ðt; rÞ; with u 0 ðt; rÞ ðk Ã ½I ext À J s A 0 ÞðtÞ À VðrÞ: Note that the mean field equation is strikingly similar to the standard age-structured system known as the McKendrick von-Foerster model in mathematical biology [38,39]. The stochastic component solves a SPDE similar to Eq 7, namely (see Methods) and the derivative dr=duj u 0 is evaluated at u 0 . This equation is now linear in the function q ξcontrary to Eq 7-and therefore an analytical solution is possible (this is done in Methods in the context of a computation of the power spectrum of activity fluctuations).
When the deterministic component exhibits multiple stable states, the LNA fails to capture transitions between these states. While the LNA is inadequate for this type of situation, it is pertinent when the network fluctuates around a unique stable equilibrium, as illustrated in Fig 5. The stochastic activity is given to first order by Eq 12. After a short transient (*10 ms), the deterministic part of the activity, A 0 (t), asymptotically reaches its steady state value A 1 (see black curve in Fig 5C and 5E). Since the neuron number is finite, fluctuations around the deterministic activity are observed, both in the transient and steady states (red curve Fig 5C  and 5E). As illustrated in panels B and D of this figure, the activity of the relaxed network is asynchronous.
The average activity in the asynchronous state, A 1 , can be computed as the time-independent solution of the infinite-size refractory density equation, Eq 13. This means that the partial derivative with respect to time is zero. After algebraic manipulations, we find that A 1 is given by (see Methods) where r 1 ðsÞ r½u 1 ðsÞ ¼ r½h 1 À VðsÞ; Finite-size spiking neural networks and h 1 ¼ I ext À J s A 1 : The equation above must be solved self-consistently for A 1 . From Fig 5A, the agreement between analytical and numerical results is excellent.

Deterministic oscillatory instability
The asynchronous state can lose stability through an oscillatory instability. Finite-size fluctuations will potentially influence the presence of oscillations. To address this issue, we performed a linearization of the deterministic field equation (Eq 13) around the steady state (see Methods). This allowed us to write the characteristic equation, whose solutions give the eigenvalues λ of the deterministic system. The characteristic equation is denoted wherekðlÞ is the Laplace transform of κ and P 1 (r) is the steady-state interspike interval Finite-size spiking neural networks probability density of the asynchronous infinite-size network, i.e. the probability density of obtaining an interval see Eq 46. Also, is the steady state of the infinite-size refractory density. The eigenvalues are thus given by the roots of CðlÞ. The time-independent solution A 1 will be stable if all eigenvalues have negative real parts. The bifurcation line-which separates the oscillatory regime from the asynchronous one-can be drawn numerically, as depicted in Fig 6A. The red line constitutes the boundary between the stability and the instability regions for the chosen I ext . The shaded area defines the region in parameter space where self-sustained oscillations are going to emerge. Well below the transition line, in the asynchronous regime, damped oscillations may occur; however, the network activity will eventually settle into finite-size fluctuations around a constant mean activity. Above the transition, network oscillations are clearly noticeable ( Fig 6C).
The existence of oscillations in fully connected, delayed inhibitory networks is well known [40], and underlies the ING (Interneuron Network Gamma) mechanism for generating gamma rhythms (see for instance [41]). The above analysis serves the purpose of delineating the oscillatory and asynchronous states in the phase diagram of Fig 6 with the help of the characteristic function. Both the phase diagram and the characteristic function will come handy in the next section where we will study the fluctuations of the asynchronous state. Note, however, that we shall not study the exact nature of the bifurcation; for the network studied in [40], the Hopf bifurcation can be either subcritical or supercritical, depending on the inhibitory coupling strength.

Fluctuations in the asynchronous regime
Asynchronous firing has been repeatedly associated with the background activity in cortical networks [42]. Finite-size effects generate fluctuations whose basic statistics are closely related Finite-size spiking neural networks to the response of the associated infinite-size dynamics. To characterize the statistical content of fluctuations we compute the power spectrum of the activity in the asynchronous regime. The power spectrum describes how the variance of A(t) is distributed over the frequency domain and accordingly helps to identify the dominant frequencies, if any.
By definition, the power spectrum is given by is the Fourier transform of the neural activity restricted to a time interval and the brackets hÁi denote an average over the noisy realizations of the network dynamics. To derive an analytical expression for this quantity, we assume that the deterministic part of the activity has reached its equilibrium state (below the threshold of instability). This means that Eq 12 becomes neglecting terms of order 1/N and above. Hence, with P x ðoÞ=N the power spectrum of the fluctuations. In Methods, P x ðoÞ is shown to be where < means to take the real part of the expression in curly brackets and P 1 (r) is the interspike interval distribution in the asynchronous state. The stability of the deterministic system-embodied in its characteristic equation CðlÞ ¼ 0-appears explicitly in the above expression. Of course, in the asynchronous regime there are no pure imaginary solutions to this equation. However, the minima of CðioÞ dictate the position of dominant frequencies, as illustrated in Fig 7. This figure also compares the power spectrum obtained from Eq 18 to the power spectrum obtained from both an average over different realizations of the fully stochastic neural network and the SPDE; it shows an excellent agreement. We carried out the comparison for more simulations than are shown here, and our results hold over a wide range of parameters. We note that there is a discrepancy between theory and numerics mainly at lower frequencies as the bifurcation is approached from below (Fig 7B). The theoretical power spectrum of Eq 18 does not capture the emergence of harmonics (for instance, the small peak around f = 400 Hz in Fig 7B). The LNA is known to fail as soon as parameters are chosen to be close to a bifurcation point. The spectra computed from the SPDE agree very well with the spectra from the whole network for all the regimes shown. Thus the SPDE can be used to provide a Finite-size spiking neural networks good picture of a spiking network's spectral properties without having to resort to the much longer full network simulations, at least when the network size is not too small.
Another important remark concerns the discrepancy between the SPDE and the full network for very small system sizes. To illustrate this effect, we show in Fig 8 comparisons    Finite-size spiking neural networks between the power spectrum of the full network and the SPDE for different sizes. The simulations presented in Fig 8 correspond to parameters under the bifurcation line (top row), close to the bifurcation line (middle row), and above the bifurcation (bottom row). The chosen parameter values within these regimes correspond to the three stars of Fig 6A. Panels A, B and C are illustrations for different number of cells. We clearly see that the SPDE gives very good results for large N, for all regimes. However, discrepancies appear when the number of cells becomes too small (of order *10; panels C). This is no surprise since our theory is expected to fail for very small networks. Indeed, the Gaussian approximation to the Poisson variable (see Methods) that is used to arrive at the SPDE (Eq 7) implicitly assumes the number of cells to be large enough. Hence the observed discrepancy when the number of cells is too small.
Another way to obtain the power spectrum would be to first compute the autocorrelation function of A ξ (t) and then use the Wiener-Khintchin theorem. However, solving Eq 14 provides a self-consistent expression (see Methods), namely Here, Q represents the response to finite-size perturbations of the activity in the past, whereas S gives the effect of the Gaussian noise η(t, r) on the stochastic part of the activity (see Eq 54 and following in Methods for expressions for Q and S). Eq 19 can be solved by taking the Fourier transform, thus forcing us to first compute the power spectrum to get the autocorrelation function. Interestingly, however, the autocorrelation function of S(t) can be obtained directly in terms of the interspike interval distribution in the asynchronous state, giving (see Methods) which corresponds to that obtained in [35].

Discussion
The origin and functional implications of variability at the network scale are ongoing questions of interest in neuroscience. There have been a number of earlier efforts to go beyond the mean-field approximation to address these questions. In the past few years the idea of studying spiking neural circuit within the framework of statistical physics to investigate finite-size effects has become a concrete research project (see e.g. [43]). In this study, an alternative method was proposed to adequately describe observable noise at the network level. We derived a stochastic partial differential equation (SPDE) describing the dynamics of the refractory density q(t, r) (Eq 7). This quantity gives the density of neurons in refractory state r at a certain time t, the refractory state being the time interval from the last spiking event. The population activity A(t) can be extracted, for instance, from the boundary condition q(t, 0). In the limit where the neuron number N goes to infinity, the standard PDE for the infinite-size refractory density is recovered.
An important point about our derivation is that we did not assume that the single-neuron spiking processes were described by a (inhomogeneous) Poisson process. In particular, as per renewal theory, the firing intensity (or hazard function) depended on the last spike, contrary to Poisson processes. Poisson random variables appeared in the course of the derivation because the number of neurons firing in a small time interval Δt is a random variable following a binomial distribution, which becomes a Poisson distribution in the limit Δt goes to zero. From there, the only assumption to arrive at Eq 7 was to approximate this Poisson distribution by a Gaussian distribution-the Gaussian approximation. This is in contrast with the work of Brunel and Hakim [18,44], which starts with deterministic single-neuron membrane equations, but Finite-size spiking neural networks approximates the spiking process of neurons embedded in the network by a Poisson process. Nonetheless, the two approaches lead to a similar description, where the number of spikes produced during a small amount of time can be approximated by a Gaussian random variable.
Our SPDE cannot be solved analytically since the stochastic term involves a square-root nonlinearity. Nonetheless, a discretization along the characteristic curves of the SPDE offers a numerical scheme that gives a very satisfying approximation to its behavior. The advantage over numerical simulations of the full network is that it overcomes some of the restrictions imposed by computation time in multiple neuron simulations. Indeed, the numerical simulation of the SPDE is independent of the number of neurons. Therefore the stochastic-field density approach appears to be a viable method to make simulations of large neural networks all the while keeping track of their intrinsic variability.
On the other hand, a system-size expansion restricted to first order in 1= ffiffiffiffi N p -the linear noise approximation-gives rise to two coupled linear SPDEs which are amenable to analytical investigations. In particular, we studied finite-size fluctuations in the asynchronous regime, for which the average network activity is constant and neurons fire asynchronously. This enabled an analytical expression for the power spectrum of these fluctuations to be obtained (Eq 18). Its structure is determined both by the characteristic function of the deterministic system (thermodynamic limit) and by the spiking properties of the network (e.g., the interspike interval density and the hazard function). The characteristic function appears in the stability analysis of the deterministic system: zeros of this function are eigenmodes of the deterministic dynamics. Similarly, via an integration along the characteristic curves, some salient aspects of the autocorrelation properties of the finite-size fluctuations in the asynchronous activity regime were computed (Eq 20), allowing comparison with the literature. We therefore note that our SPDE permits the analytical computations of both the spectral Eq (18) and the correlation Eq (20) properties of the fluctuations.
We restricted our study to fully connected networks of non-adaptive (renewal) neurons. However, it should be noted that it is possible to take care of non-renewal aspects of neural dynamics using the approximation provided by the quasi-renewal theory [37]. That theory suggests that one can replace the full adaptation potential (Eq 1) by the sum of a contribution from the last spike, À VðrÞ, and an average potential caused by other spikes over the whole past. This latter contribution can be formulated as a convolution of the population activity with a kernel. Thus, it can simply be included in an extended version of the SPDE formalism presented here.
The stochastic-field approach to finite-size effects in neural network dynamics presented herein uses a surrogate quantity q(t, r) to access the population activity. With this approach, finite-size noise explicitly appears in the main equation, Eq 7, taking the form of a Gaussian white noise in the variables t and r, divided by ffiffiffiffi N p . This equation can be integrated numerically, providing information about the latent, refractory state of neurons. This information is typically not present in other approaches that deal with finite-size effects. Moreover, and most importantly, the firing activity can thus be directly simulated by using the SPDE; such direct simulation is typically not possible with other methods (but see below). On the other hand, the linear noise approximation gives analytical expressions for the power spectrum of finite-size activity fluctuations in the asynchronous state, a feat also possible by other approaches. For instance, the linear response theory applied to leaky integrate-and-fire networks readily generates this power spectrum [20]. Also, an integral equation has recently been proposed that directly involves the population activity [35]; linearizing this equation yields the power spectrum. Our preliminary investigations and Ref. [45] suggest that all these approaches lead to analogous results. Moreover, Schwalger et al. recently used an error minimizing procedure to directly describe the activity of finite-size networks at a so-called mesoscopic level. This procedure de facto sacrifices probability conservation to obtain a nevertheless accurate activity equation in the time domain. Our approach rather involves both time and refractoriness, which is still considered the microscopic level by these authors. However, working directly with our SPDE provides some advantages. First, even though we performed a Gaussian approximation of the Poisson spiking statistics, the probability is still conserved, which is the standard expectation for biophysically interpreted neural mass models [14]. Second, as stated above, its linearization in combination with the use of the method of characteristics enabled us to obtain a stochastic differential equation-along the characteristics-and to compute spectral and correlation properties of the network activity. Finally, we note that our approach can also be extended to multiple populations, and we pointed the reader in the right direction for doing so without performing the full calculations.

Derivation of the stochastic partial differential equation
We consider a homogeneous neural population and we denote the density of neurons in refractory state r at time t by q(t, r). The number of neurons with their refractory state in [r − Δr, r) at time t, denoted m(t, r), is given by At time t, a given neuron is in refractory state r if its last spike occurred at and survived without spiking until t. Accordingly, m(t, r) represents the number of neurons a time t whose last spiket belongs to interval (t − r, t − r + Δr]. By extension, m(t + Δt, r + Δr) represents the number of neurons whose last spike also occurred in (t − r, t − r + Δr], but that survived without spiking until t + Δt. Intuitively, m(t + Δt, r + Δr) is equal to m(t, r) minus the number of these neurons that did spike in [t, t + Δt). If Δt is small enough-so that the hazard function ρ(t, r) is nearly constant on this interval-then this number is a statistically independent Poisson random variable with mean and variance equal to ρ(t, r)m(t, r)Δt [25,35]. Therefore, with a two-dimensional Gaussian white noise process η(t, r) obeying hZðt; rÞZðt 0 ; r 0 Þi ¼ dðt À t 0 Þdðr À r 0 Þ:

Finite-size spiking neural networks
The initial condition needed to solve the SPDE is obtained as follows. The total number of neurons firing in interval [t, t + Δt), here denoted n(t), is given by where A(t 0 ) is the population activity. But n(t) is also the number of neurons at time t + Δt with refractory state in [0, Δt). Hence, In the limit of small Δt we get which is the required initial condition. Moreover, n(t) must be equal to the summation of all neurons that fire in [t, t + Δt), whatever their refractory state. Therefore, we must have, symbolically, Steps similar to those used above readily yield Finally, a normalization condition for q(t, r) is obtained from the fact that, at time t, all neurons within the network have a given refractory state: Eqs 26-28-and variants thereof-are used in the main text. Let us emphasize that we do not assume that the single neurons possess Poisson statistics. The Poisson random variable arises from an application of the theory of point processes. According to this theory, a generic point process with hazard function lðtjHðtÞ; YÞ-where H is the given history of events prior to t and Θ represents any covariate, e.g. an external stimulus-generates a spike in [t, t + Δt) with probability lðtjH; YÞDt þ OðDt 2 Þ [46]. A Poisson process is obtained when the hazard function does not depend on history, in which case the hazard function identifies with the mean field itself (i.e. the firing rate). In our case, the hazard function ρ depends on the history H t fAðt 0 Þ : 0 < t 0 < tg and the external input through the potential u(t, r). Fundamentally, the number of neurons firing in [t, t + Δt) follows a binomial distribution, which becomes a Poisson distribution in the limit Δt ! 0; the underlying Poisson random variable changes as a function of time and is correlated with other Poisson variables at other times (see appendix B in [35] for a mathematical explanation).

Linear noise approximation
We start with the system-size expansion given in Eqs 11 and 12, that we rewrite here for convenience AðtÞ Replacing these two equations in the SPDE, Eq 7, yields omitting function arguments for brevity. Here, d/dt is the material derivative operator d=dt @=@t þ @=@r: Using Eqs 30 and 4, we can write the hazard function as where r 0 ðt; rÞ r½u 0 ðt; rÞ with u 0 ðt; rÞ ðk Ã ½À J s A 0 þ I ext ÞðtÞ À VðrÞ and the derivative dρ/du is evaluated at u 0 . Then, matching terms of orders Oð1Þ and OðN À 1=2 Þ, and neglecting contributions from terms of order OðN À 1 Þ and lower, we get two equations involving q 0 and q ξ : The first of these equations is the usual mean-field model (Eq 13). The second equation is a SPDE whose coefficients, source and noise terms depend on the solution of the deterministic equation. Of course, the boundary condition on q(t, r) (cf. Eqs 8 or 26) implies Finite-size spiking neural networks Furthermore, from Eq 9 and the LNA, we have, after algebraic manipulations: Together, Eqs 32-34 constitute a coupled system between the deterministic mean field (A 0 and q 0 ) and the first-order fluctuations in the LNA (A ξ and q ξ ).

Activity in the asynchronous state
Let us denote by q 1 (r) and A 1 the refractory density and activity in the asynchronous steady state, respectively. We have to solve d dr q 1 ðrÞ ¼ À r 1 ðrÞq 1 ðrÞ; ð35Þ where Eq 35 can be integrated to give where we have used the boundary condition Finally, the average asynchronous activity can be computed using the conservation property of the neural network, namely Z 1 0 q 1 ðrÞdr ¼ 1; so that

ðsÞ ds dr:
Note that the mean firing rate is implicitly given since ρ 1 depends on A 1 . Therefore, to compute the mean activity, this nonlinear equation must be solved numerically. With our choices for ρ we can push the calculation further. We have where δu has been set to its numerical value of 1 mV, hence Z r 0 r 1 ðsÞds ¼ l 0 e h 1 ½r þ tðe À r=t À 1Þ Finite-size spiking neural networks and A À 1 The change of variable t ¼ l 0 te h 1 À r=t reduces the problem to finding the solution of the equation where the lower incomplete γ function is given by and we defined s 1 tr 1 ð1Þ ¼ tl 0 exp½I ext À J s A 1 : Eq 38 is then solved self-consistently for A 1 .

Stability analysis and the characteristic equation
To study the stability of the asynchronous state, one needs the eigenvalues of the differential operator once a linearization around the steady state has been performed. We therefore consider a small perturbation and write the solution in the form q 0 ðt; rÞ ¼ q 1 ðrÞ þ εq 1 ðt; rÞ þ Oðε 2 Þ; A 0 ðtÞ ¼ A 1 þ εA 1 ðtÞ þ Oðε 2 Þ: ð39Þ Plugging these expressions into Eq 13-keeping the first order terms only-yields the partial differential equation The eigenvalues are thus given by the roots of CðlÞ. The time-independent solution A 1 will be stable if all the eigenvalues have negative real parts.

Computation of the power spectrum
We assume that the deterministic part of the activity tends toward its equilibrium state-below the instability threshold. In other words, t is large enough that On the other hand, the stochastic part of the activity now reads (see second equation in Eq 34) The power spectrum is defined by is the Fourier transform of the stochastic process A ξ restricted to 0 < t < T. Analogously, we defineq and likewise for other stochastic processes depending on t and r. After taking the Fourier transform, Eq 48 becomes (arguments of functions are omitted when no confusion may arise) where we usedq x T ðo; 0Þ ¼Ã x T ðoÞ: Also, from Eq 49 we get ð51Þ Replacingq x T in this equation by the expression that has just been obtained for it yields which gives rise to four terms that we compute separately (we denote hÁi lim T!1 1 T hÁi). First term Second term  The power spectrum is then

Autocorrelation function
In this section, we compute A ξ and (a part of) its autocorrelation function. According to the Wiener-Khintchin theorem, the autocorrelation function of A ξ is simply the Fourier transform of its power spectrum. But to allow comparison with the results in [35], we calculate a part of this autocorrelation directly. From Eq 49, we must determine q ξ (t, r). This will be done by first solving Eq 48 using the method of characteristics. Solving the first-order inhomogeneous PDE with the method of characteristics. We have to solve The ranges of r and t are [0, 1) and (−1, 1), respectively. The boundary condition is q x ðt; 0Þ ¼ A x ðtÞ: With the method of characteristics, we first transform the PDE above into an ordinary Finite-size spiking neural networks differential equation. We find a family of curves for which Along these lines, Hence, we now have to solve the ODE on the characteristic curve, namely The solution is Fðt þ t 0 ; tÞ S 1 ðtÞ dt: Replacing we get the solution for the whole (t, r) space: In the case where the hazard function is exponential, as in Eq 5, a 1 becomes a 1 ¼ 1 du

Numerical implementation
The SPDE of Eq 7 can be readily integrated using the following numerical scheme: qðt þ Dt; r þ DtÞ ¼ qðt; rÞexpfÀ r½uðt; rÞDtg À ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi qðt; rÞH½qðt; rÞf1 À expðÀ r½uðt; rÞDtÞg NDt where H(x) is the Heaviside step function, and N ð0; 1Þ is a standard normal random variable. This numerical scheme is obtained by discretizing the time evolution of q(t, r) on the characteristic curve (cf. Methods) and noting that 1 À rqdt % e À rqdt ; with e −ρdt the probability that a spike is not fired during interval dt [33]. The Heaviside function prevents negative values for q(t, r) from appearing under the square-root sign. Along the characteristic curve, the dynamics correspond to a Cox-Ingersoll-Ross stochastic differential equation [47]. The proposed numerical scheme above is thus well defined [48] and produces results in excellent agreement with simulations of the full network (see Fig 4). One way to extract the population activity is to evolve q(t, r) according to the above numerical scheme for all r ¼ iDt; 8i > 0; and compute q(t, 0)-and thus A(t)-by enforcing the conservation law (Eq 10).
Finite-size spiking neural networks