## Abstract

There has recently been a great deal of interest in inferring network connectivity from the spike trains in populations of neurons. One class of useful models that can be fit easily to spiking data is based on generalized linear point process models from statistics. Once the parameters for these models are fit, the analyst is left with a nonlinear spiking network model with delays, which in general may be very difficult to understand analytically. Here we develop mean-field methods for approximating the stimulus-driven firing rates (in both the time-varying and steady-state cases), auto- and cross-correlations, and stimulus-dependent filtering properties of these networks. These approximations are valid when the contributions of individual network coupling terms are small and, hence, the total input to a neuron is approximately gaussian. These approximations lead to deterministic ordinary differential equations that are much easier to solve and analyze than direct Monte Carlo simulation of the network activity. These approximations also provide an analytical way to evaluate the linear input-output filter of neurons and how the filters are modulated by network interactions and some stimulus feature. Finally, in the case of strong refractory effects, the mean-field approximations in the generalized linear model become inaccurate; therefore, we introduce a model that captures strong refractoriness, retains all of the easy fitting properties of the standard generalized linear model, and leads to much more accurate approximations of mean firing rates and cross-correlations that retain fine temporal behaviors.

## 1. Introduction

One of the fundamental problems in the statistical analysis of neural data is to infer the connectivity and stimulus dependence of a network of neurons given an observed sequence of extracellular spike times in the network (Brown, Kass, & Mitra, 2004). Estimating the parameters of network models from data is in general a computationally difficult problem because of the high dimensionality of the parameter space and the possible complexity of the model's objective function given an observed data set of network spike times.

One class of models that has proven quite useful is known in the statistics
literature as the *generalized linear model* (GLM; Brillinger, 1988; Chornoboy, Schramm, & Karr, 1988; see also Martignon et al., 2000; Kulkarni & Paninski, 2007; Nykamp, 2007), variants of which have been applied successfully in the hippocampus
(Harris, Csicsvari, Hirase, Dragoi, & Buzsaki, 2003; Okatan, Wilson, & Brown, 2005), motor cortex (Paninski, Fellows, Shoham, Hatsopoulos, &
Donoghue, 2004; Truccolo, Eden, Fellows,
Donoghue, & Brown, 2005), retina (Pillow
et al., 2008), and cultured cortical slice
(Rigat, de Gunst, & van Pelt, 2006). In
its simplest form, this model incorporates both stimulus-dependence terms and direct
coupling terms between each observed neuron in the network; fitting the model
parameters leads to an inhomogeneous (stimulus-driven), nonlinear coupled spiking
model with delay terms. Statistically speaking, the model is attractive due to its
explicitly probabilistic nature and because fitting the model parameters is
surprisingly simple: under certain simple conditions, the log-likelihood function
with respect to the model parameters is concave, and the maximum likelihood
estimation of the parameters can be easily performed using standard ascent methods
(Paninski, 2004; Paninski, Pillow, &
Lewi, 2007).

Once the model parameters are obtained, we are left with an obvious question: What do we do next? One of the key applications of such a network model is to better understand the input-output properties of the network. For example, we would like to be able to predict the mean firing response of the network given a novel input and to carve out the impact of the network coupling terms on this stimulus-response relationship (e.g., what is local inhibition's impact on the stimulus filtering properties of the network?). We would also like to know how the correlation properties of spike trains in the network might depend on the stimulus, and in general how correlations might encode stimulus information (Riehle, Grün, Diesmann, & Aertsen, 1997; Hatsopoulos, Ojakangas, Paninski, & Donoghue, 1998; Oram, Hatsopoulos, Richmond, & Donoghue, 2001; Nirenberg, Carcieri, Jacobs, & Latham, 2002; Schnitzer & Meister, 2003; Kohn & Smith, 2005; Pillow et al., 2008).

We can in general study these questions with direct Monte Carlo simulations. However, simulation of a large-scale probabilistic spiking network is computationally expensive, since we often need to draw many samples (i.e., run many simulated trials) in order to compute the quantities of interest to the desired precision. More importantly, direct numerical simulation often provides limited analytical insight into the mechanisms underlying the observed phenomena.

The goal of this study is to investigate how much of the behaviors of these GLM networks can be understood using standard analytical mean-field approximations (Renart, Brunel, & Wang, 2003; Wilson & Cowan, 1972; Amit & Tsodyks, 1991; Ginzburg & Sompolinsky, 1994; Hertz, Krogh, & Palmer, 1991; Kappen & Spanjers, 2000; Gerstner & Kistler, 2002; Meyer & van Vreeswijk, 2002). In particular, we develop approximations for the mean firing rates of the network given novel stimuli, as well as the auto- and cross-correlation and input-output filtering properties of these networks. These approximations are valid when the contributions of individual network coupling terms are small and lead to deterministic ordinary differential equations that are much easier to solve and analyze than direct Monte Carlo simulation of the network activity.

However, in the case of strong refractory effects, these mean-field approximations become inaccurate, since the spike history terms in the generalized linear model must be large to induce strong refractoriness, and this pushes our approximations beyond their region of accuracy. Therefore we introduce a modified model: a generalized linear model with Markovian refractoriness. This model has several advantages in this setting: it captures strong refractoriness, retains all of the easy-fitting properties of the standard generalized linear model, and leads to much more accurate mean-field approximations.

## 2. The Generalized Linear Point Process Model

*N*recurrently connected spiking neurons, modeled as a multivariate point process (Snyder & Miller, 1991). The firing rate (conditional intensity function) of neuron

*i*at time

*t*is a function of the total input,

*u*(

_{i}*t*), the neuron receives: where

*f*(.) is a smooth, nonnegative, and monotonically increasing function. For numerical simulations, we use exponential nonlinarity,

*f*(

*u*) =

*e*, but the analysis in this article is not limited to this exponential nonlinearity. The total input to the neuron is expressed as the sum of external input

^{u}*I*and recurrent input

_{i}*H*: The recurrent input is modeled as a sum of identical responses caused by each presynaptic spike and written as where

_{i}*w*(

^{j}_{i}*t*) is a coupling term from neuron

*j*(upper index) to

*i*(lower index), and

*t*

_{j,n}is the

*n*th spike time of neuron

*j*. Without loss of generality, we may express each coupling term

*w*(

^{j}_{i}*t*) as a weighted sum of exponential functions of different time constants. As we will see below, this sum-of-exponentials representation is extremely useful: we will take full advantage of the Markovian nature of the decaying exponential function. To keep notation under control, we will restrict our attention to the case of a single exponential function with time constant τ

_{i}

^{j}, where Θ is the Heaviside step function, which takes one for positive arguments and zero otherwise; generalizations to the case where

*w*(

^{j}_{i}*t*) are given by a sum of exponentials with different time constants will be left implicit and will be straightforward in all cases. If we write the output spike train of neuron

*i*as the recurrent input of equation 2.3 is rewritten as the convolution of the synaptic filter and the presynaptic spike train: with This model is a generalization of the inhomogeneous Poisson process, since past output spikes modulate the firing intensity function. If we set all the coupling weights

*J*to zero, we recover the the inhomogeneous Poisson process model. The model is also closely related to the spike-response model of Gerstner and Kistler (2002) and the inhomogeneous Markov interval process of Kass and Ventura (2001), in which the firing rate depends on only the last observed spike. Note that in the GLM, the firing intensity of a neuron may depend not only on the last spike of the neuron but on the past several spikes of this neuron (see, e.g., Paninski, 2004; Paninski et al., 2007, for further discussion). As mentioned in section 1, this model is attractive because it is fairly easy to fit to data and seems to capture the firing statistics of neural populations in a variety of experimental settings. However, in this article, we will not discuss the estimation problem but rather limit our attention to the problem of determining the network firing rate, auto- and cross-correlation functions, and input-output properties, given fixed, known model parameters.

^{j}_{i}## 3. Master Equation

*P*(

**,**

*h**t*): this is the distribution of the state variable at time

*t*. These

*h*are the solutions of the multivariate coupled Markov process, and therefore solving for

^{j}_{i}*P*(

**,**

*h**t*) will allow us to compute a number of important quantities. For example, to compute the mean firing intensity ν

_{i}(

*t*) =

*E*[

*S*(

_{i}*t*)] of neuron

*i*given some input {

*I*(

_{i}*t*)}, we may compute where

*E*[.] is the average over the spike history in an entire duration

*T*given external input {

*I*(

_{i}*t*)}. Hence, we can read the mean firing rate off directly from

*P*(

**,**

*h**t*).

*e*_{j}denoting the vector with its slot δ

_{jj′}. Note that there are no diffusion terms here; the random nature of the stochastic differential equation enters instead through the jump terms.

While this jump-advection partial differential equation (PDE) in principle completely
describes the behavior of this system, unfortunately this equation remains difficult
to handle analytically (due to the nonlinearity in the λ terms) and
numerically (due to the large dimensionality of the state variable ** h**). We will pursue a direct analysis of this PDE elsewhere; here, we will
attack this system using a different approach, based on a self-consistent expansion
of the first and second moments of

*S*(

_{i}*t*).

*h*(.)), however, we may easily solve the PDE numerically to obtain an exact solution for the mean firing rate. Figure 2 compares this exact solution versus the mean-field approximations derived in the following section. The mean-field approximation is reasonably good here but not perfect, because of the fluctuation in the recurrent input (

*J*= −1, τ = 10 ms).

## 4. Mean-Field Approximation of GLM Statistics

Now we turn to the main topic of this article. As emphasized above, we would like to
calculate the mean-firing intensities of neurons, ν_{i}(*t*) = *E*[*S _{i}*(

*t*)], or the spike-cross-correlation function between neurons, . The calculation of these statistics is difficult, however, because of the nonlinearity

*f*and the recurrent input term

*H*in the model. A typical way of approximating these recurrent effects is to replace the interaction term

*H*by its mean field (Renart et al., 2003; Wilson & Cowan, 1972; Amit & Tsodyks, 1991; Ginzburg & Sompolinsky, 1994; Hertz et al., 1991; Kappen & Spanjers, 2000; Gerstner & Kistler, 2002; Meyer & van Vreeswijk, 2002).

Let us assume that each neuron receives input from some subset, *N _{c}*, of the total population of neurons

*N*. If the synaptic strength scales with

*N*, and assuming that the population is in the “asynchronous” state, then it is known that the cross-covariance function between two neurons scales with 1/

_{c}*N*for large

_{c}*N*(Ginzburg & Sompolinsky, 1994). This means that the contribution of each network interaction term is negligible and all the neuron activities are independent in the

_{c}*N*→ ∞ limit. Hence, calculating the mean-firing intensities and cross-correlations is easy in this limit.

_{c}We first derive the mean-firing intensity in the limit of *N _{c}* → ∞ and then evaluate the finite-size effect of order
1/

*N*. Note that this mean-field solution is a good approximation not only when

_{c}*N*is very large; equivalently, this is a good approximation when the synaptic coupling

_{c}*J*is small for fixed

*N*, where the contribution of an individual synapse is small.

_{c}*H*is roughly gaussian. We can express the mean firing intensities of the neurons in our network as a function of two sufficient state variables: the mean recurrent input

*E*[

*h*] and the variance of the recurrent input Var[

*H*], where, in the third line, the distribution of the recurrent input,

*H*(

_{i}*t*), is assumed to be gaussian with mean μ

_{i}(

*t*) =

*E*[

*H*(

_{i}*t*)] and variance , and the expectation was replaced by a gaussian integral: Note for an exponential nonlinearity

*f*(

*u*) =

*e*, the gaussian integral is given by

^{u}*F*(μ, σ) = exp (μ + σ

^{2/2}).

*t*>

*t*′ (the spike variable

*S*(

_{i}*t*) is averaged given past spike history to yield λ

_{i}(

*t*)), and the third term is the correlation for

*t*′ >

*t*. In the third line, we introduced conditional firing intensity, which describes the probability of neuron

*i*firing at

*t*given the spike of

*j*at

*t*′. We also introduce an abbreviation in the following and write

*E*[.|

*S*(

_{}j*t*′)] instead of

*E*[.|spike of

*j*at

*t*′]. Similar to the approximation in equation 4.1, we assume a gaussian distribution of

*H*given a spike of neuron

*j*at time

*t*′ to find for

*t*>

*t*′ with conditional mean, , and conditional variance, . Note that equation 4.5 is valid only for

*t*>

*t*′; if

*t*<

*t*′ instead, we evaluate and use the Bayes theorem to calculate . The unconditional and conditional mean and variance of

*H*in equations 4.1 and 4.5 are evaluated as where we defined with . The evolution of , and

*A*in equation 4.7 follows simple ordinary differential equations (ODEs) and is easy to simulate, where and . Note that the differentiation of

^{kl}_{i}*A*(

^{kl}_{i}*t*) for

*k*=

*l*is a little tricky because has a delta peak at

*s*=

*s*′. While the definition of

*A*(

^{kl}_{i}*t*) in equation 4.7 is clear, we have to make sure in the evaluation of not to double-count the integration of the delta peak. Hence, in the third equation of equation 4.8, we included the delta peak along the integration over

*s*(including the peak at

*s*=

*t*) but excluded this peak from the integration along

*s*′ (excluding the peak at

*s*′ =

*t*).

### 4.1. Solution in the Mean-Field Limit.

_{i}, and cross-correlation functions, ϕ

_{ij}, in the mean-field limit of

*N*→ ∞ (which is also valid in the noncoupled limit

_{c}*J*→ 0 for finite

*N*), that is, synaptic strengths are scaled

_{c}*J ∼*1/

*N*so that the individual contribution of a spike is negligible for large

_{c}*N*; and we also assume, as discussed above, that the network is in the asynchronous state. Under these conditions, the cross-covariances between neurons are known to scale as

_{c}*Δ ϕ*

_{ij}∼ 1/

*N*for

_{c}*i*≠

*j*and

*Δ ϕ*

_{ij}∼ 1 for

*i*=

*j*(Ginzburg & Sompolinsky, 1994). Then we can easily see σ

_{i}

^{2}(

*t*)∼ 1/

*N*because

_{c}*A*1 for

_{i}^{kl}∼*k*=

*l*and

*A*1/

_{i}^{kl}∼*N*for

_{c}*k*≠ l from equation 4.7. Hence, in the limit of

*N*→ ∞, we obtain from equations 4.1 and 4.8, This self-consistent update rule corresponds to a well-known mean-field dynamics of the mean-firing intensity (Ginzburg & Sompolinsky, 1994). Similarly, the scaling of cross-covariance

_{c}*Δ ϕ*

_{ij}∼ 1/

*N*for

_{c}*i*≠

*j*implies for

*k*≠

*j*and for

*k*=

*j*from equation 4.7. Moreover the conditional covariance scales with 1/

*N*for

_{c}*k*≠

*l*≠

*j*and thus . Hence, in the limit of

*N*→ ∞, we find for

_{c}*t*>

*t*′ and, from equation 4.3, These equations are very easy to simulate, and the approximation of the mean-firing intensity is good even for a small

*N*if

_{c}*J*is small (see Figure 2A). On the other hand, the mean-field approximation is not good for large

*J*(see Figures 3 and 4B). This is because, for fixed

*N*, the fluctuation of the recurrent input increases with

_{c}*J*, making the mean-field assumption (which neglects these fluctuations) invalid. Figure 4 shows the cross-correlation functions between two neurons. Because the self-interaction term

*J ∼*1/

*N*vanishes for the large

_{c}*N*limit, the autocorrelation function in this limit can capture the Poisson peak only at

_{c}*t*=

*t*′; it does not capture any nontrivial correlation caused by the interaction term

*J*for finite

*N*(see Figure 4).

_{c}### 4.2. Estimating the Finite Size Effect.

^{i}(

*t*), and the cross-covariance function, ϕ

_{ij}(

*t*,

*t*′), in the

*N*→ ∞ limit. Ideally, once provided the true term, which is of order 1/

_{c}*N*

^{2}

_{c}different from σ

_{i}

^{2}(

*t*), we can iteratively evaluate equations 4.1 to 4.8, starting from the mean-field limit solution of equations 4.9 and 4.10 to find the precise approximation of ν

^{i}(

*t*) and ϕ

_{ij}(

*t*,

*t*′) as long as the gaussian

*H*assumption is valid. However, we need the third-order correlations to evaluate , and the equations are not closed. Here, we give up evaluating

*O*(1/

*N*

^{2}

_{c}) terms and just evaluate the finite size corrections of ν

^{i}(

*t*) and ϕ

_{ij}(

*t*,

*t*′) up to

*O*(1/

*N*) terms. Up to this order, it can be argued that , and under this approximation, we obtain a self-consistent expression for ν

_{c}^{i}(

*t*) and ϕ

_{ij}(

*t*,

*t*′). Equations 4.1 and 4.5 yield, up to the first order of 1/

*N*, where to the first order of 1/

_{c}*N*. Equation 4.11 can be evaluated using equations 4.6 to 4.8, which we repeat here for convenience: In principle, the above equations are valid even with time-dependent input, and several iterations of the above set of equations yield time-dependent cross-correlation functions. However, evaluation of the above equations is computationally expensive with time-dependent input. If one neuron is connected to

_{c}*N*surrounding neurons and time is discretized into

_{c}*T*time bins, it takes

*O*(

*N*

^{2}⋅

*N*⋅

_{c}*T*

^{2}) operations and memory to evaluate . If the input is constant in time, is a function of

*t*−

*t*′ only, and just

*O*(

*N*

^{2}⋅

*N*

_{c}⋅

*T*) operations are required to solve the equations. Therefore, we focus on calculating the finite size effect for constant input in this article.

^{1}Generally several iterations of the above set of equations are required for the convergence of mean firing intensities and cross-correlations. However, by initially setting those variables to the solution in the limit of

*N*→ ∞, we found good approximations after a couple of iterations in many cases.

_{c}The mean firing rate and cross-correlation functions calculated from equations 4.11 and 4.12 are precise for small *J* and large *N _{c}* because the recurrent input is close to gaussian in this case (see
Figure 4A). On the other hand, for large

*J*and small

*N*, the gaussian input approximation is poor, and therefore the approximation becomes worse (see Figure 4B). One key result here is that both the autocorrelations and the cross-correlations in these networks depend on the inputs

_{c}*I*(

*t*) in general, as can be seen from equations 4.6 to 4.8. Changing just the baseline (mean) of the input

*I*(

*t*) can significantly change the correlation structure in the network as seen in a number of physiological preparations (e.g., Hatsopoulos et al., 1998; Kohn & Smith, 2005), even if the coupling terms

*J*and

*w*are unchanged.

### 4.3. Linear Input-Output Filter.

The ODE derived in equation 4.9 is general and may be applied given arbitrary time-dependent input
{*I _{i}*(

*t*)}. In this section, we exploit these equations to explore the linear response properties of these spiking networks. A standard physiological technique is white noise analysis (Marmarelis & Marmarelis, 1978; Rieke, Warland, de Ruyter van Steveninck, & Bialek, 1997), in which we inject white noise input in order to study a network's input-output filtering properties. When this analysis is applied to our GLM system, the external input is given by , where

*K*is neuron

_{i}*i*'s linear stimulus filter and ξ

_{i}(

*t*) is the white noise stimulus with mean and small variance . We keep the input variance small here in order to treat the input as almost constant and develop a consistent linear response theory. Note that 〈.〉 describes the average over the fluctuating stimulus ξ

_{i}(

*t*) here. We use

*δ X*=

*X*− 〈X〉 for any function

*X*of input stimulus and neglect

*O*(

*δ I*

^{2}) terms in the following calculations.

*N*→ ∞) limit. The mean-firing intensity of a neuron is given by the following self-consistent equations: In particular, for a constant stimulus, ξ

_{c}_{i}(

*t*) = ξ

_{i}

^{(0)}, we obtain from equation 2.4. Note that

*X*

^{(0)}describes responses at constant stimulus ξ

_{i}

^{(0)}, for example, . We can solve the zeroth-order terms in equation 4.14 self-consistently by finding the roots of the equations . In Figure 5, we compare this self-consistent mean-field solution with numerical simulation of a single neuron as a function of constant input,

*I*

^{(0)}, and self-inhibition strength

*J*. We see that the approximation is good for small values of self-interaction

*J*, but that the error increases with ∣

*J*∣ and

*I*

^{(0)}, as expected given the analysis in the previous section.

*δ ξ*, yields Now we can rewrite equation 4.15 by introducing vector notation as where we defined and . Therefore, the Fourier transformation of equation 4.16 yields with a gain function . Note that

**is an identity matrix here (not the input vector). We should note that a mathematically identical linear response filter was derived for a related model in Ginzburg and Sompolinsky (1994). One difference is that in their model, spinlike binary variables flip depending on input level with a fixed time constant. This corresponds, in our model, to the case that all the synaptic time constants τ**

*I*_{i}

^{j}are identical. Hence, the linear response filter in equation 4.17 is a natural generalization of their result.

*K*is set to be a single exponential function with time constant τ

_{I}, for simplicity, the Fourier transformations of the synaptic filter and the input filter are and , respectively. In this case, we obtain, after a Fourier inverse transformation by a Cauchy integral, where and . We see that the linear response may be described as a combination of the stimulus-filter term

*K*along with an additional term due to the recurrent filter

*w*. Of course, in the limit

*J*→ 0, we recover

*G*→

*f*′

*K*. Note also that the input-output filter changes with the baseline input through

*f*′. In the case of the exponential nonlinearity,

*f*, and

*J*<0, the larger the input baseline, the sharper the linear filter. This is because the effective strength of of the spike history effect is modulated by the input-output nonlinearity and the baseline input level (see, e.g., Bryant & Segundo, 1976, for some similar effects in an in vitro physiological preparation).

*z, t*) with unit variance. Input to the neuron

*i*is described by , where are spatial filters with means

*z*

_{1}= 1 and

*z*

_{2}= −1. For simplicity, we assume that two neurons receive the same baseline of input,

*I*

^{(0)}

_{1}=

*I*

^{(0)}

_{2}=

*I*

^{(0)}with an identical temporal stimulus filter, , and have symmetric synaptic interactions

*J*

_{12}=

*J*

_{21}=

*J*= −1,

*J*

_{11}=

*J*

_{22}= 0, and . The derivative of the nonlinear function is described by for

*i*= 1, 2. Proceeding as above, we can calculate the spatiotemporal input-output filter of neuron

*i*as . Similar to the previous single-neuron case, the Fourier transformation of the temporal stimulus filter is , while the synaptic filter is , which yield the Fourier transformation of the temporal input-output filter, After a Fourier inverse transformation using Cauchy integration, we find Figure 6B compares the analytically derived input-output filter of equation 4.20 with numerical simulation for different input baseline

*I*

^{(0)}= 2 and

*I*

^{(0)}= 4. Figure 6B, also shows the spatiotemporal input-output filter of neuron 1, , for the two input baselines. Under large input baseline conditions, we see stronger interactions between the two neurons and, hence, significant differences in the effective input-output filtering properties of the neuron in the two conditions. Thus, as emphasized above, the filtering properties of the network can change in a stimulus-dependent way, without any changes in the system parameters

*K*,

*J*, or

*w*. (See, e.g., Pillow et al., 2008, for some related results concerning the effects of interneuronal coupling terms on the filtering properties of populations of primate retinal ganglion neurons.)

## 5. The Generalized Linear Point Process Model with Markov Refractoriness

In section 4, we saw that we could derive simple
approximate expressions for many important quantities (firing rate, correlations) in
the basic generalized linear point process model as long as the individual synaptic
coupling terms ∣ *J* ∣ are not so large that the
accuracy of the mean-field expansion is compromised (see Figures 3 and 4).

In particular, the small ∣ *J* ∣ condition is likely
acceptable for the multineuronal coupling terms, since in many cases, the inferred
coupling parameters—in motor cortex (Paninski et al., 2004; Truccolo et al., 2005) and retina (Pillow et al., 2008), for example—have been empirically found to be small
(though larger network effects are found in hippocampus; see Harris et al., 2003; Okatan et al., 2005). Similarly, we might expect the
“slow” self-inhibition terms in the GLM (the terms responsible for
adaptation in the firing rate, for example) to be relatively small as well (Pillow,
Paninski, Uzzell, Simoncelli, & Chichilnisky, 2005; Truccolo et al., 2005).

However, it is clear that the generalized linear model with small weights ∣ *J* ∣ will not be able to account for the strong
refractoriness that is a fundamental feature of neural point-process dynamics at
fine timescales: in the GL model, strong, brief inhibition is produced by large,
negatively weighted, and sharply decaying history effects *w ^{i}_{i}*(

*t*), and these large

*w*(

^{i}_{i}*t*) terms spoil our mean-field approximation. Thus, it is clear that we need to adapt the methods introduced above in order to handle strong refractory effects.

One possible solution would be to take a discrete-time approach: instead of modeling
neural responses as point processes in continuous time, we could simply model each
cell's spike count within a discrete-time bin as a Bernoulli random variable whose
rate ν^{i}(*t*) is an increasing function of *u _{i}*(

*t*), very much as in the continuous-time model described above (Chornoboy et al., 1988; Okatan et al., 2005). This discrete-time Bernoulli model shares much of the ease of fitting (including concavity of the log likelihood; Escola & Paninski, 2007) as the continuous-time model. In fact, the discrete-time model converges in a natural sense to the continuous-time model as the time binwidth

*dt*becomes small. The advantage of this discrete-time formalism is that the firing rate ν

^{i}(

*t*) can never exceed 1/

*dt*(since the Bernoulli probabilities cannot exceed one), and therefore, in a crude sense, the discrete-time model has a refractory effect of length

*dt*. A mean-field analysis can be developed for this discrete-time model, which exactly mirrors that developed above in the continuous-time model. We simply need to exchange our ordinary differential equations for discrete-time difference equations.

*f*(.) and

*u*(

_{i}*t*) are defined as in section 4 and

*I*(.) is the indicator function that takes the value 1 if the argument is true. We have introduced an auxiliary refractory variable

*x*(

_{i}*t*), which takes a discrete state from 1 to

*M*. We assume that this variable

*x*(

_{i}*t*) is itself Markovian, with transition rates that is, when

*x*(

_{i}*t*) is in state

*M*, it will transition to state 1 with each spike, and then transitions occur from state

*m*to state

*m*+ 1 with rate 1/τ

_{r}until

*x*(

_{i}*t*) has reached the spiking state

*M*once again. Refractoriness in this model is enforced because the neuron is silent whenever

*x*(

_{i}*t*) is in one of the postspike quiescent states

*x*(

_{i}*t*)<

*M*. It is easy to see that this is a kind of inhomogeneous renewal model (the delay τ

_{d}required for the state variable

*x*to move from state 1 to first reach the active state

_{i}*M*is a sum of (

*M*−1) independent and identically (i.i.d.) distributed exponential random variables of mean τ

_{r}, and so τ

_{d}is an i.i.d. gamma variable with parameters (

*M*− 1, 1/τ

_{r})) and may be considered a special case of the “inhomogeneous Markov interval” model introduced by Kass and Ventura (2001). By adjusting the number of states

*M*and the rate 1/τ

_{r}, we may adjust the refractory properties of the model: for example,

*M*= 1 implies that we have an inhomogeneous Poisson model, while

*M*= 2 provides an exponentially decaying relative refractory effect, and if we let

*M*be large, with τ

_{r}scaling like τ

_{r}

*∼ τ*/

*M*, we obtain an absolute refractory effect of length τ.

Finally, it is easy to show that the additional Markov term in the definition of the
rate λ_{i}(*t*) does not have a negative impact on the estimation of the
GL model parameters *J* given spiking data; as discussed in Paninski
(2004), maximizing the likelihood in this
model (or constructing an EM algorithm, as in Escola & Paninski, 2007) requires that we maximize a nonnegatively weighted
version of the standard point-process log likelihood, and this weighted log
likelihood retains all of the concavity properties of the original GL model. Thus,
this new model requires just a single concave optimization and is as easy to fit as
the standard GL point process model. The advantage is that a strong relative
refractory effect is intrinsic to the model and does not need to be enforced by a
large, brief, negative *w ^{i}_{i}*(

*t*) term. Instead, we may use small adjustments to the

*J*terms to fine-tune the model to match the short-time details of the observed interspike interval density.

^{i}_{i}^{2}As we will see below, this leads to much more accurate mean-field approximations of the firing rates and correlation functions in these networks.

*i*is in state

*x*(

_{i}*t*) = 1,…,

*M*. Note that the boldface characters denote vectors. The dynamics of

*p*_{i}is described by where has the same elements as

*W*(

_{i}*t*) of equation 5.2 except for which describes the probability of neuron

*i*firing at time

*t*given

*x*(

_{i}*t*) =

*M*. Note that we used and wrote the last component of

*p*_{i}(

*t*) as [

*p*_{i}(

*t*)]

_{M}=

*p*(

_{i}*t*).

*i*being in state

*m*given a spike of neuron

*j*at time

*t*′. In particular, the last component is written as . The evolution of is described by a similar equation, where the transition matrix has the same elements as

*W*(

_{i}*t*) of equation 5.2 except for

Before turning to “mean-field” approximation, it is worth noting that
we may solve for the firing rates and correlations exactly in the special case of no
synaptic couplings, *J* = 0. (The analogous case in the
standard GL model is the inhomogeneous Poisson case, which is of course trivial.) We
can proceed by exploiting either the renewal nature of the model (Gerstner &
Kistler, 2002), which leads to convolution or
infinite-sum formulas for the firing rate, or the Markovian structure, which leads
to somewhat more intuitive ordinary differential equations. We pursue the second
approach here.

*I*, it is easy to solve the fixed point of

_{i}

*p*_{i}(

*t*); we find for the last component that .

*J*= 0 case, since the cross-coupling terms are set to zero.) In case of no couplings,

*J*= 0, the conditional firing intensity is given by for

*t*>

*t*′. We find those corner elements of as for

*t*>

*t*′. This means that for

*t*>

*t*′. Hence, follows the same differential equation as

*p*_{i}(

*t*) in this

*J*= 0 case but starting from the initial condition because the state

*x*is reset to state 1 just after a given spike of neuron

_{i}*i*.

In the following sections, we apply mean-field approximations for the firing rates
and correlations in the nonzero *J* case.

### 5.1. Mean-Field Approximation of the GL Model with Markov Refractoriness.

*H*(

_{i}*t*) given

*x*(

_{i}*t*) =

*M*to find with conditional mean and conditional variance . Also, assuming a gaussian distribution of recurrent input

*H*(

_{i}*t*) given

*x*(

_{i}*t*) =

*M*and a spike of neuron

*j*at time

*t*′, we find for

*t*>

*t*′ with conditional mean and variance . The dynamics of

*p*(

_{i}*t*) and are described by equations 5.3 and 5.5, respectively.

### 5.2. Solution in the Mean-Field Limit.

*N*→ ∞, where

_{c}*J ∼*1/

*N*and assuming an asynchronous state so that

_{c}*Δ ϕ*

_{ij}∼ 1/

*N*. Because all the cross-covariance functions vanish and in the limit, we find, from equations 5.11 to 5.13, for

_{c}*t*>

*t*′, while the cross-correlation function is given by equation 4.3. For the corner terms of the transition matrices, we find from equation 5.4 and for

*t*>

*t*′ from equation 5.6. Hence, the evolution of

*p*_{i}(

*t*) and is described by for

*t*>

*t*′, where the initial condition of is given by equation 5.10. Finally, from equation 5.15, the evolution of is written as an ODE, where we can apply equation 5.10 to evaluate in this limit. Figure 7 plots the mean firing intensity calculated from the above equations for various input, and the mean-field approximation provides good approximation of the time dependent firing intensities. Figure 8 compares the above approximation (

*N*→ ∞ limit) with the mean-field approximation, including the finite size effect that we discuss in the next section. In the mean-field limit, all the cross-correlations vanish and the autocorrelation functions are approximated by including the Markov refractory effect but not the self-interaction terms

_{c}*J*. We can also calculate the linar input-output filter discussed in section 4.3 with this Markov refractoriness (see the appendix).

^{i}_{i}### 5.3. Estimating the Finite Size Effect.

Similar to the case without Markov refractoriness, we evaluate the mean firing
intensity and cross-covariance functions up to 1/*N _{c}* terms. We find, as in the case without Markov refractoriness, for

*i*≠

*j*, for

*i*≠

*k*, and the third-order covariances, such as , are of order 1/

*N*

^{2}

_{c}under the asynchronous state if the contributions of individual synaptic weights are small,

*J ∼*1/

*N*. This implies .

_{c}

*p*_{i}(

*t*) and are described by equations 5.3 and 5.5, respectively. We now want to calculate μ

_{i}(

*t*), , and σ

_{i}

^{2}(

*t*) to the first order of 1/

*N*. First, direct differentiation of in equation 5.15 yields an ODE update equation (see the appendix): with . The origin of the second term on the right-hand side of equation 5.19 is due to the correlation of spike variable

_{c}*S*(

_{k}*s*) and the refractory variable

*x*(

_{i}*t*). Next, we need to evaluate to the zeroth order of 1/

*N*for

_{c}*k*=

*i*or

*k*=

*j*because

*J ∼*1/

*N*makes the contribution of them ∼ 1/

_{c}*N*, and we need to evaluate to the first order of 1/

_{c}*N*for

_{c}*k*≠

*i*and

*k*≠

*j*case. Except for the case

*i*=

*j*=

*k*, we find (see the appendix) to the order described above. The precise evaluation of this term for

*i*=

*j*=

*k*is much harder but, to a good approximation, equation 5.20 holds (see the appendix). Note that the approximation of the

*i*=

*j*=

*k*term affects the order 1/

*N*term of the autocovariance function, , but affects only the

_{c}*O*(1/

*N*

^{2}

_{c}) terms of the mean firing intensities and cross-covariance functions.

*O*(1/

*N*

^{2}

_{c}) on the quantities in equation 5.13, we find (see the appendix) Altogether, using equation 4.12, we obtain the following update equations to evaluate the finite size effect up to order 1/

*N*: and the dynamics of

_{c}

*p*_{i}(

*t*) and follows equations 5.3 and 5.5. As is the case without the Markov refractoriness, calculating this finite size effect is computationally expensive for time-dependent input, although in principle, this equation should work under that case as well. We plot the mean firing intensities of two neurons and cross-correlations in Figure 8. The finite size effect well captures the cross-correlation functions. It also successfully captures a peak in the cross-correlation function if two neurons are not directly connected but receive a common input from a third neuron (see Figure 8). Again, if the values of the above equations are initially set to the solution in the limit of

*N*→ ∞ given by equations 5.16 to 5.18, only a couple of iterations of the above equations are sufficient for the evaluation of cross-correlation functions.

_{c}## 6. Discussion

We have introduced mean-field methods for analyzing the dynamics of a coupled
population of neurons whose activity may be well approximated as the output of a
generalized linear point process model. Our approximations for the mean time-varying
firing rate and correlations in the population are exact in the mean-field limit
(*N _{c}* → ∞ under

*J ∼*1/

*N*scaling), though we have found numerically that the finite size correction to the mean-field equations is useful even for physiologically relevant connectivity strengths. The approximations may be computed by solving a set of coupled ordinary differential equations. This approach is much more computationally tractable than direct Monte Carlo sampling and may lead to greater analytical insight into the behavior of the network. In addition, we have introduced a new model, the generalized linear point process model with Markovian refractoriness, that captures strong refractoriness, retains all of the easy-fitting properties of the standard generalized linear model, and whose firing rate dynamics are much more amenable to mean-field analysis than the standard modells. Note that most of our illustrations involved the simulation of only a couple of mutually connected neurons; this small-

_{c}*N*case can be considered a kind of worst-case analysis, since with more asynchronous neurons (larger

_{c}*N*), the central limit theorem becomes progressively applicable, the distribution of recurrent input approaches the assumed gaussian distribution, and our approximations become more accurate.

_{c}This kind of mean-field analysis—or more generally, approaches for reducing the complexity of the dynamics of networks of noisy, nonlinear neurons—has a long history in computational neuroscience, as reviewed (Hertz et al., 1991; Gerstner & Kistler, 2002; Renart et al., 2003). While the point-process models we have discussed here are somewhat distinct from the noisy integrate-and-fire-type models that have been analyzed most extensively in this literature (e.g., Mattia & Del Giudice, 2002; Shriki, Hansel, & Sompolinsky, 2003; Moreno-Bote & Parga, 2004, 2006; Fourcaud-Trocme, Hansel, van Vreeswijk, & Brunel, 2003; Chizhov & Graham, 2007; Doiron, Lindner, Longtin, Maler, & Bastian, 2004; Lindner, Doiron, & Longtin, 2005), it is more important to stress the differences in the motivation of this work and previous work. Our main goal here was to provide analytical tools that an experimentalist who has fit a point process model to his or her data (as in, e.g., Paninski, 2004; Truccolo et al., 2005; Okatan et al., 2005; Pillow et al., 2008) may use to understand the behavior of the network or single-cell models that have just been inferred. In particular, we were interested in predicting, for example, the mean firing rate of a specific GLM network to a novel arbitrary dynamic input stimulus. This contrasts with the literature cited above, which has for the most part focused on the mean firing rates, first-order response, and fixed-point stability properties behavior of idealized, infinitely large populations (Mattia & Del Giudice, 2002, is an exception here) with homogeneous membrane and connectivity properties. In addition, most analytical studies of the input-output property of model neurons with a dynamical input start from a Fokker-Planck formalism and rely on a direct numerical simulation of the partial differential equation (Nykamp & Tranchina, 2000; Knight, Omurtag, & Sirovich, 2000; Chizhov & Graham, 2007; this direct approach is feasible only in the case that the system dynamics may be reduced to a state space of dimension at most two or so, and therefore does not apply in the networks studied here), or on linear response theory that assumes a small time-dependent component relative to the baseline component (Shriki et al., 2003; Fourcaud-Trocme et al., 2003; Doiron et al., 2004; Lindner et al., 2005), or on some quasi-stationary assumption (Knight et al., 2000; Mattia & Del Giudice, 2002; Fourcaud-Trocme & Brunel, 2005; Moreno-Bote & Parga, 2004, 2006) restricting the applicable input stimulus to slowly changing one or two simple step stimuli with sufficiently long gaps between steps that the network may reach equilibrium. Thus, it is difficult to apply the insights developed in these previous analyses directly to obtain a quantitative prediction of the dynamic firing rates of a specific nonhomogeneous, nonsymmetric network.

We derived, in this article, the finite size correction to the mean-field equation
that captures cross-correlation functions between neurons. This finite size effect
originates from the x*J _{ij} ∼* 1/

*N*scaling of synaptic strengths that guarantees fluctuation of inputs under asynchronous states. Hence, the finite size effect described in this article is the contribution of small fluctuation in the input around the mean. Another kind of finite size effect has also been discussed in the literature. In a sparsely connected network, correlations between two neurons disappear if

_{c}*N*is sufficiently smaller than the total number of neurons,

_{c}*N*(Derrida, Gardner, & Zippelius, 1987). The finite

*N*effect has been evaluated using stochastic Fokker-Planck equations (Brunel & Hakim, 1999; Mattia & Del Giudice, 2002). One should note, however, that the evaluation of this finite size effect for an experimentally estimated network structure is generally not straightforward because the evaluation of a Fokker-Planck equation with many state-space variables is computationally hard. Finally,

_{c}*J*1/

_{ij}∼*N*is not the unique way to scale synapses. Under the balanced input assumption, yields order one fluctuation of input even in the limit of

_{c}*N*→ ∞ (Sompolinksy, Crisanti, & Sommers, 1988; van Vreeswijk & Sompolinsky, 1996, 1998). However, nontrivial cross-correlation functions cannot be captured under the standard asynchronous assumption of mean-field analysis in the

_{c}*N*→ ∞ limit. The calculation of the finite size effect of a network with is not the scope of this article.

_{c}Two works that bear a stronger mathematical resemblance to ours are Ginzburg and
Sompolinsky (1994) and Meyer and van
Vreeswijk (2002). These authors applied
mean-field approximations to coupled model neurons to study the mean firing rates
and autocorrelation and cross-correlation functions in the network. However, the
two-state neuron model discussed in Ginzburg and Sompolinsky (1994) (where the neurons flip between
“active” and “inactive” firing states according to a
Markov process with some finite rate constant) is somewhat distinct from the
point-process models we have treated here (where the neuron remains in the
“active” state for a negligible time—i.e., the neuron spikes
instantaneously), and we have not been able to translate these earlier results to
the problems of interest in this article. Renewal point spiking neuorns are analyzed
in Meyer and van Vreeswijk (2002), and their
results are closer to ours. However, according to their refractory model, one has to
consider infinitely many refractory states in the continuous time limit, whereas our
spiking neuron model has *M* refractory state; *M* can
be as small as 2 while keeping the order 1 strong refractory effect. Accordingly,
their cross-correlations should be evaluated by integral equations that are
computationally (and conceptually) somewhat more involved, while the methods
developed here require us to evaluate just a couple of differential equations.

Finally, the effect of common input on the cross-correlation function was previously
studied in a related model (Nykamp, 2007),
using an expansion of the output *f*(.) nonlinearity arround *J* = 0. One major difference between the analysis we have
presented here and this previous work is that a simple expansion around *J* = 0 does not lead to a good approximation of the
firing rate, even in the mean-field limit of *N _{c}* → ∞, because the recurrent input changes the baseline firing
rate (cf. Figure 5); thus, it is much more
accurate to expand around the zeroth-order firing rate given by the roots of
equation 4.13, as discussed in
section 4.3. (On the other hand, our
mean-field approach does require a gaussian approximation that is known to be
inaccurate in the case of large

*J*terms.) It would be interesting to explore whether the methods developed here could help lead to more accurate inference of the common-input effects discussed in Nykamp (2007).

There are many possible applications of this mean-field method for problems that require fast evaluation of mean firing rates and cross-correlation functions. One example is to evaluate the information coded by spiking of a recurrently connected network about the input stimulus. This kind of information calculation usually requires averaging over recent spike history (Toyoizumi, Aihara, & Amari, 2006), but the gaussian approximation of the input described here could greatly ease the computational complexity to evaluate the information of stimulus coded by the network. We hope to pursue this direction in future work.

## Appendix: Mathematical Details

In this appendix, we collect the details of the expansion analyses summarized in the main text.

### A.1. Calculation of , , and *B*^{kl}_{i}.

*B*

### A.2. Approximation of , and *B*^{kl}_{i}(*t*).

*B*

*N*terms of the mean firing intensity and cross-covariance functions—we need to evaluate , , and to the first order of 1/

_{c}*N*. Assuming that the synaptic strengths scale as

_{c}*J ∼*1/

*N*, we need to evaluate to the first order of 1/

_{c}*N*, to the zeroth order of 1/

_{c}*N*if

_{c}*k*=

*i*or

*k*=

*j*because

*J ∼*1/

*N*but to the first order of 1/

_{c}*N*otherwise (in the following, we divide these into five cases and consider each case separately), and to the first order of 1/

_{c}*N*if

_{c}*i*,

*j*,

*k*are all different, and to the zeroth order of 1/

*N*if

_{c}*k*=

*l*,

*k*=

*i*, or

*l*=

*i*.

Next, to evaluate , we use if *i*, *k*, *l* are all different; for *i* ≠ *k*. As we
discussed at the beginning of this section, in order to evaluate to the first order of 1/*N _{c}*, we have to consider the number of combinations of indexes. For
example, there are

*N*terms of and

_{c}*N*

^{2}

_{c}terms of . Hence, we consider the following five possibilities for combining indexes:

- •
- •
- •
- •
- •

*i*receives extremely strong input, the conditional firing intensity is almost zero if

*s*and

*t*′ are as close as the refractory time constant of the neuron (

*M*− 1)τ

_{r}, so the contribution of this interval on the integral is small. If

*s*and

*t*′ are separated more than the refractory time constant, the interaction between the spike at

*s*and

*t*′ becomes small, that is, , and becomes small. Hence, to a good approximation, The validity of this intuitive approximation should be examined by numerical simulations. Note that contributes to the order 1/

*N*term of the autocovariance function

_{c}*Δ ϕ*

_{ii}(

*t, t*′) but contributes only to order 1/

*N*

^{2}

_{c}terms of the mean firing intensity ν

^{i}(

*t*) and cross-covariance functions with

*i*≠

*j*.

*B*^{kl}

_{i}as where we used and . We need to evaluate

*B*^{kl}

_{i}up to the first order of 1/

*N*if

_{c}*i*,

*k*,

*l*are all different and up to the zeroth order of 1/

*N*if

_{c}*k*=

*l*,

*k*=

*i*, or

*l*=

*i*. Up to these orders,

*B*^{kl}

_{i}can be approximated as Note that for

*i*≠

*k*, if

*i*,

*k*,

*l*are different, and for

*i*≠

*k*.

### A.3. Calculation of the Linear Filter in the Generalized Linear Model with Markov Refractoriness.

*J*= 0. In this case, the linear response is described by The calculation of is straightforward thanks to the bidiagonal nature of . After direct matrix inversion, we find with and . Hence, the linear response is simplified as In particular, if

^{k}_{i}*M*= 2, we find that the Fourier inverse transformation of the gain function has an analytical form, that is, using , we find Note that in the limit of τ

_{r}→ 0, the trivial gain function is . Hence, we can see that the suppressive kernel is added to this instantaneous gain function due to the Markov refractoriness. This linear response for

*J*= 0 case is a special case of a more general linear response for a renewal neuron (Gerstner & Kistler, 2002). Note that in equation A.18, we discussed a more general case and considered the interaction between the refractory effect and spike interaction effect of a network of recurrently connected neurons.

^{k}_{i}## Acknowledgments

We thank E. Shea-Brown for helpful conversations. T.T. is supported by the Robert Leet and Clara Guthrie Patterson Trust Postdoctoral Fellowship, Bank of America, Trustee. L.P. is supported by an NSF CAREER award, an Alfred P. Sloan Research Fellowship, the McKnight Scholar award, and NEI grant EY018003.

## Notes

Note that the computational complexity to calculate the finite size effect is not *N*^{2} but *N*^{2} ⋅ *N _{c}* because each synaptic time constant,

*τ*, is distinct. Hence, we needed to evaluate each term separately. However, if all the synaptic constants onto each neuron are constant, that is,

_{i}^{k}*τ*= τ

_{i}^{k}_{i}, we can directly evaluate without evaluating

*α*separately, and the computational complexity to evaluate the finite size effect reduces to

_{ij}^{k}*O*(

*N*

^{2}⋅

*T*).

In this article, we mostly use *M* = 3 (except for Figure 7B, where the effect of abrupt
refractoriness is studied), because this is the smallest value of *M* that requires the following mean-field formulation (we do
not need a vector formulation for an *M* = 2 case because
there is only one degree of freedom that corresponds to the probability of being
in the active state). It is worth noting that Escola and Paninski (2007) discuss methods for additionally
estimating an optimal transition matrix *W _{i}*(

*t*) via an EM algorithm. This provides another method for adjusting the short-time details of the model's responses. In addition, we may extend many of the mean-field methods developed below to the case of more general rate matrices

*W*(

_{i}*t*). However, for simplicity, we will not pursue these extensions here.