## Abstract

High firing irregularity is a hallmark of cortical neurons in vivo, and modeling studies suggest a balance of excitation and inhibition is necessary to explain this high irregularity. Such a balance must be generated, at least partly, from local interconnected networks of excitatory and inhibitory neurons, but the details of the local network structure are largely unknown. The dynamics of the neural activity depends on the local network structure; this in turn suggests the possibility of estimating network structure from the dynamics of the firing statistics. Here we report a new method to estimate properties of the local cortical network from the instantaneous firing rate and irregularity (CV_{2}) under the assumption that recorded neurons are a part of a randomly connected sparse network. The firing irregularity, measured in monkey motor cortex, exhibits two features; many neurons show relatively stable firing irregularity in time and across different task conditions; the time-averaged CV_{2} is widely distributed from quasi-regular to irregular (CV_{2} = 0.3–1.0). For each recorded neuron, we estimate the three parameters of a local network [balance of local excitation-inhibition, number of recurrent connections per neuron, and excitatory postsynaptic potential (EPSP) size] that best describe the dynamics of the measured firing rates and irregularities. Our analysis shows that optimal parameter sets form a two-dimensional manifold in the three-dimensional parameter space that is confined for most of the neurons to the inhibition-dominated region. High irregularity neurons tend to be more strongly connected to the local network, either in terms of larger EPSP and inhibitory PSP size or larger number of recurrent connections, compared with the low irregularity neurons, for a given excitatory/inhibitory balance. Incorporating either synaptic short-term depression or conductance-based synapses leads many low CV_{2} neurons to move to the excitation-dominated region as well as to an increase of EPSP size.

## INTRODUCTION

Firing properties of single neurons are determined both by intrinsic membrane characteristics and by their connectivity matrix within their local network. Whereas the electrical properties of a neuronal membrane are rather well understood (Hodgkin and Huxley 1952; Koch 1999), it is currently impossible to obtain the full connectivity matrix of a local network in the cortex. Even if the connectivity matrix of a specific area of a neural tissue could be obtained experimentally, it would still be desirable to represent it in the form of a set of probabilistic or descriptive rules to provide a clear insight of the network architecture. In any case, it is of interest to devise a method that provides information about the relationship between the parameters defining the connectivity and the recorded activity of neurons.

Cortical neurons have been described to have a firing irregularity close to a Poisson process (Compte et al. 2003; Shinomoto et al. 2003; Softky and Koch 1993; but see Maimon and Assad 2009; Shinomoto et al. 2009). Several theoretical studies have argued that a balance of excitation and inhibition can explain this high irregularity (Amit and Brunel 1997; Brunel 2000; Burkitt 2001; Gerstein and Mandelbrot 1964; Shadlen and Newsome 1998; Tsodyks and Sejnowski 1995; van Vreeswijk and Sompolinsky 1996). More recently, direct measurements of conductance changes during UP states in vitro and in vivo have revealed proportional changes of excitatory and inhibitory conductances (Haider et al. 2006; Shu et al. 2003). These lines of evidence support the balanced network hypothesis in which excitatory and inhibitory inputs to neurons are balanced. This balance should be, at least in part, generated by local connectivity of the cerebral cortex, but the properties of this connectivity are not fully understood yet.

Some of the previously mentioned theoretical papers derived the statistics of firing (firing rate, firing irregularities) of single neurons in a network as a function of the parameters defining network connectivity (Brunel 2000; Burkitt 2001). In this paper, we propose to investigate the reverse relationship: estimate parameters of the local circuit of a given neuron from its statistics of firing as obtained from extracellular electrophysiological recordings. The problem of finding local network parameters from the firing statistics is an inverse problem, and it is difficult to solve it in general. From theoretical and simulation studies of leaky integrate-and-fire (LIF) neurons, it is known that the balance between excitation and inhibition is a crucial parameter that determines the behavior of firing rate and firing irregularity (see e.g., Brunel 2000; Shadlen and Koch 1998). If a network is purely driven by excitatory inputs, the firing irregularity typically decreases as the firing rate increases (Fig. 1*B*). However, if the recurrent inhibition is balanced with or even dominates excitation, firing irregularity tends to be stable and decoupled from the firing rate changes (Fig. 1*C*). Such a weak modulation of the firing irregularity is also observed in vivo, in area MT and V1 (Shimokawa and Shinomoto 2009). Motivated by these findings, we set out to statistically estimate the parameter regions that can best explain the observed dynamics of the firing rate and irregularity. The method we propose here relies on the instantaneous measure of the firing rate ν and an instantaneous measure of the firing irregularity, such as the CV_{2} metric (Holt et al. 1996) as a function of time. This method is composed of the following three steps. *1*) Recording and analysis: the instantaneous firing rate and the instantaneous firing irregularity are measured within a sliding window of duration Δ*T* ms. The firing irregularity is measured by the metric CV_{2} and compared with the theory. *2*) Theory: given a particular single neuron model and a model of its local network, we predict how the output firing rate ν and the CV_{2} depend on the external input, either numerically or analytically. Changes in external inputs lead to changes in single neuron firing statistics (ν, CV_{2}) and draw a trajectory in ν − CV_{2} plane. This trajectory depends strongly on network parameters, such as the ones that control the balance between excitation and inhibition (Fig. 1, *B* and *C*). *3*) Fitting the model to the data: we calculate the distance (cost function) between the experimental data points obtained from the electrophysiological recordings and the model trajectory in the ν − CV_{2} plane. The best fit parameters are searched. We numerically searched the global minima in a discretized parameter space.

We applied this strategy to the spike data of single neurons recorded from the motor cortex of two behaving Rhesus monkeys as described in Roux et al. (2003, 2006). Analysis of these recordings reveals that firing irregularity is relatively weakly modulated in time compared with the firing rate, and its time average over the task period is nearly independent from task conditions, such as the direction of the arm movement. The distribution of CV_{2}s over the population ranges from 0.3 to >1 and is possibly bimodal, which is suggested by previous results from other cortical areas (Shinomoto et al. 2003, 2009). We then attempted to extract some information about parameters of the local circuits of the primary motor cortex. Our parameter estimation method suggests the following. *1*) Irregular spiking neurons operate in balanced or inhibition-dominated network. *2*) Quasi-regular spiking neurons can operate in an excitation-dominated region if the number of recurrent connections per neurons is very small (∼10), but as the recurrent connectivity increases, the ratio of excitation-inhibition shifts to an inhibition-dominated region. *3*) The estimated excitatory postsynaptic potential (EPSP) sizes and the ratio or inhibitory PSP (IPSP)/EPSP of irregular spiking neurons are larger than those of quasi-regular neurons if the number of presynaptic neurons is the same. As the number of local connections per neuron increases, the estimated EPSP size becomes smaller. *4*) Incorporating either synaptic short-term depression or conductance-based synapses leads to a large number of neurons (especially low irregularity neurons) to move to the excitation-dominated region and to an increase in EPSP amplitudes. Preliminary results were presented in an abstract (Hamaguchi et al. 2008).

## METHODS

### Electrophysiology and tasks

Data for this analysis were obtained from two earlier studies (Roux et al. 2003, 2006). Two male Rhesus monkeys participated in this study (data from *monkey R* are partly described in Roux et al. 2003), and single neuron spiking data from *monkey O* were recorded on the same electrodes with the local field potentials (LFPs) as described in Roux et al. (2006). Detailed descriptions of the experimental techniques were described in Roux et al. (2003, 2006). Care and treatment of the animals during all stages of the experiments conformed to the European and French government regulations.

##### EXPERIMENTAL TASKS.

The two monkeys were trained to perform two types of tasks: a choice-reaction time task (chRT) and a self-paced movement task (Self). In front of the monkey, a vertical panel with three touch-sensitive light-emitting diodes was set. The diodes were horizontally aligned. First, the center diode was lit in yellow. The monkey had to touch it to initiate the trial. Then the following two tasks were randomly performed in blocks:

##### CHOICE-REACTION TASK (chRT).

After a fixed delay (500 ms), the two peripheral diodes were lit simultaneously as the preparatory signal (PS), one in red and the other in green. The position of the color was randomly chosen. An auditory, directionally noninformative signal followed after either a short or a long delay as the response signal (RS). If the RS occurred after a short delay, the monkey had to touch the red diode, and if it occurred after a long delay, the green one should be touched. (*n* = 89 neurons; *monkey O*: *n* = 58, *monkey R*: *n* = 31).

##### SELF-PACED MOVEMENT TASK (SELF).

After a fixed delay (500 ms), only one of the two peripheral diodes was lit in red or green as PS. The position and color were randomly chosen. In this task, there was no response signal, and the monkey had to initiate the movements by himself. If the color was red, the monkey had to initiate the movement after a correctly estimated short delay. For the green diode, the monkey had to initiate the movement after a long delay. (*n* = 75 neurons; *monkey O*: *n* = 42, *monkey R*: *n* = 33).

In both chRT and Self tasks, same delay conditions were used, but the values of the delays were different for the two monkeys. In *monkey O*, two combinations of delays were used: 600–1,200 and 1,000−1,400 ms (Roux et al. 2006) and in *monkey R* only 600–1,200 ms (Roux et al. 2003), for the short and long delays, respectively. A trial was considered as correct when the movement was performed to the correct target and initiated within a temporal window of 300 ms starting at the end of the delay. Otherwise, that trial was classified as an error trial and was not rewarded. We analyzed only successful trials in this paper. For each task, movement onset (MVT) was defined when the hand releases the center yellow diode, and we aligned the time axis of each trial to the MVT.

To record extracellularly multiple single-neuron activities, a multielectrode, computer-controlled microdrive (*monkey R*: Reitboeck system, Thomas Recording, Giessen, Germany; *monkey O*: MT-EPS, AlphaOmega, Nazareth, Israel) was used to transdurally insert up to seven microelectrodes (*monkey R*: quartz insulated platinum-tungsten electrodes, 80 μm OD, impedance:: 2–5 MΩ at 1,000 Hz; *monkey O*: epoxy insulated tungsten electrodes, Frederick Haer, 0.5–1.2 MΩ at 1,000 Hz). The inter-electrode distance was between 230 and 700 μm; however, because the electrodes were driven independently from each other, their distance in depth may have varied for each recording. From each electrode, electrical signals were amplified with a gain of 5,000 to 10,000 (MCP+, AlphaOmega) and high-pass filtered by using in-house hardware (active filtering) at 300 Hz for selecting action potentials (spikes). Spikes of only one single neuron per electrode were then isolated by using a window discriminator in *monkey R* or between one and three neurons per electrode (although rarely >2 were discriminated in practice) by using on-line spike sorting with a template matching algorithm (MSD, AlphaOmega) in *monkey O*. Sorted spike data along with behavioral events (occurrences of signals and performance of the animal) were stored for off-line analysis with a time resolution of 1 kHz.

### Irregularity measures

Several measures of spike time irregularity have been proposed so far in different contexts. The most widely used measure, the coefficient of variation (CV) of interspike intervals (Softky and Koch 1993), is defined as
*x*〉 is the sample average of *x*. The CV is equal to 1 for a stationary Poisson spike train and takes values of 0 for a completely regular spike train. If the spike train is not stationary, i.e., the firing rate changes in time, the CV tends to overestimate the irregularity of the underlying process (Nawrot et al. 2008). For example, if the firing rate of a Poisson spike train changes in time, the CV computed from finite time windows becomes >1 around the time of the firing rate changes (Fig. 1*A*). Several alternative measures of firing irregularity that are more robust to rate modulations have been proposed recently: CV_{2} (Holt et al. 1996), LV (Shinomoto et al. 2005), IR (Davies et al. 2006), SI (Miura et al. 2006), and CV in operational time (Nawrot et al. 2008). CV_{2}, LV, IR, and SI are calculated based on neighboring ISIs, and they become zero for completely regular spike trains. According to their original definition, only CV_{2} and LV take values of 1 for a Poisson spike train, but IR and SI can also be normalized such as to be unity for the Poisson case.

In this paper, we have chosen the CV_{2} (Holt et al. 1996) as a local measure of spike time irregularity because it provides a less biased estimate of irregularity for a small, finite sample size (Ponce-Alvarez et al. 2010). For a more detailed comparison of these measures and their performance, see Ponce-Alvarez et al. (2010). The CV_{2} is computed in the following equations
*t _{k}* is the

*k*th spike from a cell and CV

_{2}

^{(k)}is the value of the CV

_{2}assigned to the

*k*th spike.

### Data analysis

We calculated instantaneous firing rate, ν (*t*), the instantaneous firing irregularity CV_{2}(*t*) (Holt et al. 1996), and their SEs [σ* _{n}*(

*t*), σ

_{CV2}(

*t*)] as a function of time, using a sliding window. The instantaneous firing rate at time

*t*, ν(

*t*), is calculated by taking the number of spikes divided by the length of the sliding window around time

*t*, including all trials. The SE of the firing rate σ

*(*

_{n}*t*) is the intertrial SD of the firing rate divided by the square root of the number of trials.

The instantaneous firing irregularity at time *t* is calculated as follows: a CV_{2}^{(k)} value calculated from the preceding and the following ISIs is assigned to *k*th spike in the sliding window (Fig. 2*A*), using *Eq. 3*. The first and last spikes of each trial are not taken into account in the analysis because they do not have preceding or following ISIs. The instantaneous firing irregularity CV_{2}(*t*) is the mean of CV_{2}^{(k)} values within a sliding window averaged over all the trials (*Eq. 4*). The SE of CV_{2}(*t*), σ_{CV2}(*t*), is defined as the SD of CV_{2}^{(k)} values within the sliding window over the trials, divided by square root of the number of spikes. Here we approximated the distribution of the sample mean of CV_{2}^{(k)} as a Gaussian distribution.

Whereas the use of the CV_{2} measure provides a much more robust way of assessing single neuron irregularity than the CV measure, we have noticed two issues on the use of CV_{2}. One is the validity of the Gaussian approximation of the sample mean of CV_{2} that determines the confidence limit of CV_{2} in our paper. The other is the bias toward higher values during epochs of transient very low-firing rate; even though CV_{2} is less susceptible to the changes of firing rate compared with CV, there is still a bias in CV_{2} estimates when the firing rate changes abruptly especially when firing rate drops off to near zero rate. We now address these two issues.

##### CONFIDENCE LIMIT OF CV_{2}.

The first question is how to determine the confidence limits of the measured CV_{2}(*t*) values. From the central limit theorem, the distribution of the sample mean from independently and identically distributed variables converges to Gaussian. Once a sample mean is obtained, the true mean is expected to be distributed according to a *t*-distribution around the sample mean. The *t*-distribution also converges to a Gaussian for large sample sizes. For CV_{2} measures, however, two neighboring CV_{2} values are correlated due to the shared ISI, and the average of the correlated variables does not generally converge to its true mean even for large sample sizes. To remove the effect of within sample correlation, one can average every other CV_{2} value, but this down-sampling may increase the estimation error of the true CV_{2} value. To study the properties of CV_{2} statistics, we numerically computed the distribution of sample means of CV_{2} from a gamma distributed ISI and confirmed that averaging all CV_{2} values actually yields to lower SE than taking every other CV_{2} value (Supplementary Fig. S1*B*),^{1} and that the assumption of the Gaussian distribution of the sample mean of CV_{2} is valid if there are large numbers of samples (Supplementary Fig. S2). The details of how to obtain confidence limits of the CV_{2} statistics are given in the supplementary materials. We set the threshold to 20 spike counts per sliding window above which we can safely approximate the sample mean distribution as Gaussian. We define a data point as the couple of {ν(t), CV_{2}(*t*)} calculated within one sliding window around time *t*. The data points containing <20 spike counts are not used in the parameter estimation part, and are not presented in the figures.

##### SUPPRESSING THE BIAS TOWARD HIGH CV_{2}.

The second question is how to suppress the upward bias of estimated CV_{2}s observed at abrupt firing rate changes. We find that the spikes that are nearest neighbors in a time window with low spike counts tend to have high CV_{2} values when this time window follows and/or precedes windows with much larger spike counts. This is caused by the abnormally long ISIs in this low-firing rate region. One way to correct this overestimation error is to remove spikes at the border of low to high spike count (see also Davies et al. 2006). This bias-suppression method is very effective and drastically reduces the overestimation error (Supplementary Fig. S3). The details of the bias-suppression method are described in the supplementary materials. The threshold of low to high spike count was set to 20 spike counts within the sliding window over trials. After applying this bias-suppression method, the firing irregularity at time *t*, CV_{2}(*t*), is calculated as the average of CV_{2}^{(k)} within the sliding window, including all trials.

To summarize, the firing irregularity is calculated as follows: after applying the bias-suppression method described in the preceding text, the CV_{2}(*t*) and σ_{cv2}(*t*) are calculated as the sample mean and conventional SE from the CV_{2}^{(k)} values within the sliding window over all the trials. The removal of spikes in the bias-suppression method is only applied to the CV_{2} statistics, and the firing rates are not affected.

In this paper, we used Δ*T* = 100 ms time window, and data were always aligned to MVT. Only spikes occurring during periods from MVT −1,100 to MVT +300 ms, for short delay trials and MVT −1,500 to MVT +300 ms, for long delay trials were taken into account. There are four conditions (short or long delays and left or right arm movements) for each of the two types of tasks (chRT and Self). Therefore the activity of one neuron can be recorded in eight different conditions at most. We do not focus on the description of the neural representation of the movement or the stimuli but only on an estimation of the parameters of the local network in which each neuron is embedded. We treat each condition independently and estimate the network parameters for each condition separately. In the parameter estimation part, we analytically calculate the CV_{2} of the LIF neuron model driven by a white noise (Brunel 2000; Siegert 1951; Tuckwell 1988) by using numerical inverse Laplace transform of the analytically obtained ISI distribution in the frequency domain.

### Neuron and network model

We consider a sparsely and randomly connected network of LIF neurons (Amit and Brunel 1997; Brunel 2000) (Fig. 3) as a model of a local cortical circuit in the primary motor cortex that is activated during the arm-reaching task. The membrane potential dynamics of neuron *i* (*i* = {1, 2, …, *N*}) is described as follows
*C* is the membrane capacitance, *g*_{L} is the leak conductance, *V*_{L} is the leak voltage, and *D* is the synaptic conduction delay. When *V _{i}* >

*V*

_{L}, the

*i*th neuron emits a spike and

*V*is reset to

_{i}*V*

_{r}and fixed during the absolute refractory period τ

_{ref}. The synaptic efficacy,

*J*is defined as

_{ij}*J*=

_{ij}*J*if the synapse is excitatory, or

*J*= −

_{ij}*gJ*if the synapse is inhibitory. Each neuron also receives excitatory inputs from external populations. We referred to a set of neurons in other areas of the brain that project their axons to the local network as the external population. The synaptic efficacy is

*J*=

_{ij}*J*

_{ext}where the

*j*th neuron is now in an external population. For simplicity, we set

*J*

_{ext}=

*J*in this paper. In this model, the synaptic current is described as an instantaneous current (δ-function).

The whole network is composed of *N*_{E} = 0.8*N* excitatory neurons and *N*_{I} = 0.2*N* inhibitory neurons, making the ratio of the excitatory and inhibitory neuron, β = *N*_{I}/*N*_{E} = 0.25. We also set the number of recurrent, excitatory (inhibitory) presynaptic neurons per postsynaptic neuron as *C*_{E} (*C*_{I} = β*C*_{E}).

### Self-consistent analysis of firing rate and CV

For simplicity, we assume homogeneous neuron properties for both excitatory and inhibitory neurons. Because the connectivity is the same for both neuron types, the statistical properties of the inputs are the same. As a result, both neuron types have the same firing rate. Assuming that the correlations in the inputs to different neurons in the network can be neglected, and the jump of EPSP(IPSP) is much smaller than the distance to the threshold, then one can rewrite the input of *Eq. 6* by using a diffusion approximation (Amit and Brunel 1997; Brunel 2000)
_{m} = *C/g*_{L} is the membrane time constant, 〈d*W _{i}*(

*t*)〉 = 0 and 〈d

*W*(

_{i}*t*)d

*W*(

_{j}*t′*)〉 = δ

*δ(*

_{ij}*t*−

*t′*). The mean input μ(

*t*) is given as

_{l}(

*t*) is the average input from local connections and μ

_{ext}(

*t*) is the average input from outside of the network. The firing rate ν(

*t*) is the measured spike emission rate, representing the firing rate of the local cortical circuits. We will use measured firing rate as ν(

*t*). The number of afferent connections from the external population to the network is

*C*

_{ext}. The firing rate of the external population is ν

_{ext.}The variance of the input σ

^{2}(

*t*) is

*t*). Note that ν

_{ext}enter the equations only as a product of

*C*

_{ext}and ν

_{ext}. Therefore there is only a single variable describing the total external input rate—the sum of all firing rates of the external neurons projecting to the neurons the spikes of which we are measuring. Therefore from now on, we simply describe

*C*

_{ext}ν

_{ext}as ν′

_{ext}.

To calculate the firing rate and CV_{2} of the local network, we assume that within the sliding window, a network state is very close to the stationary state (quasi-stationarity). This assumption becomes valid if the time scale of the firing rate change is much longer than the membrane time constant. The quasi-stationary firing rate and the CV_{2} value can be calculated by the following procedure: solving the self-consistent equation of the firing rate (Amit and Brunel 1997; Brunel 2000), numerically calculating the ISI distribution by using the inverse Laplace transform of the analytical expression of the ISI distribution in the frequency domain (Sugiyama et al. 1970), and computing the mean of the CV_{2} values.

The self-consistent equation of the firing rate can be obtained by solving the following equation
*f* is the transfer function indicating the output firing rate of a LIF neuron as a function of ν, ν′_{ext}, and Π, and given as
*y*_{th} = (*V*_{th} − *E _{syn})*/σ and

*y*

_{r}= (

*V*

_{r}−

*E*

_{syn})/σ,

*E*

_{syn}=

*V*

_{L}+μ, τ = τ

*. Note that μ and σ given in*

_{m}*Eqs. 8*and

*11*are the functions of firing rates ν, ν′

_{ext}and network parameters Π. The error function, erf(

*x*) is defined as erf(

*x*) = 2/

*Eq. 14*is determined by external population rate V′

_{ext}and network parameters Π, we denote the solution of

*Eq. 14*as

_{2}, we first calculate the ISI distribution from the numerical inversion of the analytical expression of the ISI distribution in the frequency domain (Sugiyama et al. 1970). The analytical expression of the ISI distribution in the frequency domain is

*p*

_{ISI}(ω|

*y*

_{th},

*y*

_{r}) is the Laplace transform of the ISI distribution for a given stationary Gaussian noise input of mean μ and SD σ.

*U*(ω,

*y*) is the combination of confluent hypergeometric functions (Brunel 2000)

*M*(λ,

*a, b*) is the confluent hypergeometric function of the first kind (see Abramovicz and Stegun 1970). Note that the ISI distribution

*p*

_{ISI}(

*t | y*

_{th},

*y*

_{r}) depends only on mean input μ and SD σ, but μ and σ are actually functions of the external firing rate ν′

_{ext}, the firing rate ν*, and the network parameters Π. Therefore we first solve the

*Eq. 16*to obtain ν* for computing μ and σ, and then we compute

*p*

_{ISI}(

*t | y*

_{th},

*y*

_{r}). From this, we calculate CV

_{2}from the numerical integration of the following equation

_{2}is a function of the network parameters and the external input, we denote CV*

_{2}= ψ(ν′

_{ext}, Π).

Now the single neuron parameters describing single neurons are τ* _{m}*, τ

_{ref},

*D, V*

_{L},

*V*

_{th},

*V*

_{r}, while network parameters are

*g, J, C*

_{E}, and ν′

_{ext}(=

*C*

_{ext}ν

_{ext}). Single neuron parameters are taken from La Camera et al. (2006) and Rauch et al. (2003) where the parameters of a single LIF neuron were fitted from the response of pyramidal cells to a noisy current injection. Here we set τ

*= 30 ms, τ*

_{m}_{ref}= 2 ms,

*V*

_{th}= −50 mV,

*V*

_{reset}= −55 mV, and V

_{L}= −60 mV. Now the remaining free parameters are

*g, J, C*

_{E}, and ν′

_{ext}. From now on, we denote a set of network parameters to be estimated as Π = {

*J, g, C*

_{E}}. The external firing rate ν

_{ext}is also one of the parameters controlling the activity, but this parameter can dynamically change during the task. We therefore left ν′

_{ext}as a dynamical variable and will not estimate it as a fixed parameter. In fact, we can marginalize the likelihood of Π by the probability

*p*(ν′

_{ext}).

### Theory for estimating network parameters from single neuron statistics

The next step is to search for the network parameters that best describe the data. For a given data set (ν(*t*_{i}), CV_{2}(*t*_{i})), we want to maximize a posterior probability *P* (Π|{ν, CV_{2}}) where Π = {*J, C*_{E}, *g*} is a set of the network parameters. From the Bayes rule, the posterior is
*p* (Π) = *c* (constant) within the searched parameter ranges to avoid making bias in the estimation of the network parameters. By choosing the sampling time *t _{i}* to be integer multiples of the sliding window size Δ

*T*= 100 ms, we can neglect the correlation within the sample mean of rate {ν(

*t*)} and within the sample mean of firing irregularity {CV

_{i}_{2}(

*t*)} where

_{i}*t*=

_{i}*i*Δ

*T*and

*i =*{1, 2, … }. Then we can decouple the probability

*p*({ν, CV

_{2}}|Π) into the product of

*p*[ν(

*t*), CV

_{i}_{2}(

*t*)|Π]. Next we assume that the estimation errors made on ν(

_{i}*t*) and CV

_{i}_{2}(

*t*) are approximately independent. As a consequence, we can decouple

_{i}*p*[ν(

*t*), CV

_{i}_{2}(

*t*)|Π] into

_{i}*p*[ν(

*t*)|Π]

_{i}*p*[CV

_{2}(

*t*)|Π]. Now the probability distribution

_{i}*p*({ν, CV

_{2}}|Π) can be decomposed into the products of the independent distributions of

*p*[ν(

*t*)|Π] and

_{i}*p*[CV

_{2}(

*t*)|Π]. Then, the posterior distribution is

_{i}_{ν}(

*t*) and σ

_{i}_{CV2}(

*t*) are SE of the rate and CV

_{i}_{2}defined at time

*t*. By using the saddle-point approximation (Daniels 1980) and taking log of posterior, we can approximate the log-posterior as

_{i}_{ext}gives the minimum of

*Eq. 23*. Here we assumed an uninformative prior for ν

_{ext}. Intuitively, the saddle-point approximation is equivalent to search for the closest point between [ν(

*t*), CV(

_{i}*t*)] and a point on a theoretical curve [ϕ (ν

_{i}_{ext}, Π), ψ(ν

_{ext}, Π); Fig. 2

*B*]. The distances are normalized by σ

*(*

_{n}*t*) and σ

_{i}_{CV2}(

*t*) for each data point. Therefore the cost for the parameter Π is

_{i}*C*

_{E}= [10, 100, 1,000],

*g*= [0, 0.25, 0.5, …, 8], and

*J*= [0.01, 0.02, …, 0.39, 0.40, 0.50, 0.60, …, 0.90, 1.00]. Because we have used a limited parameter range, this can be viewed as a prior distribution of the parameter space. For example, EPSPs size <0.01 mV are not considered here because such values are far smaller than the noise level in intracellular recordings. The performance of the estimation depends on the number of data points and their distribution over the ν – CV

_{2}plane. The larger the number of data points {ν(

*t*), CV(

_{i}*t*)} and the wider their distribution in the rate − CV

_{i}_{2}space, the better we can estimate the network parameters. Therefore we choose neurons and task conditions which have ≥10 data points, and the difference between the minimum and maximum of the firing rate is ≥10 Hz. The number of trials is >30 [

*n*= 507 conditions (

*monkey O*:

*n*= 303,

*monkey R*:

*n*= 204) from

*n*= 164 neurons (

*monkey O*:

*n*= 100,

*monkey R*:

*n*= 64)].

### Synaptic properties

To investigate how the choice of the model, such as synaptic properties, affects the estimate of network parameters, we also used other models with synapses with short-term depression (STD) and conductance-based synapses and separately estimated the network parameters in these two cases.

##### SHORT-TERM DEPRESSION.

STD is a major feature of synapses between cortical pyramidal cells (e.g., Varela et al. 1997). A simple phenomenological model has been developed by Tsodyks and Markram (1997) that provided a very good fit of their data. In the STD synapse model, the EPSP size *J _{ij}* in

*Eq. 6*is replaced by a dynamical variable that obeys

*x*(

_{j}*t*) is the fraction of vesicles available at time

*t*, and

*u x*(

_{j}*t*) is the fraction of vesicles released by a single action potential.

*x*(

_{j}*t*) decreases by a factor

*u x*at each spike and recovers exponentially with the time constant τ

_{j}_{STD}. Only excitatory synapses (both external and recurrent) are modeled according to

*Eqs. 25*and

*26*, while inhibitory synapses are kept fixed. We used

*u*= 0.5, τ

_{STD}= 500 ms. The STD level of afferent synapses now depends on the external input rate ν

_{ext}, thus the number of afferent synapses

*C*

_{ext}needs to be explicitly modeled here. Here we assumed that each neuron receives same number of afferent and recurrent excitatory synapses (

*C*

_{ext}=

*C*

_{E,}) (see e.g., Braitenberg and Schüz 1998). To analyze the network steady states in the presence of STD, we use the results of Romani et al. (2006). We first calculate the moments of

*x*in a stationary state by using Laplace transform of the ISIs

*p*

_{ISI}(

*w*|

*y*

_{th},

*y*

_{t})

*Eq. 15*to calculate the stationary firing rate

*x*

_{ext}〉 is the average fraction of vesicles available for a Poisson firing neuron at ν

_{ext}Hz and is obtained by using

*p*

_{ISI}(

*n*/τ

*) = ν*

_{STD}_{ext}/(ν

_{ext}+

*n*/τ

_{STD}) in

*Eqs. 27*and

*28*.

##### CONDUCTANCE-BASED SYNAPSES.

We have also investigated a network with conductance-based synapses. In the conductance-based synapse model, the PSP size *J _{ij}* in

*Eq. 6*is replaced by

*a*

_{E/I}is a nondimensional parameter (ratio between synaptic conductance and leak conductance) that control the size of excitatory (

*a*

_{E}) or inhibitory (

*a*

_{I}) postsynaptic potentials, and

*V*

_{E/I}are the excitatory (

*V*

_{E}) or inhibitory (

*V*

_{I}) reversal potentials. Here we used

*V*

_{E}= 0 mV,

*V*

_{I}= −70 mV. We calculate the stationary firing rate of the LIF network with the conductance based synapses by using the effective time constant approximation (Meffin et al. 2004). The following parameters are used in

*Eq. 15*

## RESULTS

Here we propose a novel method to estimate parameters of a local network in which a neuron is embedded, from the combined dynamics of firing rate and firing irregularity. Our method is valid if both the rate and the firing irregularity can be regarded as quasi-stationary within the sliding window for each sampling time *t _{i}*. We apply this method to the single-neuron spike data recorded in the motor cortex of two behaving monkeys (

*monkey O*,

*n*= 100 neurons;

*monkey R*,

*n*= 64 neurons). First, we analyze the properties the firing irregularity of our spike trains. This analysis reveals two features of the spiking activity in our set of motor cortical neurons: many neurons show temporally stable firing irregularity; the firing irregularity over the population of neurons, as measured using CV

_{2}, are widely distributed from 0.3 to 1. Second, we estimate the network parameters from the dynamics of rate and firing irregularities. For most analyzed neurons (78.5%), the dynamics of rate and irregularity is consistent with the neuron being part of an inhibition-dominated network.

### Firing properties of neurons

To study how the firing rate and the firing irregularity change during the task, their means are measured within sliding windows. For the firing irregularities, spike trains are preprocessed with the bias-suppression method; sliding windows which contain <20 spikes are discarded and not used for the data analysis (see methods). Figure 4*A* shows a typical example of the dynamics of rate and firing irregularity of a single neuron recorded in *monkey O* during all successful trials in one behavioral condition. Figure 4*A, top*, shows the spike raster plot aligned to MVT. Both mean firing rate and CV_{2} are presented in Fig. 4*A, bottom*. The firing rate (thick solid line) clearly shows a ramping up activity preceding and following movement onset. On the other hand, the firing irregularity does not seem to change much, and the temporal modulation of CV_{2} (thin solid line) does not strongly deviate from the confidence limit of 2σ_{CV2} region (shaded region). To compare the significance of the change level, we computed the psychophysical measure *d′* (Green and Swets 1966) of the highest and lowest value of the rate and the CV_{2}, where *d′* = |*x*_{H} − *x*_{L}|/[var(*x*_{H}) + var(*x*_{L})]^{1/2}. Here *x*_{H(L)} is equivalent to the highest (lowest) value of ν(*t*) or CV_{2}(*t*) in each trial, and window. var(*x*_{H(L)}) is the trial-to-trial variance of the firing rate and CV_{2} within the 100 ms window around the highest (lowest) value. The larger *d′*, the better one can discriminate the membership of a sample, i.e., whether a sample of mean rate (or mean CV_{2}) is derived from the lowest or highest value group. Figure 4*B* shows the histogram of *d′* for the firing rate (solid line) and the CV_{2} (dashed line) for all task conditions and neurons recorded from both monkeys. For most of the neurons, the *d′* for the firing rate is higher than that for the CV_{2}. This indicates that most of the neurons show large firing rate changes, but the firing irregularity is rather weakly modulated.

##### BROAD DISTRIBUTION OF CV_{2}.

Figure 4*E* shows the distribution of mean CV_{2} values averaged over four behavioral conditions across the population of neurons of the two monkeys (*n* = 164). Spike irregularities are widely distributed from quasi-regular (CV_{2} ∼0.3) to highly irregular (CV_{2} > 1). This holds true for the distribution of each monkey. In all recorded neurons, mean CV_{2} values are essentially independent of task conditions (Fig. 4*C*). Neurons recorded in different tasks also showed very similar values of firing irregularity in different tasks (Fig. 4*D*). Hence each neuron seems to be characterized by a CV_{2} value that is independent on the context of the task.

A possible explanation for the wide distribution of CV_{2}s is that distinct neuron classes (e.g., fast spiking interneurons and pyramidal neurons) could have different CV_{2} distributions. Because these neuron types are expected to have also distinct mean rates, we investigated the correlation between mean firing rate and firing irregularity across neurons. A rate-CV_{2} scatter plot is shown in Fig. 4*F*. Each point represents for a single neuron an average of time-averaged firing rate and CV_{2} from four task conditions. The correlation coefficient is 0.05, and the least squared error fit gives CV_{2} = 0.7 + 0.001ν. This indicates that in a population of neurons there is negligible correlation between the mean values of firing rate and CV_{2}.

### Estimating network parameters

In this section, we estimated the network parameters that can describe the observed dynamics of firing rate and firing irregularity. Basic assumptions are as follows; quasi-stationarity of the spike statistics during the sliding window (100 ms), the network is a sparse, randomly connected network of excitatory and inhibitory neurons; and neurons can be described as LIF neurons. The third assumption is justified by previous studies that have shown that an input-output relationship of neocortical neurons can well be described as LIF neurons (Arsiero et al. 2007; La Camera et al. 2006; Rauch et al. 2003) even though additional parameters such as adaptation improve the fitting. Here we use the standard current-based LIF model and fixed the single neuron model parameters as described in methods.

The random, sparsely connected network is characterized by three parameters. *1*) The inhibition/excitation balance *g* (IPSP/EPSP amplitude ratio). When *g* = 4 = 1*/*β, it gives the exact balance between the local, recurrent excitatory and inhibitory inputs, because we have assumed the ratio of the total number of excitatory and inhibitory neurons in a local network is 1/4. *2*) The EPSP amplitude *J* in mV. *3*) The number of recurrent excitatory connections per neuron *C*_{E} (the number of recurrent inhibitory synapses per neuron *C*_{I} = β*C*_{E} is always proportional to *C*_{E}).

A typical example of the estimated parameters and of the activity of a neuron exhibiting a high firing irregularity (CV_{2} > 0.9) is shown in Fig. 5*A*. Figure 5*A, top*, shows the dynamics of firing rate and firing irregularity. Figure 5*A, bottom*, shows the distribution of data points in the ν-CV_{2} plane. Contour plots (circles) around the data points represent error bars of [0.5, 1, 2] SE. The theoretical prediction lines from the estimated parameters are also overlaid. We estimated the network parameters for this neuron by using the method described in *Theory for estimating network parameters from single neuron statistics*. The estimated parameter that gives the global minimum of cost is (*g, J, C*_{E}) = (6.9, 0.27, 100) (thick gray line, cost per data: 1.343). However, we noticed that very different parameters also give very similar cost values, such as (*g, J, C*_{E}) = (5.8, 0.1, 1000) (dashed line, cost per data point: 1.345), (*g, J, C*_{E}) = (7.2, 0.8, 10) (dash-dotted line, cost per data point: 1.346) as shown in Fig. 5*A*. These parameter sets show similar cost values and the predicted theoretical lines are almost identical. This indicates that there are wide regions in parameter space that produce an essentially undistinguishable firing behavior. In other words, the firing rate-irregularity relationship is redundantly represented in the parameters of *g, J, C*_{E}.

Another typical example of neurons exhibiting quasi-regular spike trains (CV_{2} < 0.6) is shown in Fig. 5*B*. The estimated parameter that gives the global minimum is (*g, J, C*_{E}) = (5.4, 0.03, 1,000) (thick gray line, cost per data: 0.636). However, we also found that other parameters, e.g., (*g, J, C*_{E}) = (6.0, 0.08, 100), (dashed line, cost per data: 0.660) and (*g, J, C*_{E}) = (8.0, 0.1, 10), (dash-dotted line, cost per data: 0.714), also show a very similar firing rate-irregularity relationship. The parameter space is also degenerate in low CV_{2} neurons.

##### STRUCTURE OF THE SPACE OF PARAMETERS THAT MINIMIZES THE COST FUNCTION.

Why can more than one parameter set generate similar firing rate-CV_{2} dynamics? How are these parameters related to each other? To understand these questions, we studied the constraints on the parameters that can generate similar output firing statistics. The firing statistics in our model depends on mean μ and SD σ of the inputs to the neuron. We therefore ask the question, what are the network parameters that can generate given values of μ and σ? To answer this question, we rewrite *Eqs. 8* and *11* as
*C*_{E} is large, then to get a finite variance we need the EPSP amplitude *J* to scale roughly as 1/√*C*_{E} (van Vreeswijk and Sompolinsky 1996, 1998). Then to get a finite mean we need total excitatory and inhibitory inputs to roughly balance. Balancing excitatory and inhibitory inputs is only possible if the IPSP amplitude is strong enough
*g* > 4 because we set β = 0.25. This provides a first constraint on parameter space. This “balance” condition leads to the equation ν′_{ext} = *C*_{E}ν(*g*β − 1). Inserting this expression in *Eq. 37*, we find the following relationship between the input variance σ^{2} and the firing rate ν
^{2} and ν. For each neuron, this relationship should be characterized by a single constant *A*. Given a particular value of *A*, this equation provides us with a second constraint on parameter space. Different values of parameters *J, g*, and *C*_{E} that give the same value of *A* should give the same firing statistics. Note that this relationship describes a two-dimensional manifold in the three-dimensional (*J, g, C*_{E}) parameter space. Does this two-dimensional manifold fit with the one obtained from the fitting procedure? To answer this question, we plotted *Eq. 40* in *g* > 4 region for a fixed value of *A* (*A* = 3,025 for Fig. 6*A, A* = 225 for *B*) that fits the global minimum shown as an open circle. Interestingly, this equation obtained from simple balance arguments seems to fit the contours of the cost function very well, especially for *C*_{E} > = 100. For lower *C*_{E} (*C*_{E} = 10), the exact balance condition line does not fit the lowest cost region very well. This is consistent with the fact that the analytical argument is obtained in the large *C*_{E} limit and therefore cannot be expected to reproduce quantitatively the results of the fitting procedure in the low *C*_{E} region. In any case, this analysis shows that, to obtain a certain firing rate and firing irregularity, there are infinite numbers of network architectures with different values of *g, J*, and *C*_{E}, but these parameters are constrained in a specific way.

The constraint on the network parameters can be easily understood from the *Eq. 40*. This equation implies that for a given value of *C*_{E}, *J* is approximately inversely proportional to *g*. It also implies that when *C*_{E} increases, the magnitude of *J* scales as the square root of *C _{E}* (

*J*∝ 1/√

*C*

_{E}). Therefore as

*C*

_{E}increases, the low cost region tends to be confined to the lower values of

*J*region. The dependence of the best-fit parameters on the firing irregularity of the neurons can be also easily understood from

*Eq. 39*. In general, the higher

*A*, the larger becomes the variance of the input for a given same firing rate. Therefore the higher the constant

*A*, the higher the firing irregularity (CV

_{2}). If two neurons have the same level of EPSP size and IPSP/EPSP ratio but have different levels of firing irregularity, the neuron with the lower CV

_{2}value should have a smaller recurrent connectivity

*C*

_{E}because the lower the CV

_{2}, the smaller is the constant

*A*.

##### DISTRIBUTION OF ESTIMATED PARAMETERS OVER THE POPULATION.

We focused on two typical neurons in our data set, one irregular, the other more regular. We now show the distribution of the estimated parameters across the population of neurons. To understand how the best-fit parameters are distributed, we plot the parameter that gives a global minimum for each condition in Fig. 7. Estimated parameters are colored based on their mean CV_{2} to clarify the relationship between the firing irregularity and the estimated parameters.

When the number of recurrent connections per neuron, *C*_{E}, is small (*C*_{E} = 10), the estimated parameters *g* for irregular spiking neurons (CV_{2} > 0.9) are concentrated in the inhibition-dominated region (*g* > 4), but for quasi-regular spiking neurons, parameters providing a good fit extend in both the excitation and inhibition dominated region depending on EPSP size.

When the number of recurrent connections becomes large (*C*_{E} = 100), the local cortical circuits need to be inhibition dominated networks to account for the dynamics mainly for high CV_{2} neurons and the EPSP size *J* approaches more realistic values (*J* ≈ 0.2). For a neuron with a CV_{2} lower than 0.6, the estimated balance level (*g*) distributes mainly in inhibition-dominated but can still be slightly extended to excitation-dominated region. The estimated EPSP size *J* of quasi-regular neurons tends to be lower than that of irregular spiking neurons.

As the number of recurrent connections becomes large (*C*_{E} = 1,000), the EPSP size becomes smaller for a similar level of spike irregularity, and for almost all the range of irregularity, the estimated balance level is confined to the inhibition-dominated region. For all *C*_{E}, 78.5% (322/410) of estimated parameters the cost of which is <5 are classified as locally inhibition-dominated (*g* > 4). This result is consistent with the first constraint (*Eq. 38*) on the parameter space. The conclusion of inhibition-dominance does not change for a different set of single neuron parameters (Supplementary Fig. S4).

##### EFFECT OF THE SYNAPTIC PROPERTIES ON THE ESTIMATION RESULTS.

We found that the landscape of the cost-function is degenerate, but the estimated network parameter is mostly confined in the inhibition-dominated region. The higher CV_{2} neurons tend to have larger EPSP size *J* for a fixed number of presynaptic neurons *C*_{E}. The neuron model we have used was the simple current-based LIF neuron that allows us to clearly understand the estimation process. However, real synaptic currents have rich biophysical properties, such as voltage-dependence and short-term plasticity. To understand how our results are affected by incorporating additional biophysical synaptic features to our model, we applied our estimation method to two types of LIF neuron models—with synapses with STD and with conductance-based synapses—and studied their effect on our estimation results.

#### STD in excitatory synapses.

One significant feature of the neocortical neurons is their activity-dependent short-term synaptic plasticity (Thomson et al. 1993). Because short-term plasticity in recurrent synapses nonlinearly modulates the input statistics, its contribution to the firing statistics is not straight forward. Here we used our parameter estimation method to the LIF neuron model with excitatory STD synapses. In the STD model, the average EPSP size is *Ju*〈*x*〉 where *J* is the “absolute” synaptic strength, *u* is the fraction of released vesicles from the available resources, and 〈*x*〉 is the average fraction of available vesicles (see methods). 〈*x*〉 in turn depends on the firing rate of presynaptic neurons and is <0.2 for the parameters chosen here and if the firing rate is more than 10 Hz. Parameter values of *u* = 0.5 and 〈*x*〉 ∼ 0.2 make the effective EPSP amplitude *Ju*〈*x*〉 ∼10 times smaller than the absolute synaptic strength *J*. This indicates that, for a synapse with STD to achieve the same EPSP size as a synapse with no STD, the absolute value of EPSP size *J* must be nearly 10 times larger than that of a synapse with no STD. This also means that the parameter *g* is no longer a good measure for the balance between excitation and inhibition. This balance is now *g/u*〈*x*〉, where 〈*x*〉 is firing rate dependent. We therefore plotted the best-fit parameters in the *Ju*〈*x*〉 – *g/u*〈*x*〉 plane in Fig. 8*A*, as well as in the *Ju* – *g/u*〈*x*〉 plane in *B*. For each neuron, we computed 〈*x*〉 from its median firing rate. The effective EPSP size during the task can be approximated by *Ju*〈*x*〉, but we also used *Ju*, which corresponds to the EPSP size in absence of prior presynaptic activity and is therefore directly comparable to experimental values observed in slice conditions in which neurons are essentially silent.

The distribution of the estimated parameters shows significant differences with the distribution of the estimated parameters in the network with no STD (compare Figs. 7 and 8, *A* and *B*). First, all neurons are no longer confined to the inhibition-dominated region—rather, neurons tends to be now clustered around the balance point. Neurons with high CV_{2} still tend to be in the inhibition-dominated region, while many neurons with low CV_{2} are in the excitation-dominated region. Second, best-fit parameters tend to have much larger EPSP amplitudes compared with the case with no STD (compare Figs. 7 and 8*B*) due to the fact that in an active network STD reduces the effective EPSP size by a significant factor. A feature that is qualitatively unchanged is that high CV_{2} neurons tend to have higher *J*, larger *g* values.

#### Conductance-based synapses.

We turn now to another biophysical property of synapses, the voltage dependence of the synaptic currents. Modeling studies have shown that the LIF neuron model with conductance-based synapses does not require a strict balance between excitation and inhibition (Meffin et al. 2004) to achieve high irregularity of firing. We used conductance-based synapses in LIF neuron models and estimated network parameters for the same set of neuron data (see the parameters in methods). The estimated parameters showed a qualitatively similar distribution as the current-based LIF neuron network (Fig. 8*C*). The estimated parameters are shown in the effective EPSP size [*a*_{E}(*V*_{r} − *V*_{E})] and the effective ratio of IPSP to EPSP [*g* = (*a*_{I}|V_{r} − *V*_{I}|)/*a*_{E}(*V*_{r} − *V*_{E})] measured at the reset potential V_{r}. The high CV_{2} neurons tend to have higher effective EPSP values and higher *g* values for a fixed *C*_{E}. Therefore our estimation results are consistent in both the current- and conductance-based models. In Fig. 8*D*, we also showed the distribution of the effective time constants of each neuron computed from its median firing rate. The high CV_{2} neurons tend to have lower effective membrane time constant (compare Fig. 8, *C* and *D*) because they have larger EPSP and IPSP sizes (see *Eq. 33*), and given that there is no significant correlation between mean firing rate and CV_{2} (Fig. 4*F*).

#### Stochastic synapses.

Another prominent property of synapses is the variability of the amplitude of postsynaptic potentials from spike to spike. This variability can be quantified by the CV (SD divided by the mean). For connections between cortical pyramidal cells, the CV is typically of order 0.5 (Markram et al. 1997). How does this variability affect our results? It turns out that we can answer this question in a very simple way, by evaluating how the variability affects the variance of the fluctuations due to random synaptic inputs. With stochastic synapses, the variance becomes
_{syn} is the CV of PSPs, which is assumed here for simplicity to be the same for all types of synapses. The next question is the impact of this variability on the best-fit parameters. This question can be answered by inspecting *Eq. 40*. To satisfy this equation for a given value of *A*, we see that an increase in CV_{syn} should lead to a corresponding decrease in *J* and/or *g*. In the case of *g*, the constraint *g* > 1/β still holds, so *g* cannot be decreased below this value. In practice, for realistic values of CV_{syn}∼0.5, the decrease in typical PSP strengths is small (∼10%).

#### Finite synaptic decay times.

In the model studied here, postsynaptic currents are assumed to be delta functions. Including a finite decay time of postsynaptic currents complicates significantly the calculation of both firing rates and CV of the ISIs as a function of the statistics of the input. In fact, only expressions in the limit in which synaptic decay time is either very short or very large compared with the membrane time constant are known (see e.g., Brunel and Sergi 1998; Fourcaud and Brunel 2002; Moreno et al. 2002). However, we can understand the qualitative effect of such time constants by looking at the results of these previous studies. In practice, incorporating a finite decay time constant leads to an effective decrease in the SD of the noise (see e.g., Fourcaud and Brunel 2002). So the result of incorporating a synaptic decay time constant should be the opposite of incorporating stochasticity of synaptic transmission. For example, for a fixed value of *g* and *C*_{E}, the estimated EPSP amplitude *J* should increase as the synaptic decay time constant increases.

## DISCUSSION

In this paper, we analyzed single-unit recordings from motor cortex of two behaving monkeys. These data show two interesting features of the firing statistics of these neurons: in many neurons, the firing irregularity, as measured by CV_{2}, fluctuates very little in time and across task conditions even though the firing rates exhibit pronounced modulations and the distribution of the CV_{2} values of all cells is rather broad, ranging from 0.3 to >1. We then explored what a specific temporal rate-CV_{2} relationship tells us about the structure of the local network in which the neuron is embedded. We introduced a method to find the parameters of a randomly connected network of excitatory and inhibitory LIF neurons that minimize the distance between data and model. This method indicates that most recorded cells in motor cortex live in local networks that are dominated by inhibition when considering a simple current-based model. This result can be understood with the following intuitive argument: if rate modulations are simply due to changes in external excitatory inputs, then one would expect that the CV_{2} decreases monotonically with rate. This is because when the external inputs are low, firing occurs due to fluctuations in the inputs and is approximately Poisson with a CV_{2} close to 1. On the other hand, when the external inputs are high, firing becomes more deterministic with a CV_{2} that decreases with input strength. In fact, in cortical slices, the CV of both pyramidal cells and interneurons decreases monotonically with the mean injected current, when the variance of the noise is fixed (La Camera et al. 2006; Rauch et al. 2003). This argument no longer holds true in networks dominated by inhibition because the increase in external inputs tends to be counterbalanced by an increase in the firing rate of inhibitory neurons. In such networks, an increase in firing rate can be induced by a combination of an increase in the mean inputs and an increase in the variance of the fluctuations (see also Renart et al. 2006). Whereas the increase in the mean inputs tends to decrease the CV_{2}, the increase in the variance has the opposite effect of decreasing the CV_{2}. Therefore there exists a combination of the two effects that leads to an approximately constant firing irregularity as the rate changes (Miura et al. 2007).

### Weak modulation in firing irregularity

High irregularity is a hallmark of electrophysiological recordings in cortex in vivo (Compte et al. 2003; Shinomoto et al. 2003; Softky and Koch 1993). Shimokawa and Shinomoto (2009) have studied the instantaneous firing irregularity based on Bayesian estimation and concluded that the firing irregularity fluctuates only weakly for neurons in cortical areas V1 and MT compared with those in the LGN of the thalamus. Here we also find that the CV_{2} is approximately constant for most cells throughout the whole trial. Our estimation method then suggests this weak modulation in firing irregularity could be due to recurrent dynamics of a local network dominated by inhibition. This is consistent with a recent in vitro study that showed that a balanced input can decouple the firing rate from the firing irregularity (Miura et al. 2007). In other words, when the excitatory and inhibitory conductances are proportionally changed, firing irregularity can be fixed but the firing rate can be modulated independently (see also Ponce-Alvarez et al. 2010).

There are, however, alternative explanations to the weakly modulated irregularity. Barbieri and Brunel (2007) showed that the CV of a simple LIF neuron becomes weakly dependent on the rate in a wide range of rates when the reset potential is very close to threshold and the membrane time constant is short. With these parameter settings, an approximately constant CV could be explained by single cell properties, rather than by network properties. The CV versus rate relationship of cortical cells measured in vitro shows, however, a pronounced decrease of CV as a function of rate, which seems at odds with the single cell scenario (La Camera et al. 2006; Rauch et al. 2003).

### Wide distribution of irregularity measures

The results of our analysis suggest that the distribution of the CV_{2} is wide, ranging from 0.3 to >1. This is consistent with the analysis of spiking activities recorded in other areas (Shinomoto et al. 2003, 2005, 2009). Again, one might think about two possible explanations to account for this wide distribution; network or single cell. If all recorded cells have the same characteristics, our analysis suggests that there is a wide range of coupling strengths in local networks of motor cortex and that neurons with low CV_{2} could be part of networks that are relatively weakly coupled (either with small numbers of recurrent connections per neuron or with very weak EPSP/IPSP amplitudes). Another “network” explanation is provided by Shinomoto et al. (2005), who showed that in area TE layer V–VI neurons tend to have a lower firing irregularity than layer II–III neurons. In both scenarios, the network connectivity would be the source of the observed wide distribution of firing irregularities. Another possibility is that there are several classes of neurons with distinct intrinsic properties, each of which being characterized by a distinct CV.

### Comparison of best-fit parameters with estimates from the literature

Our method provides constraints on parameters characterizing the network architecture: *J* (the EPSP size), *C*_{E} (the number of excitatory recurrent connections per neuron), and *g* (the IPSP/EPSP ratio). The strongest constraint is on *g*: our method predicts that for most recorded neurons, *g* > 4 (inhibition dominance) when using a simple current-based model. The second constraint can give us an upper bound on *J* given a particular value of *C*_{E}. The next question is whether the best-fit parameters are compatible with known experimental data?

The most robust finding of our estimation method is that the IPSP to EPSP ratio should be large enough to compensate for the fact that there are less inhibitory neurons than (excitatory) pyramidal neurons in cortex. This is in agreement with the fact that inhibitory synapses are located mainly on the cell body and, thus the size of their synaptic potential is less attenuated than that of excitatory synapses that are mainly located on the dendrites (Braitenberg and Schüz 1998). Furthermore, in our network interneurons are constrained to have the same firing rate as pyramidal cells. If interneurons have higher firing rates (La Camera et al. 2006; McCormick et al. 1985; Wilson et al. 1994), then the constraint on *g* can be alleviated (Amit and Brunel 1997). Incorporating STD, however, leads to a large number of neurons (especially the low CV ones) to move to the excitation-dominated region.

Our estimate of the EPSP size depends strongly on the number, *C*_{E}, of excitatory recurrent connections per neurons and is not very constrained. This is because for a fixed IPSP to EPSP ratio, the same firing rate and firing irregularity can be achieved for widely different values of *J* and *C*_{E}, provided *Eq. 40* is satisfied. In a model with no synaptic STD, for realistic values of *C*_{E} (≥1,000) (Braitenberg and Schüz 1998), our average EPSP sizes are typically small, on the order of 0.1- 0.2 mV, for high CV neurons, and even smaller (<0.1 mV) for low CV neurons (see Fig. 7). Inserting STD in the model leads to much higher EPSP sizes because STD leads to strong decrease of the effective EPSP size. In vitro studies reported average EPSP sizes in the range [0.1–1 mV] in neocortex (Holmgren et al. 2003; Markram et al. 1997; Mason et al. 1991; Sjöström et al. 2001). EPSP sizes have been reported to be in the order of 0.1–0.6 mV in monkey motor cortex in vivo (Matsumura et al. 1996). These numbers are consistent with our estimates provided STD is included in the model.

We have investigated a network of LIF neurons because of analytical tractability. Interestingly, both the *f-I* curve and *CV-I* curve of cortical pyramidal cells can be well fitted by the corresponding curves of integrate-and-fire neurons in vitro (Rauch et al. 2003). Exponential integrate-and-fire (EIF) neurons (Fourcaud-Trocmé et al. 2003) have been shown to provide a better fit of the timing of spikes of cortical pyramidal cells on injection of random fluctuating inputs (Badel et al. 2008), but both the *f-I* curve and *CV-I* curve of EIF neurons are qualitatively similar to those of LIF neurons (Fourcaud-Trocmé et al. 2003; Brunel, unpublished observations), so we would not expect any qualitative difference between LIF and EIF neurons as far as the results reported here are concerned.

### Functional implications

The surprising result that a high proportion of the analyzed neurons from motor cortex of the behaving monkey seem to be embedded in a network dominated by inhibition is in agreement with the notion that during movement preparation (neuronal activity was recorded mainly during movement preparation) the corticospinal motor output must be inhibited to prevent its premature recruitment. This preparatory activity goes in line with a wide spread oscillatory activity in the beta range of the LFP in motor cortical areas (MacKay 2005; Murthy and Fetz 1992), which stops abruptly with movement onset. Interestingly, delay-related oscillatory activity is often associated with an inhibitory behavior as well as expectancy, an activity that closely corresponds to the withholding of the movement during the delay period (for a review, see MacKay 2005; see also Riehle et al. 2006). It is only at the end of the delay that the inhibition of corticospinal neurons would be lifted and a local wave of excitatory input would briskly ignite a descending motor volley, corresponding to the small networks of rather regular spiking neurons.

## GRANTS

K. Hamaguchi was supported by Japan Society for the Promotion of Science Research Fellowship for Young Scientists. A. Riehle and N. Brunel were supported by the French Government (ANR-NEUR-05-045-01).

## DISCLOSURES

No conflicts of interest, financial or otherwise, are declared by the author(s).

## ACKNOWLEDGMENTS

Data were collected in the lab of A. Riehle mainly by Sébastien Roux and Franck Grammont. We thank M. Martin for animal care, B. E. Kilavik and A. Ponce-Alvarez for many helpful discussions; and G. Mongillo for clarifications about the short-term depression model.

Present address of K. Hamaguchi: Mooney Lab, Bryan Research Bld., Duke University Medical Center, Research Dr., Durham, NC 27710.

## Footnotes

↵1 The online version of this article contains supplemental data.

- Copyright © 2011 The American Physiological Society