## Abstract

The local field potential (LFP) captures different neural processes, including integrative synaptic dynamics that cannot be observed by measuring only the spiking activity of small populations. Therefore, investigating how LFP power is modulated by external stimuli can offer important insights into sensory neural representations. However, gaining such insight requires developing data-driven computational models that can identify and disambiguate the neural contributions to the LFP. Here, we investigated how networks of excitatory and inhibitory integrate-and-fire neurons responding to time-dependent inputs can be used to interpret sensory modulations of LFP spectra. We computed analytically from such models the LFP spectra and the information that they convey about input and used these analytical expressions to fit the model to LFPs recorded in V1 of anesthetized macaques (*Macaca mulatta*) during the presentation of color movies. Our expressions explain 60%–98% of the variance of the LFP spectrum shape and its dependency upon movie scenes and we achieved this with realistic values for the best-fit parameters. In particular, synaptic best-fit parameters were compatible with experimental measurements and the predictions of firing rates, based only on the fit of LFP data, correlated with the multiunit spike rate recorded from the same location. Moreover, the parameters characterizing the input to the network across different movie scenes correlated with cross-scene changes of several image features. Our findings suggest that analytical descriptions of spiking neuron networks may become a crucial tool for the interpretation of field recordings.

- data-driven models
- gamma oscillations
- neural coding
- primary visual cortex
- recurrent networks
- slow oscillations

## Introduction

The local field potential (LFP) is obtained by low-pass filtering extracellular field potentials and is a signal sensitive to key integrative synaptic processes that cannot be obtained by only measuring the spiking activity of small populations (Buzsáki et al., 2012). Investigating how LFPs are modulated by external stimuli can thus offer additional insights into neural representations of the external world beyond those offered by spike trains only (Buzsáki et al., 2012; Einevoll et al., 2013). However, the LFP captures multiple stimulus-dependent neuronal processes over a wide range of temporal scales, so it is inherently more ambiguous and difficult to interpret than spikes. The interpretation of experimentally recorded LFPs needs network models that can disambiguate the contributions from different neural processes and identify quantitatively the key aspects and parameters of network dynamics.

An example of the richness of network representation of stimuli revealed by LFPs comes from studies of visual cortex. The LFP power in primary visual cortex (V1) carries visual information in the gamma (50–100 Hz) and high-gamma (>100 Hz) bands (Henrie and Shapley, 2005; Belitski et al., 2008), as well as in the low-frequency (<12 Hz) range (Belitski et al., 2008). Oscillations in the gamma band are a common feature of networks composed of excitatory and inhibitory neurons (Wilson and Cowan, 1972; Tsodyks et al., 1997; Li 2001; Börgers and Kopell, 2003; Brunel and Wang, 2003). In networks of noisy and/or randomly connected spiking neurons, fast network oscillations are accompanied by highly irregular single-cell firing and the amplitude of such oscillations can be strongly modulated by external inputs (Brunel and Hakim, 1999; Brunel, 2000; Brunel and Wang, 2003). Moreover, such networks, when stimulated by a time-dependent inputs (mimicking thalamic inputs to V1 during naturalistic stimulation), can generate by entrainment low-frequency LFPs that reproduce well the stimulus information and the spike-LFP relationships observed experimentally (Mazzoni et al., 2008).

Here, we study the dynamics of a network of excitatory and inhibitory integrate-and-fire neurons in the presence of time-dependent inputs and compute a novel analytical expression of the LFP spectra and the information that they convey about stimuli as a function of both the intrinsic parameters (neuronal and synaptic) and of the extrinsic parameters characterizing the dynamic input to the network. Unlike numerical simulations, these analytical formulae permit a systematic fit of models to real data to estimate hidden network parameters from LFP recordings. We fitted these analytical expressions to experimental V1 LFPs recorded in response to movie stimulation (Belitski et al., 2008) and evaluated how well these models explain LFP spectra and the stimulus information they carry. We also estimated how well the fitted parameters predicted observed multiunit activity and how well the estimated cross-movie-scene changes in the parameters characterizing the input to the network correlated with cross-scene changes of movie image features.

## Materials and Methods

#### Ethics statement

The data analyzed here have already been analyzed in part in previous studies (Belitski et al., 2008; Belitski et al., 2010). Recordings were obtained from the visual cortex of adult rhesus monkeys (*Macaca mulatta*) using the procedures described herein. All procedures were approved by local authorities (Regierungspraesidium Tuebingen) and were in full compliance with the guidelines of the European Community (EUVD 86/609/EEC) and with the recommendations of the Weatherall report on the use of nonhuman primates in research (Weatherall, 2006). Before the experiments, a form-fitting head post and recording chamber were implanted under aseptic surgical conditions and general balanced anesthesia (Logothetis et al., 2010). As a prophylactic measure, antibiotics (enrofloxacin, Baytril) and analgesics (flunixin, Finadyne veterinarian) were administered for 35 d postoperatively. The animals were socially (group-) housed in an enriched environment under daily veterinary supervision and their weight and food and water intake were monitored.

#### Electrophysiological recordings

We analyzed five different recording sessions from V1 of three different adult male macaque monkeys. These recordings were obtained while the animals were anesthetized (remifentanyl, 1 μg/kg/min), muscle-relaxed (mivacurium, 5 mg/kg/h), and ventilated (end-tidal CO_{2} 33 mmHg, oxygen saturation >95%). Body temperature was kept constant and lactated Ringer's solution supplied (10 ml/kg/h). Vital signs (SpO_{2}, ECG, blood pressure, end-tidal CO_{2}) were monitored continuously. Neuronal activity was recorded from opercular V1 (foval and parafoveal representations) using microelectrodes (300,800 kOhms; FHC). Signals were high-pass filtered (1 Hz; digital two-pole Butterworth filter), amplified using an Alpha Omega amplifier system and digitized at 20.83 kHz. The extracellular field potential was filtered to extract multiunit activity (MUA) and LFPs using standard signal conditioning techniques as described in detail previously (Belitski et al., 2008). Spikes were identified within MUA with a threshold for spike detection set at 3.5 SDs of the overall MUA signal. A spike was recognized as such only if the last spike occurred more than 1 ms earlier. For the present analysis, we did not separate single and multiunits. Binocular visual stimuli were presented at a resolution of 640 × 480 pixels (field of view: 30 × 23 degrees, 24 bit true color, 60 Hz refresh) using a fiberoptic system (Avotec). Stimuli consisted of naturalistic complex and commercially available movies (30 Hz frame rate), from which 120–210 s sequences were presented and repeated 30–50 times. For all of the videos used as stimuli, we computed from RGB values the biluminance value of every pixel (*x*,*y*) in the screen for every frame, as described by Balters et al. (2008). We determined for each electrode a 1° × 1° receptive field (RF) using reverse correlation (Ringach and Shapley, 2004) between the MUA and the time contrast (i.e., the interframe luminance difference) of each pixel in the screen. Details of this calculation have been described previously (Ince et al., 2012). For each RF and each instant of time, we computed five visual features: (1) luminance, the average luminance over all the pixels in the RF; (2) spatial contrast, the Michelson contrast (difference between maximal and minimal luminance in the RF divided by their average); (3) temporal contrast, average over all of the pixels in the RF of the absolute value of the luminance difference relative to the previous frame; (4) orientation, the direction of gradient of contrast (Kayser et al., 2003); and (5) color, the red–green color axis extracted from the RGB values.

#### Model

In the following section, we define the structure and the equations of our model of network of spiking neurons and define and name all of its relevant parameters. Actual parameter values used in the simulations will be specified when appropriate.

We consider a fully connected network of leaky integrate and fire (LIF) neurons, with *N*_{E} excitatory neurons and *N*_{I} inhibitory neurons. Each neuron, *i*, receives two sources of input current: an external current *I*_{ext,i} and a recurrent current *I*_{rec,i}, due to spikes emitted by other neurons of the network. The external current is assumed to represent all inputs from outside of the local circuit (i.e., feedforward inputs from LGN, long-range horizontal connections within V1, and feedback connections), whereas the recurrent current is assumed to represent inputs from the local network itself.

The membrane potential of a neuron *i* evolves in the subthreshold regime according to the following formula:
where τ_{m,i} is the membrane time constant and the currents are expressed in units of potential. If the voltage reaches the threshold *V*_{θ,i}, a spike is emitted and the voltage is reset to a value *V*_{r,i} during an interval τ_{rp,i}.

Spikes arriving at synapses from excitatory neurons produce AMPA-receptor-mediated currents, *I*_{A}, whereas spikes arriving from inhibitory neurons produce GABA-receptor-mediated currents, *I*_{G}. Because the network is fully connected, all neurons belonging to the same population *a* (*a* = *E*, *I*) receive the same recurrent synaptic input as follows:
The evolution of *I*_{A,a} = *s*_{A,a} and *I*_{G,a} = *s*_{G,a} can be described using the auxiliary variables *x*_{A,a} and *x*_{G,a} as follows:
where τ_{rA,a} (τ_{rG,a}) is the rise time of the synaptic AMPA (GABA) currents, τ_{dA,a} (τ_{dG,a}) is the decay time of the synaptic AMPA (GABA) currents, *t*_{j,k} is the time of the kth spike arriving at synapse *j*, and τ_{L} is the latency of the postsynaptic currents. *J*_{a,E} (*J*_{a,I}) is the individual synaptic efficacy from excitatory (inhibitory) neuron to neuron belonging to population *a*, which scales as the inverse of the number of neurons in the population. These equations give rise to excitatory/inhibitory postsynaptic currents (EPSCs/IPSCs), the time course of which is described by a delayed difference of exponentials, similar to experimentally recorded EPSCs and IPSCs (Destexhe et al., 1998).

#### External input to the network

In addition to the recurrent current, each neuron receives an external current which describes the input to the network as follows:
where μ_{ext,a}(*t*) is the average external input into population *a* at time *t*, σ_{ext,a} is the amplitude of the fluctuations around μ_{ext,a}(*t*), and η_{i}(*t*) is white noise uncorrelated from neuron to neuron (<η_{i}(*t*)> = 0 and <η_{i}(*t*)η_{j}(*t*′)> = δ_{ij}δ(*t* − *t'*)). The external input μ_{ext,a}(*t*) evolves according to an Ornstein-Uhlenbeck (OU) process as follows:
Where τ_{OU} is the characteristic time constant of the OU process, μ_{OU,a} is the mean value of the input to population *a*, σ_{OU} represents the fluctuations around μ_{OU,a}, and *z*(*t*) is a Gaussian white noise. An OU process is a stochastic process in continuous time, which is a simple generalization of a random walk. The variance of the process remains finite in the long time limit because it contains a drift term that tends to bring back the process toward its mean. It is the simplest stationary stochastic process in continuous time that has a nonzero correlation time constant and is often used in physics, biology, mathematics, and finance to model fluctuations with a power spectrum that decays monotonically with frequency. In fact, the power spectrum of such a process decays as the inverse of the square of the frequency beyond a frequency that is inversely proportional to its time constant (see Equation 41). Therefore, most power is concentrated at low frequencies. This is a very natural choice for a model of external inputs to V1 during presentation of a natural movie because time changes of visual features in natural movies have been reported, both in general (Wong and Atick, 1995) and for the particular movies used here (Montemurro et al., 2008), to have the most power at low frequencies, with the power decreasing approximately inversely proportionally to the square of the frequency. Moreover, previous work showed that the power spectrum of the multiunit firing in LGNs recorded in response to the same movies used here also has the most power at very low frequencies (Mazzoni et al., 2010) and that the low-frequency activity of model recurrent networks of excitatory and inhibitory neurons entrains to the low-frequency fluctuations of this geniculate firing (Mazzoni et al., 2008). Finally, massed cortical activity in sensory regions (measured with LFPs or MEG) entrains to the time course of slow variations in natural stimuli (Chandrasekaran et al., 2010; Gross et al., 2013). As we will see later, it turns out that, despite its simplicity, the OU process provides a good fit of experimentally recorded LFPs. Importantly, the calculations performed in the next section do not depend on the process being an OU process and can be easily generalized to arbitrary stochastic processes with a given spectrum.

#### Computation of LFPs from the model

Different approaches have been used to model LFPs. Lindén et al. (2011) used a biophysical modeling approach with detailed, spatially extended models of cortical neurons and their spatial distribution in the cortical network. Because our model does not take into account the spatial organization of cortical neurons, we used the simpler approach of Mazzoni et al. (2008) to describe the generation of the LFP. This approach is based on the fact that the major contribution to the LFP generation is given by pyramidal cells because their apical dendrites are arranged in approximate open field configurations (Protopapas et al., 1998). In this perspective, the dendrites of pyramidal cells are seen as dipole-like structures in which currents flow in the cell through apical excitatory synaptic contacts and flow out through basal inhibitory contacts (Leung, 1991). The LFPs can then be modeled as the sum on pyramidal cells of the absolute values of the AMPA and GABA currents and the external current: *LFP* = |*I*_{A,E}| + |*I*_{G,E}| + *I*_{ext,E}. This LFP proxy, though simplified, was shown to be sufficient to reproduce quantitatively several important properties of cortical field potentials, including spectral shapes and spectral information content (Mazzoni et al., 2008), cross-frequency and spike-field relationships (Mazzoni et al., 2010), and LFP phase of firing information content (Mazzoni et al., 2011).

#### Numerical integration methods

The equations for the membrane potential *V*_{i}(*t*) (Equation 1) plus conditions for spike emission and refractoriness with recurrent currents *I*_{rec,a} (Equations 3–6), external currents *I*_{ext,a} (Equations 7–8) are integrated using a Runge–Kutta method with Gaussian white noise with time step *dt* (Honeycutt, 1992). The mean value of the external currents, μ_{ext,a} (*a* = *E*, *I*), was calculated using mean-field equations to control the values of the firing rates ν_{a}, as follows:
The LFP was calculated as described in the “Model” section and, because the sample frequency of the experimental data was equal to 1 kHz, it was updated every 1 ms. The parameters used in the simulations are as follows: *N*_{E} = 1600, *N*_{I} = 400, *J*_{EE} = 0.05, *J*_{IE} = 0.085, *J*_{EI} = 0.15, *J*_{II} = 0.255, τ_{m,E} = 20 ms, τ_{m,I} = 10 ms, *V*_{θ} = 18 mV, *V*_{r} = 11 mV, τ_{rp,E} = 2 ms, τ_{rp,I} = 1 ms, τ_{rA,E} = τ_{rA,I} = 0.5 ms, τ_{dA,E} = τ_{dA,I} = 2 ms, τ_{rG,E} = τ_{rG,I} = 0.5 ms, τ_{dG,E} = τ_{dGI} = 5 ms, τ_{L,I} = τ_{L,E} = 1 ms, σ_{ext,I} = σ_{ext,E} = 5 mV, ν_{E} = 3 *Hz*, ν_{I} = 12 *Hz*. Most of these parameters are consistent with anatomical and electrophysiological data: the ratio between inhibitory and excitatory neurons (Braitenberg and Schütz, 1991; Sahara et al., 2012), membrane time constants (McCormick et al., 1985), synaptic time constants (Markram et al., 1997; Xiang et al., 1998; Zhou and Hablitz, 1998; Bartos et al., 2001), amplitude of the external noise (Destexhe and Paré, 1999; Anderson et al., 2000; Haider et al., 2013), and mean firing rates (Wilson et al., 1994; Hromádka et al., 2008). The simulations lasted *T* = 2 s and we used a time step of *dt* = 0.05 ms.

### Computation of power spectra and of their variance

The LFP power spectrum was calculated, for both simulation and data, using Welch windowing and sampling frequency of 1 kHz. With this method, the variance of spectral estimation is known to be well approximated by the expression *l*_{f} is the mean of the power spectrum at frequency *f* and *K* is the number of windows used to evaluate it (Welch, 1967). This expression of the variance of the LFP power at frequency *f* was used in the analytical calculation of the information about the stimuli carried by the LFP power (see Computation of Information about stimulus carried by the LFP power, below).

#### Computation of Information about stimulus carried by the LFP power

To characterize how well the LFP power *l*_{f} computed using a certain window of length *T* at a given frequency *f* encodes about which scene of the movie (of a set *S*) was being presented, we computed the mutual information *I*(*S*; *l*_{f}) between the power at frequency *f* and the set of scenes, as follows:
where *P*(*s*) is the probability of scene *s*, *P*(*l*_{f}) the probability of the frequency *f* to have power *l*_{f} over all trials and all stimuli, and *P*(*l*_{f}|*s*) the probability of *l*_{f} to be observed when scene *s* is presented.

In practice, information in Equation 11 can be estimated for both model and real LFP time series by using known asymptotic properties of spectral analysis, as follows. For a finite interval of time, the variance of the power is given by σ^{2} = 11*l*^{2}/9*K*, where *l* is the mean value of the LFP power and *K* is the number of windows (Welch, 1967). We assume that the distribution of the LFP power, for each scene *s*, is close to a Gaussian with mean *l*_{s} and variance σ_{s}^{2}. We can thus rewrite Equation 11 as follows:
Using a change of variables, *x* = *l*_{s} + *z*σ_{s}, and imposing σ_{s}^{2} = *l*_{s}^{2}/*F* and σ_{t}^{2} = *l*_{t}^{2}/*F* with *F* = 9*K*/11, we obtain the following:
where:
The mean value of LFP power at a given frequency and scene can be evaluated empirically from real recordings and can be computed analytically using Equation 37 when considering model responses.

Note that the accuracy of the Gaussian approximation in these data was previously reported (by comparison between Gaussian estimates and direct nonparametric estimates of mutual information) to be accurate within the range of a few percentage points (Magri et al., 2009; Magri et al., 2012). We refer the reader to these previous studies for further details. In brief, the reason of the good accuracy of the Gaussian approximation to information rests on the fact that power spectra estimates are asymptotically Gaussian and that, under the conditions investigated here, the SD of true spectra is very well described by their asymptotic behavior predicted by the Gaussian estimates (see Fig. 2*B*).

#### Procedures for fitting data to models and for quantifying goodness of fit

To fit the experimental data of the LFP power using the analytical expression (Equation 37), we used the following fitting procedure. For each scene, the LFP power of the observed data was averaged across trials to obtain for each frequency the mean value and the SE. Then, to find the best values of the model parameters given the observed data, we minimized the reduced χ^{2} as follows:
where *O*_{f} and *SE*(*O*_{f}) are the mean and the SE, respectively, of the observed LFP power at frequency *f*, and *E*_{f}(*p*_{1}… *p*_{N}) is the value of the LFP power expected from the theory at frequency *f* and with parameters (*p*_{1}… *p*_{N}), and *DF* is the number of degrees of freedom ([(number of frequencies) − (number of free parameters) − 1]).

The analytical expression of the LFP power depends on the parameters set (*p*_{1}… *p*_{N}) in a complex way. This fact implies that different parameter sets generate very similar power spectra and that the landscape of function (Equation 15) presents several local minima. Therefore, to explore the parameters' space extensively, we minimized Equation 15 using a standard simplex method routine (Press et al., 1992), starting from several different initial conditions. Moreover, as explained in section “Calculation of the LFP Power and Information,” the parameter set must be in a region in which the asynchronous state is stable. Therefore, for each parameter set that was found as outcome of a run of the minimization algorithm, we used Equation 21 to determine whether the parameter set belonged to the stability region of the asynchronous state; if not, then that parameter set was discarded. To validate our fitting procedure, we applied it to synthetic LFP data for which the ground-truth values of the network parameters was known. First, we computed the mean of the LFP power using the analytical expression derived below in the Results section (Equation 37) and its variance using the Welch formula. Second, we randomly generated 40 different LFP spectra (matching the 40 trials of the experimental data) with a Gaussian distribution with the analytically computed mean and variance. Finally, we tested the minimization algorithm on these synthetic LFP spectra. To investigate the stability of the fitting procedure, we repeatedly ran the algorithm starting from 50 different initial conditions that were randomly drawn in a region of the parameter space around the “true” values (more precisely, if *x*_{i} is the true value of parameter *i*, the initial conditions were drawn randomly in the region [0.1*x*_{i}, 1.9*x*_{i}]). At the end of this procedure, we had, for each initial condition, the values of the χ^{2} (Equation 15) together with an associated estimated parameter set. To quantify the performance of the minimization algorithm, we calculated the Euclidean distance between the true and estimated parameters.

The parameter set (*p*_{1}, …, *p*_{N}), consists of two different types of parameters. The first type are the parameters that vary when the external input changes and can therefore be different in each scene; that is, the firing rates of the excitatory and inhibitory neurons (ν_{E}, ν_{I}), the amplitude of the external noise on excitatory and inhibitory neurons (σ_{E}, σ_{I}), and the parameters that describe the input (σ_{OU}, τ_{OU}). The second type are the parameters that are intrinsic to the network and that, for a given channel, are supposed to be the same for all the scenes; that is, the synaptic efficacies (*J*_{EE}, *J*_{EI}, *J*_{IE}, *J*_{II}).

We estimated the fraction of LFP data variance explained by the model using a cross- validated procedure. For each recording channel, we split the trials randomly into two nonoverlapping sets of equal size: the “training set” and the “test set.” We then computed the best-fit model parameters by minimizing χ^{2} on the training set and quantified the fraction of variance, *R*^{2}, of the LFP power in the test sample that could be explained by the model with the parameter that best fitted the training sample. *R*^{2} is defined as the correlation between the mean values over the test trials of the observed LFP power and the corresponding values of the fit obtained using the training sample. We calculated two types of fractions of variance: the first, *R*_{s}^{2}, quantifies the fraction of variance across frequencies for a given scene and is given as follows:
Where *O*_{f} and *L*_{f} represent, respectively, the mean value over the test trials of the observed LFP power at frequency *f* and the corresponding value of the fit obtained using the training sample, whereas *Ō* and *L̄* are the respective means across frequencies. *R*_{s}^{2} gives a measure of the covariance across all frequencies, in a given scene, between the mean value of the observed LFP power and the value of the fit.

The second value, *R*_{f}^{2}, quantifies the fraction of variance across scenes for a given frequency *f* and is defined as follows:
Where *O*_{s} and *L*_{s} represent, respectively, the mean value over the test trials of the observed LFP power of the scene *s* at frequency *f* and the corresponding value of the fit obtained using the training sample of the same scene *s*; in this case, *Ō* and *L̄* are the respective means across scenes. *R*_{f}^{2} gives a measure of the covariance at a given frequency across all the scenes between the mean value of the observed LFP power and the value of the fit.

## Results

### Calculation of the LFP power and information

In this section, we derive the analytical expression of the rate and LFP power spectra of the network model, as well as the information about the stimulus carried by the LFP power spectrum.

We consider a fully connected network of *N*_{E} excitatory and *N*_{I} inhibitory LIF neurons as described in the “Model” section with a dynamical external input given by Equation 7. The calculation was performed in three steps: We first computed the stationary state of the network in the large *N* limit using a mean-field approach. We then checked the stability of this stationary state using linear stability analysis. Finally, we applied linear response theory to study the response of a finite-size network to a time-dependent input.

### Computation of the stationary state of the network in the limit of large network size

A network such as the one described in the “Model” section can be characterized using a mean-field approach (Brunel and Hakim, 1999; Brunel, 2000; Brunel and Hansel, 2006; Ledoux and Brunel, 2011). Using this approach, its stationary state(s) can be characterized analytically as a function of network parameters. In an infinitely large network, the stationary state of such a network is characterized by constant average firing rates and asynchronous firing of the neurons (hence the term “asynchronous state”). The mean firing rates of population *a* = *E*, *I* are given by the following:
Where the static transfer function of single neurons is given by the following:
This static transfer function gives the firing rate of a neuron of population *a* = *E*, *I* as a function of its mean inputs μ_{a0}, for a given level of fluctuations in the inputs σ_{a0}. In other words, it represents the *f–I* curve of neurons in the network. Note that transfer functions of LIF neurons have been used successfully to fit *f–I* curves recorded in cortical pyramidal cells in the presence of noise (Rauch et al., 2003).

### Stability analysis of the stationary state

Equations 18–20 hold if the asynchronous state is stable. However, the stability of the asynchronous state is not guaranteed for any choice of the parameters. The stability of the asynchronous state can be assessed using a standard linear stability analysis (Brunel and Hakim, 1999; Brunel and Hansel, 2006; Ledoux and Brunel, 2011). This consists of perturbing the stationary (asynchronous) state with a small exponential perturbation, ε*e*^{λ} (where ε ≪ 1 and λ = α + *i*ω). The eigenvalues λ that govern the fate of such perturbations are the solutions of the following characteristic equation:
Where *A*_{ab} = *J*_{ab}*N*_{a}*S*_{ab}, in which *N*_{a}(λ), the neuronal filter, and *S*_{ab}(λ), the synaptic filter, are given by the following:
Where
Where Γ are gamma functions and *M* are confluent hypergeometric functions (Abramowitz and Stegun, 1970), whereas:
The neuronal filter describes how the instantaneous firing rate of a neuron responds to a perturbation in its synaptic inputs, whereas the synaptic filter describes how the synaptic currents respond to a pertubation in presynaptic firing rates (Brunel et al., 2001; Fourcaud and Brunel, 2002; Brunel and Wang, 2003; Brunel and Hansel, 2006; Ledoux and Brunel, 2011).

If the real part of all eigenvalues λ (solutions to Equation 21) is negative, then fluctuations are damped and the asynchronous state remains stable. If at least one eigenvalue has a positive real part, fluctuations are amplified. In that case, there is a pair of complex eigenvalues with positive real part and the network displays synchronous oscillations with a frequency determined by the imaginary part of the corresponding eigenvalues. In this study, the parameters of the network will be always chosen in the region in which the asynchronous state remains stable to small perturbation.

### Linear response theory: response of a finite-size network to a time-dependent input

As mentioned previously, in the large *N* limit, if the network is in the asynchronous stable state and receives a constant input, the average firing rates of the excitatory and inhibitory populations, given by Equation 18, are constant and so is the LFP. A finite network receiving time-dependent inputs is subject to two distinct sources of fluctuations: an external source and an intrinsic source. The external source of fluctuations in the average activity of the network is given by the external dynamic input, which is modeled as an OU process with a slow varying time constant. The intrinsic source of fluctuations is related to the fact that the underlying network we deal with is of finite size. In a network of finite size, Equation 18 must be corrected to take into account finite-size fluctuations. As we will see in the following section, these two sources of fluctuations affect the LFP power in different frequency regions. The slow time-varying external input is reflected in the low-frequency region of the LFP power spectrum (as observed in Mazzoni et al., 2008), whereas the finite-size fluctuations are amplified in the gamma-frequency region of the LFP power spectrum. This latter amplification effect is due to a network phenomenon that has been studied extensively using linear response theory (Brunel and Hakim, 1999; Mattia and Del Giudice, 2004; Lindner et al., 2005). Finite-size effects provoke small perturbations in the global network activity, the decay of which is governed by the eigenvalues of the system (i.e., the solutions to Equation 21). In E-I networks, the eigenvalues with the largest real part have a nonzero imaginary part, leading to damped oscillations at a frequency that is typically in the gamma range (Brunel and Wang, 2003; Geisler et al., 2005; Ledoux and Brunel, 2011), and therefore to a peak in the LFP power spectrum in that frequency range.

The first unavoidable source of fluctuations around stationary firing rates are finite-size effects. Finite-size effects can be evaluated approximating the total activity of a network of size *N* by a Poisson process of rate *N*ν, where ν is the firing rate of individual neurons. When *N* is large, the resulting fluctuations around the mean rate ν are well described by a white noise term of variance ν/*N* (Brunel and Hakim, 1999). Therefore, we can write the expression of the average instantaneous firing rates of excitatory and inhibitory populations as follows:
Where ν_{a0} (*a* = *E*, *I*) is a constant term, δν_{a}(*t*) represents the fluctuations due to the dynamic external input, and the third term represents the finite-size fluctuations in which η* _{a}*(

*t*) is a Gaussian white noise.

The second source of fluctuations are due to time-dependent external inputs to the whole network. These external inputs can be described in general as follows:
Where μ_{ext,0} are the average inputs and δμ_{ext}(*t*) describes the fluctuations in inputs seen by all neurons in the network. In the following, we compute the response of the network to arbitrary δμ_{ext}(*t*). We will then specialize our calculations to the OU process described in the Materials and Methods section, Equations 7–8.

Using Equations 25–27 and assuming that the excitatory and inhibitory populations receive the same external input, we can write the total input to the excitatory and inhibitory populations in terms of synaptic currents *S*_{ab}(*t*) as follows:
Then, to compute the instantaneous firing rate, we use a Volterra series up to the first order, as follows:
Where *K*_{μa,0} (*t* − *t*′) is a kernel (also called impulse response) that describes how an input at time *t*′ affects the firing rate at time *t*.

Putting Equations 28–30 together, we obtain the following:
We now perform a Fourier transform of Equations 31 and 32. The Fourier transforms of the excitatory and inhibitory firing rates, δν_{E}(ω) and δν_{I}(ω) are given by the following:
Where *N*_{a}(ω) and of *S*_{ab}(ω) are the neuronal and synaptic filters, given by Equations 22–24, in which λ = *i*ω, and η̇_{a} is the Fourier transform of

From the above equations, we can calculate the power spectra of the average firing rate of excitatory and inhibitory neurons, which are given by the following:
Where *Y*(ω) = (1 − *A*_{EE}(ω)) (1 + *A*_{II}(ω)) + *A*_{IE}(ω)*A*_{EI}(ω).

We can now compute the power spectrum of the LFP, using the fact that the LFP is given by the sum of the absolute value of the currents on excitatory neurons. We find that the power spectrum is given by the following:
Equations 37–39 show that the LFP power is composed by two terms: *L*_{L,FS}, which depends on finite-size fluctuations and determines the shape of the power spectrum in the gamma-frequency region, and *L*_{L,ext}, which depends on the power of the external input and has a major effect on the low-frequency region of the spectrum.

As mentioned in the Materials and Methods, we chose to describe the external input with an OU process (Equation 8). If the input is governed by an OU process, then the instantaneous value of the current is given by the following:
The Fourier transform of Equation 40 is as follows:
Therefore, |δμ_{OU}(ω)|^{2} behaves as a low-pass filter with cutoff frequency equal to 1/2πτ_{OU} and a decay proportional to 1/ω^{2} at high frequency.

From Equation 37, one can derive the behavior of the LFP power at high frequency. Because |δμ_{OU}(ω)|^{2} behaves like a low-pass filter, the term in Equation 37 that explicitly depends on the external input, *L*_{L,ext}, becomes negligible at high frequency. The term *L*_{L,FS} becomes in the high frequency limit, ωτ_{m} ≫ 1, as follows:
Equation 42 shows that the LFP power at high frequency is proportional to the firing rates of excitatory and inhibitory neurons. This behavior is due to the fact that the total firing rate of both excitatory and inhibitory populations can be well approximated by Poisson processes because they are a superposition of essentially uncorrelated spike trains of individual neurons, which are themselves close to Poisson. It is known that the high-frequency limit of a power spectrum of a Poisson spike train is proportional to the total number of spikes in the window in which the spectrogram is computed. Because spikes are events narrowly localized in time, their Fourier transform has a constant tail of large power at high frequencies; therefore, the high-frequency tail of the spectrum of a sequence of spikes distributed according to a Poisson process is approximately proportional to the total number of spikes (Bair et al., 1994). Because the LFP is a sum of two components that are both a convolution of the spike times with the relevant synaptic response functions, it is expected that the high-frequency end of the LFP spectrum, too, will be proportional to the total number of spikes fired by the excitatory neurons divided by a factor ω^{4}, which is due to synaptic filtering (see Equation 24).

To test our analytical results, we performed simulations of a fully connected network of *N*_{E} excitatory and *N*_{I} inhibitory LIF neurons as described in the “Numerical Integration Methods” section and with an external input given by Equation 7. The external input (Fig. 1, top) modulates both the activity of the excitatory and inhibitory populations (see rasters and instantaneous average firing rates of excitatory and inhibitory populations in Fig. 1, middle) and LFP (Fig. 1, bottom).

In Figure 2*A*, we plotted the LFP power with different values of the amplitude of the fluctuations of the OU process describing the external inputs. As expected, the external inputs mainly affect the power in the low-frequency range. The agreement between simulation and theory is very good for small σ_{OU} and degrades progressively when σ_{OU} increases, as expected from a calculation using linear response theory.

The agreement between the theoretical curves and the simulations is good, not only for the mean LFP power, but also for the SD of the LFP power (Fig. 2*B*). The theoretical LFP power SD was calculated using the Welch formula for the variance, σ^{2} = 11*L*^{2}/9*K*, in which K is the number of windows used to calculate the simulated LFP power, whereas the SD of the simulated LFP power was calculated with both the canonical SD and the Welch formula. A further validation of our analytical approach is given by the agreement between the information carried by the theoretical and simulated LFP powers calculated using Equation 13 (Fig. 2*C*). The information in Figure 2*C* was calculated using three LFP power spectra (theoretical and simulated) with different values of the parameters that control the external input (μ_{OU} = 20 mV and σ_{OU} = 3 mV, μ_{OU} = 25 mV and σ_{OU} = 5 mV, and μ_{OU} = 30 mV and σ_{OU} = 0 mV). As expected, the information displays two well separated peaks: a peak in the very low-frequency region due to the modulation of the external input and a second peak in the gamma frequency. The peak in the gamma-frequency region is due to changes in the amplitude of the finite-size term in Equation 38 caused by modulation of the average firing rates of the excitatory and inhibitory populations by the external input. Even if the modulation of the peaks in the LFP power spectrum is due to variation of the external input, their origin is intrinsically different: the peak at gamma frequency is a network phenomenon due to amplification of finite-size fluctuations that depend on the average firing rate, which in turn depend on the amount of external input received by the network. The peak at low frequency is a direct reflection of the temporal variations of the external input.

### Fitting the model to synthetic data obtained from network models of known characteristics

We then investigated whether the fitting procedure was able to recover with a reasonable accuracy the parameter values of the model network from LFP spectra obtained through our analytical calculations. This is an important validation step, because our analytical work shows that the LFP power depends on network parameters in a complex, nonlinear way. For example, it is therefore conceivable that the presence of local minima or of degeneracies in the cost function could prevent accurate parameter estimation.

We explored this issue by applying our fitting procedures to synthetic LFP data for which the ground-truth value of the network parameters was known (see Materials and Methods for details). At the end of this procedure we had, for several initial conditions, the value of the χ^{2} (Equation 15) and the best sets of parameters. For each set of parameters, we calculated the Euclidean distance from the true values and plotted (Fig. 3*A*) the χ^{2} values as a function of the Euclidean distance. Figure 3*A* shows that the goodness of the fit critically depends on the initial condition and also that χ^{2} values with *p* < 0.05 (see insert in Fig. 3*A*) correspond to sets of parameters that can be far from the real set. In Figure 3, *B* and *C*, we show the distribution of the parameters for all the cases in which the χ^{2} had a *p* < 0.05. Even if the distributions of the parameters are quite broad, we found that the medians of the distributions are a very good approximation of the real parameter values from which the data were generated.

Our purpose will be to fit the experimental LFP power spectra for all the scenes of the movie with a fixed set of synaptic efficacies because the underlying network is not changing in time, and with the other parameters different for each scene, to capture the variability due to the input dynamics. We therefore made a second test, as the one described in the preceding paragraph, in which we fixed the synaptic efficacies at their median values. The outcome is shown in Figure 4, *A* and *B*. Even if low χ^{2} values correspond in many cases to sets of parameters distant from the real set and the distribution of the parameters are rather broad, the median of each distribution is still quite close to the real values of the parameters.

### Fitting the model to V1 LFP power spectra during stimulation with natural movies

We then investigated how well the analytical expressions (Equations 37–39) that we derived for the LFP spectra could fit LFPs recorded in monkey primary visual cortex during binocular stimulation with movies. Making use of the validation made on synthetic data, we performed the fit of the experimental data using the following procedure. First, for each channel and each scene of the movie, we ran the fitting algorithm with all the 10 parameters free to vary scene-by-scene and then calculated the median values of the synaptic efficacies of all the scenes that had χ^{2} with *p* < 0.05. Second, we repeated the fit of the LFP power spectrum of all of the scene, fixing the synaptic efficacies at the median values found in the first step. In this case, for each scene, the fit was repeated starting from 20 different initial conditions. The final values of the parameters were taken as the median of the corresponding distributions calculated using the cases in which the χ^{2} had *p* < 0.05. At the end of each minimization for both steps, a check of the stability of the asynchronous state was made, as explained in the “Procedures for Fitting Data to Models and for Quantifying Goodness of Fit.” This procedure is depicted schematically in Figure 5.

### Goodness of fit and fraction of explained variance

The procedure explained in the previous section provides excellent fits of a subset of LFP power spectra. An example of good fit is shown in Figure 6*A* (blue curve). The percentage of scenes that were well fitted using the criterion χ^{2} < χ^{2}(*p* = 0.05) varies, from channel to channel, between 7% and 86% with a distribution of mean equal to 28%, which is displayed Figure 6*B*, left. Figure 6*A* shows two other examples of fits that also capture reasonably well the real LFP power spectra, even though the condition χ^{2} < χ^{2}(*p* = 0.05) was not satisfied. The purple curve in Figure 6*A* represents a fit for which χ^{2}(*p* = 0.05) < χ^{2} < χ^{2}(*p* = 0.001). Using this criterion, χ^{2} < χ^{2}(*p* = 0.001), we obtained a broader distribution of well fitted scenes (Fig. 6*B*, middle); in this case, the percentage of well fitted scenes ranged from 11% to 95%, with a mean of 43%. The experimental LFP power spectra decreased with frequency, with either a shoulder or a broad peak in the gamma range. These features could be very well explained by the model. In some recordings, additional narrow peaks were seen at 30 Hz and harmonics (likely reflecting locking to the 30 Hz monitor refresh rate) and/or at 12 Hz (possibly reflecting either a signal artifact or entrainment to some temporal components of the movie or of it presentation system; see as an example the green curve in Fig. 6*A*). Clearly, these additional narrow peaks were not explained by the model (in which we did not make any attempt to include such kinds of sources of variation) and artifactually decreased the quality of the fit in those recordings.

We then validated our fitting procedure through cross-validation (see Materials and Methods). In brief, the data were divided in training and test sets. We then performed the fitting procedure on the training set and quantified the goodness of fit using the fraction variance of the LFP power computed from the test set that can be explained by the model that best fits the training set. We calculated the fraction of variance both across frequencies as a function of the scene, *R*_{s}^{2} (Equation 16), and across scenes as a function of the frequency, *R*_{f}^{2} (Equation 17). *R*_{s}^{2} is a measure of the covariance across all frequencies, in a given scene, between the mean value of the observed LFP power and the value of the fit, whereas *R*_{f}^{2} is a measure of the covariance at a given frequency across all the scenes between the mean value of the observed LFP power and the value of the fit. When using the criterion *R*_{s}^{2} > 0.95, the percentage of well fitted scenes strongly increased as shown by the distribution displayed in Figure 6*B*, right, the percentage of scenes with *R*_{s}^{2} > 0.95 ranged from 93% to 100% depending on the channel and had a mean of 98%.

Figure 7*A* shows the χ^{2} distribution over all scenes and all recordings (Fig. 7*A*, left), the distribution of *R*_{s}^{2} over all scenes and all recordings (Fig. 7*A*, middle), and *R*_{f}^{2} as a function of frequency, averaged over all recordings (Fig. 7*A*, right). The χ^{2} distribution is peaked at χ^{2} = 1.15, but presents a considerable amount of scenes with χ^{2} > χ^{2}(*p* = 0.001). Conversely, the distribution of *R*_{s}^{2} peaked at 0.98 and the percentage of fits with *R*_{s}^{2} < 0.95 is <15%. Figure 7*A*, right, shows that, on average, the frequency region that is better captured by the fitting procedure is in the high-frequency range (>50 Hz). Figure 7, *B*, *C*, and *D*, show the χ^{2} distribution, the *R*_{s}^{2} distribution, and *R*_{f}^{2} versus frequency of three different channels with different degrees of goodness of fit. Figure 7*B* displays a recording in which the χ^{2} distribution was almost entirely below the value χ^{2} (*p* = 0.001) (Fig. 7*B*, left), meaning that the majority of scenes was well fitted (see as an example the purple curve in Fig. 6*A*); therefore, the *R*_{s}^{2} values are all close to 1 and >0.95 (Fig. 7*B*, middle). Figure 7*C* displays an intermediate case in which both *R*_{s}^{2} and *R*_{f}^{2} are high, but part of the χ^{2} distribution falls above the value χ^{2} (*p* = 0.001); (see as an example of such a fit the blue curve in Fig. 6*A*). Last, Figure 7*D* shows one of the worst cases in which the χ^{2} distribution (Fig. 7*D*, left) is almost entirely above the value χ^{2} (*p* = 0.001); conversely, the fractions of variances explained by the model are at values that are comparable to the best fitted scenes (cf. Fig. 7*B*,*C*, middle and right, and see as an example of such a fit the green curve in Fig. 6*A*). Note the narrow troughs in the *R*_{f}^{2} versus *f* plots at a few frequencies (12, 30, and 60 Hz) that are likely due to the artifacts mentioned in this section, two paragraphs above.

### Best-fit parameters

The distributions of synaptic efficacies found performing the first step of the fitting procedure is shown in Figure 8*A*; in light gray are shown the distributions of the synaptic efficacies for all sessions and channels (37 recordings); in dark gray are the distributions of the synaptic efficacies for all of the channels in a single recording site. As expected, the whole distribution of the synaptic efficacies, which is obtained using the recordings from several multielectrodes placed in different portions of the cortex, is much broader than the distribution of synaptic efficacies obtained using the recordings from a single multielectrode. In Figure 8*B* are displayed the distributions of the scene-dependent parameters obtained from the fit of the LFP power of every scene of a single recording (dark bars) and of every scene of all the recordings (light bars) once the synaptic efficacies of each recording were set. For each scene, each single parameter was determined taking the median of all the values with χ^{2} < χ^{2}(*p* = 0.05), which was obtained by running the fitting algorithm from many different initial conditions. All parameters varied considerably from scene to scene and were in the range of biologically plausible values. The temporal evolution of the external input is reflected in the variability of the parameters that characterize the OU process, σ_{OU} and τ_{OU}. The variations of external input in turn modulate the excitatory and inhibitory firing rates, ν_{E} and ν_{I}, and the amplitudes of the external noise, σ_{ext,E} and σ_{ext,I}.

### Comparison of model predictions of spike rates with experimental MUA

As a further confirmation of the validity of the analytical approach presented here, we determined whether the parameters that we obtained from the fitting procedure were correlated with other experimental observables independently of the LFP. As first, we determined whether the MUA obtained by high-pass filtering the same extracellular signal used to extract the LFP trace was correlated with the firing rates of the excitatory and inhibitory populations that we obtained from the fits. The average MUA was calculated for all the scenes of a recording and correlated with the corresponding values of ν_{E}, ν_{I}, and the total firing rate given by *f*_{E}ν_{E} + *f*_{I}ν_{I} in which *f*_{E} = 0.8(*f*_{I} = 0.2) is the fraction of excitatory(inhibitory) neurons. We found significant correlations in approximately half of the recordings; the level of correlation and the fraction of channels in which a correlation was found varied depending on the criterion used to select or reject the quality of the fit. Because the MUA of a single channel is thought to reflect the activity of a few neurons, we also looked at the correlations with the total MUA (TMUA) of a multielectrode obtained as the average of the MUA of all its channels. The fraction of recordings in which a significant correlation was found when using the selective criterion χ^{2} < χ^{2}(*p* = 0.001) are shown in Figure 9, *A*, *B*, and *C*, top. The mean amount of correlation between firing rates and experimental MUA and TMUA is ∼0.5 for excitatory, inhibitory, and total firing rates, but the percentage of recordings in which it was found is generally higher for inhibitory firing rates (45% with MUA and 62% with TMUA) and total firing rates (54% with MUA and 51% with TMUA) than for excitatory firing rates (35% with MUA and 24% with TMUA). The level of correlations in each recordings varied considerably (from 0.2 to 0.9); some examples are given in the second and third rows of Figure 9, *A*, *B*, and *C*.

### Comparison of model parameters and parameters extracted from the movie

We then investigated whether correlations were present between the parameters found through the fitting procedure and the features of the movie: that is, the spatial and temporal contrast, the luminance, the orientation, and the color. Because the evolution of the average values of the features along the movie should modulate the firing rates of the neurons, we looked at the correlations between the firing rates extracted from the fits and the average values of the features of the movie. The highest values of correlation (∼0.4) were found between the firing rates and both the spatial and temporal contrast, as reported in Table 1. The fraction of recordings in which those correlations were found was higher for the spatial contrast (∼30%) than for the temporal contrast (∼15%), as shown in Figure 9, *A*, *B*, and *C*, top. Some examples of recordings in which those correlations were found are displayed in the bottom panels of Figure 9, *A*, *B*, and *C*.

Figure 9, *A*, *B*, and *C*, top, also show that the correlations between firing rates and the other features were present in a small subset of recordings. This is not surprising because, for this particular kind of naturalistic stimuli, the correlations between firing rates and luminance, color, and orientation have been shown to be mainly a consequence of the correlations between these features and spatial and temporal contrast (Ince et al., 2012).

We then confronted the features of the movie with the parameters that, in the model, characterize the external input: the amplitude of the fluctuations, σ_{OU}, and the time constant, τ_{OU}, of the OU process. We found correlations in the range [0.3–0.4] between σ_{OU} and the SD of the different features of the movie; the subsets of recordings with significant correlation were of variable size depending on the feature, from <5% for orientation up to ∼40% for temporal contrast (Fig. 9*D*). The correlations between the time constant of the OU process, τ_{OU}, and the time constant of the different features (calculated by fitting with an exponential the autocorrelogram of the features) was found in a very small subset of recordings (∼5%), with values ranging between 0.2 and 0.7.

A summary of all the mean amounts of correlation and the relative percentage of recordings in which the correlation was found are reported in Table 1 for the three different criterions; as a general trend, when the strictness of the criterion decreases, correlations and relative percentages increase. Table 1 also shows that a significant correlation was found between the amount of external noise to the excitatory and inhibitory neurons, σ_{ext,E} and σ_{ext,I}, and both the MUAs and the features of the movie. The levels of correlations between the σ_{ext} and MUAs were comparable to those found for the firing rates, but the fraction of channels in which a significant correlation was found was smaller. Because, in general, the amount of noise grows proportionally to the firing rate, the presence of those correlations is not surprising, and the correlations between the σs and the features could just be seen as a consequence of the relationship between noise and firing rates. We also performed a nonparametric permutation test to further assess the significance of the correlations. The null hypothesis distribution of correlation values was obtained by computing correlations between the signals taken in randomly paired trials. The *p* < 0.05 significance level of the absolute value of correlation obtained through 10,000 such reshuffles was 0.19.

### Information carried by LFP spectrum

Last, we compared the information carried by the experimental LFP power spectra and the information carried by the LFP power spectra obtained by the fit. The information content carried by the experimental data was calculated using the toolbox developed by Magri et al. (2009), whereas the information carried by the LFP power fits was calculated using the expression derived in the section “Computation of Information about Stimulus Carried by the LFP Power.” Both experimental and theoretical information were calculated for each recording using only the scenes with χ^{2} < χ^{2}(*p* = 0.001). Figure 10*A* compares the information content of the experimental LFP power across all recordings with the mean of the theoretical information calculated from the fits, whereas in Figure 10*B* are shown some examples for single recordings. The agreement between experimental and theoretical information is quite good in a broad range of frequencies. The discrepancy in the range 10–30 Hz could be due both to the presence of anomalous peaks in some of the recordings and to the fact that, in that frequency region, the SEs of the data are often larger that the values used in the theoretical calculation of the information (Welch approximation).

## Discussion

In this study, we considered the dynamics of a network of excitatory and inhibitory integrate-and-fire neurons and computed analytically its global firing rate in the presence of time-dependent inputs driving the whole network. This calculation then allowed us to compute an LFP—the sum of the average excitatory and inhibitory synaptic currents impinging on excitatory neurons of the network.

We then investigated whether the analytical formula for the LFP could be used to fit electrophysiological recordings from the primary visual cortex of anesthetized macaques upon presentation of a movie. As in Belitski et al. (2008), the movie was segmented in 100 2-s-long scenes and the power spectra of the LFP (both average and SD over repeated trials) were computed separately in each scene. All LFPs were then fitted simultaneously using the model with a single set of coupling strengths characterizing network connectivity, but different parameters characterizing the input in each scene. The fits were shown to be surprisingly good given the simplicity of the model. Furthermore, mean firing rates extracted from the fitting procedure were shown to be significantly correlated with MUA extracted from the same electrode and also with parameters characterizing the input to the network.

Note that we could have used the average value of the MUA of each scene to further constrain our fitting; we chose not to use this information to have an *a posteriori* test of the reliability of the values found through the fitting procedure. The fact that the *a posteriori* test was positive, meaning that significant correlations between fitted and experimental parameters were found, is of key importance because it shows that our model, despite its simplicity, is able to capture the relevant variables that govern the dynamics of the cortical microcircuit under study.

These results represent, to our knowledge, the first time an LFP computed analytically from a network model of spiking neurons has been shown to fit real data quantitatively. Not only does the model reproduce single LFP spectra, but it also captures well how LFP spectra are modulated by external inputs. The fitting procedure gave biophysically reasonable values for all parameters; synaptic strengths were in the range of experimentally observed values *in vitro*, though typically lower than the average reported values (Markram et al., 1997; Sjöström et al., 2001; Holmgren et al., 2003; Lefort et al., 2009). This could be due to the fact that, in active networks, synaptic efficacies tend to be reduced due to shunting effects. Synaptic efficacy estimates from *in vivo* experiments are indeed typically lower than *in vitro* measurements (Matsumura et al., 1996). In addition, in our calculations, synaptic efficacies enter the analytical formulas through their product with the number of relevant presynaptic neurons. In a sparsely connected network, these numbers would be smaller and therefore lead to larger synaptic efficacies. Mean firing rates were very similar to typical recorded average firing rates in cortex (Shafi et al., 2007; Hromádka et al., 2008; O'Connor et al., 2010). The amplitude of the external noise were again in the range of experimentally observed values (Destexhe and Paré, 1999; Anderson et al., 2000).

### Comparison with previous theoretical work

Our analytical calculations extend previous theoretical results on the dynamics of networks of integrate-and-fire neurons (Brunel and Hakim, 1999; Brunel, 2000; Mattia and Del Giudice, 2002; Mattia and Del Giudice, 2004; Lindner et al., 2005; Ledoux and Brunel, 2011; Buice and Chow, 2013). From the theoretical point of view, the main novelty of the present study consists of implementing simultaneously nonstationary external inputs and finite-size effects; computing analytically the power spectrum of firing rate and synaptic currents; and finally, comparing an LFP measure to real data.

Our study provides a principled way to “invert” LFP recordings by estimating the underlying parameters of the network of spiking neurons that generates the observed behavior. Our attempt to fit network models to field potential data bears similarities to work described previously by Moran et al. (Moran et al., 2007; Moran et al., 2009; Moran et al., 2011). These investigators used neural mass models (also called firing rate models) to compute MEG/EEG/LFP measures and used these measures to fit electrophysiological data using dynamical causal modeling. The main advance of our approach with respect to this previous work is that our equations for the network's LFP spectra arise from a model that incorporates the complete dynamics of its constituent excitatory and inhibitory spiking neurons, whereas neural mass models are described by phenomenological equations for the mean activity of whole populations of neurons. From the mathematical point of view, the main difference between the mass model LFP spectra derived by Moran and colleagues and our model lies in the neuronal transfer function: for networks of LIF neurons, the transfer function is given by Equation 22, whereas for a neural mass model, it is given by a much simpler low-pass filter (see also Ledoux and Brunel, 2011). This leads to two important qualitative differences between the two models. First, contrary to neural mass models, the shape of this transfer function depends in LIF neurons on both the mean external inputs and the level of external noise. This input dependence of the transfer function leads to the input dependence of the LFP spectrum. Second, for sufficiently weak noise, LIF transfer functions may have resonances at a frequency equal to the single neuron firing rate, which are not observed in neural mass models. These equations could easily be generalized to any spiking neuron model for which the neural transfer function can be computed, using, for example, the method of Richardson (2007).

Other investigators have attempted to account for power spectra of LFPs that qualitatively resemble experimentally recorded LFPs in V1. Kang et al., 2010 showed that an excitatory–inhibitory rate model of V1 can generate broadband peaks in the gamma band due to a resonance in the excitatory–inhibitory network. This regime is in practice very similar to the regime we studied here, but it was investigated in a much simpler rate model and no attempt was made to fit their model to the data. Battaglia and Hansel, 2011 showed that the chaotic dynamics of a multilayered cortical network model can lead to broadband gamma peaks in the LFP. This is a qualitatively different regime compared with ours. It would be interesting to investigate whether the chaotic scenario could also reproduce quantitatively experimental LFPs. However, this might represent a formidable challenge due to the lack of analytical expressions for LFPs in that regime.

From the neurophysiological point of view, the advantage offered by our analytical computation of mesoscopic observable quantities from networks of individual interacting neurons is that it offers a crucial and still missing theoretical step for the model-based interpretation of the dynamics of cortical microcircuits made of recurrently interacting neurons. An intriguing result that we found when inverting LFP recordings in terms of best-fit model parameters was that the LFP spectrum was particularly sensitive to the parameters defining the activity of inhibitory neurons. It has been long thought that the LFP does not capture much the activity of inhibitory neurons because the dendrosomatic dipoles emanated by these neurons cancel out due to geometric open-field configuration (Murakami and Okada, 2006). Our finding implies that, although the inhibitory cells do not contribute much directly to the mean extracellular potential, recurrent-network-model-based analysis approaches can infer indirectly their activity from their contribution to the recurrent microcircuit dynamics. In this respect, an important future research direction is to apply this model-based approach to analyzing stimulus designs that are thought to decouple the activity of excitatory and inhibitory neurons across stimulus conditions (Gieselmann and Thiele, 2008; Ray and Maunsell, 2011) to understand how different stimulus paradigms manipulate the effective ratio between excitation and inhibition.

### Implications for cortical dynamics

One of the implications of our work is that the spectral power of gamma-band activity in visual cortex can be reproduced reasonably well by a recurrent network model with dynamics in an asynchronous state. This suggests that visually evoked gamma-band activity is due to a resonant phenomenon; the model is close to a Hopf bifurcation leading to oscillatory population activity in the gamma range. These results are consistent with work from Shapley's group (Burns et al., 2010; Kang et al., 2010; Burns et al., 2011; Xing et al., 2012). Note, however, that the vicinity of a Hopf bifurcation means that the system might easily switch to a stronger oscillatory state provided the network state changes (i.e., from anesthetized to awake) and/or different external inputs are used. In particular, power spectra from awake monkeys in response to near optimal stimuli sometimes exhibit much stronger peaks at gamma frequency, consistent with a truly oscillatory regime (Gieselmann and Thiele, 2008).

We found that, in our model, the LFP power in the high-frequency limit is proportional to firing rates, as shown by Equation 42. Interestingly, in primary visual cortex, gamma-band oscillations (40–80 Hz) may under some stimulation conditions decorrelate from spike rate variations, whereas high-gamma activity (80 Hz or more) always tightly follows spiking activity (Gieselmann and Thiele, 2008; Ray and Maunsell, 2011). By careful time-frequency analysis of the extracellular signal, Ray and Maunsell (2011) showed that part of the reason that high-gamma power and spiking activity are always tightly linked is the spectral leakage of spike shapes into power in the high-gamma range. Our results, however, are obtained using model LFPs that are made only of synaptic currents and do not contain any contamination of spikes. This suggests that part of the observed link of high-gamma and spiking activity may arise from fundamental constraints on network dynamics, rather than simply on spectral leakage of spike shapes.

## Footnotes

This work was supported by the Compagnia di San Paolo, of the Max Planck Society, and of the SI-CODE project of the Future and Emerging Technologies (FET) program within the Seventh Framework Programme for Research of the European Commission (under FET-Open Grant FP7284553).

The authors declare no competing financial interests.

- Correspondence should be addressed to Francesca Barbieri, Unit of Neuroscience Information and Complexity, Centre National de la Recherche Scientifique Unité Propre de Recherche-3293, 1 Avenue de la Terrasse, Gif-sur-Yvette (Paris) 91198, France. barbieri{at}unic.cnrs-gif.fr