Chaotic Balanced State in a Model of Cortical Circuits

C. van Vreeswijk
H. Sompolinsky
Racah Institute of Physics and Center for Neural Computation, Hebrew University, Jerusalem, 91904 Israel

The nature and origin of the temporal irregularity in the electrical activity of cortical neurons in vivo are not well understood. We consider the hypothesis that this irregularity is due to a balance of excitatory and inhibitory currents into the cortical cells. We study a network model with excitatory and inhibitory populations of simple binary units. The internal feedback is mediated by relatively large synaptic strengths, so that the magnitude of the total excitatory and inhibitory feedback is much larger than the neuronal threshold. The connectivity is random and sparse. The mean number of connections per unit is large, though small compared to the total number of cells in the network. The netw ork also receives a large, temporally regular input from external sources. We present an analytical solution of the mean-field theory of this model, which is exact in the limit of large network size. This theory reveals a new cooperative stationary state of large networks, which w e term a balanced state. In this state, a balance between the excitatory and inhibitory inputs emerges dynamically for a wide range of parameters, resulting in a net input whose temporal fluctuations are of the same order as its mean. The internal synaptic inputs act as a strong negative feedback, which linearizes the population responses to the external drive despite the strong nonlinearity of the individual cells. This feedback also greatly stabilizes the system's state and enables it to track a time-dependent input on time scales much shorter than the time constant of a single cell. The spatiotemporal statistics of the balanced state are calculated. It is shown that the autocorrelations decay on a short time scale, yielding an approximate Poissonian temporal statistics. The activity levels of single cells are broadly distributed, and their distribution exhibits a skewed shape with a long power-law tail. The chaotic nature of the balanced state is revealed by showing that the evolution of the microscopic state of the network is extremely sensitive to small deviations in its initial conditions. The balanced state generated by the sparse, strong connections is an asynchronous chaotic state. It is accompanied by weak spatial cross-correlations, the strength of which vanishes in the limit of large network size. This is in contrast to the synchronized chaotic states exhibited by more conventional network models with high connectivity of weak synapses.

1 Introduction

The firing patterns of neurons in the cortex of intact animals often exhibit a strong degree of temporal irregularity. This can be seen by the broad interspike interval histograms (ISI) of cortical neurons, which are typically close to those generated by a Poisson process with a short refractory period (Abeles, 1991; Bair, Koch, Newsome, & Britten, 1994; Burns & Webb, 1976; Douglas, Martin, & Whitteridge, 1991; Softky & Koch, 1993). The irregular neuronal dynamics is also manifested in intracellular recordings of the membrane potential, which exhibit strong temporal fluctuations. One of the long-standing problems in cortical dynamics is understanding the origin of this irregularity and its computational implications (Douglas & Martin, 1991; Ferster & Jagadeesh, 1992). In vitro experiments show that cortical neurons fire in a relatively regular fashion when they are injected with a constant current. Thus, the irregularity of the in vivo neuronal activity must be due to fluctuations in their synaptic input (Holt, Softky, Koch, & Douglas, 1996; Mainen & Sejnowski, 1995). These fluctuations may be due to variations in the intensity of the sensory stimuli or may result from the stochastic action of synapses. However, since cortical cells have thousands of synaptic contacts, one would expect that the summation of the synaptic inputs at the soma averages out most of the fluctuations in the synaptic input and yields a membrane potential with only a small residual fluctuation. This is a particularly difficult issue in conditions where the cortex is vigorously active so that the cell receives many synaptic inputs within a single integration time constant (Holt et al., 1996; Softky & Koch, 1993). One possible resolution of this problem is to assume that the fluctuating synaptic inputs are substantially correlated and therefore are not averaged out. Indeed, the spike trains of pairs of neurons in cortex and in thalamus are often correlated in a relatively narrow time scale (of the order of 10 msec ) (Abeles, 1991; Gray & Singer, 1989; Perkel, Gerstein, & Moore, 1967a, 1967b; Vaadia et al., 1995). However, the observed size of these correlations indicates that in general, only a small fraction of the neuronal activity is tightly correlated. Another possibility, addressed in this article, is that although the inputs to a cell are only weakly correlated, the cell is sensitive to the residual correlations in the somatic potential.

Several mechanisms that generate enhanced sensitivity of a cell to small fluctuations in its potential have been explored (Bell, Mainen, Tsodyks, & Sejnowski, 1994; Ermentrout & Gutkin, in press; Gerstein & Mandelbrot, 1964; Shadlen & Newsome, 1994, 1995; Softky, 1995; Troyer & Miller, 1997). One possibility is that the excitatory and inhibitory inputs to a cortical cell are balanced in a way that leaves the cell's average potential close to threshold, and its firing pattern is therefore susceptible to small fluctuations. An interesting question is what might be the mechanism that leads to this balance. An interesting recent study (Tsodyks & Sejnowski, 1995) explored the possible involvement of local cortical dynamics in balancing excitation and inhibition. This numerical study invoked a strong stochasticity in the synaptic action in the form of a large failure probability. In a related study (A mit & Brunel, 1997a, 1997b), the variability in the network activity is at least partially due to fluctuating external inputs to the local network. In addition, neither study properly addresses important issues concerning the behavior of the models and the robustness of their variability as the network size is scaled up.

In this article, we investigate the hypothesis that the intrinsic deterministic dynamics of local cortical networks is sufficient to generate strong variability in the neuronal firing patterns. Neuronal dynamics is highly nonlinear; hence it may seem natural to expect that neuronal networks with deterministic dynamics will exhibit chaotic behavior. However, studies of simple models of large networks with a high degree of connectivity (Abbott & van Vreeswijk, 1993; Gerstner & van Hemmen, 1993; Grannan, Kleinfeld, & Sompolinsky, 1992; Hansel, Mato, & Meunier, 1995; van Vreeswijk, 1996; Wilson &Cowan, 1972; Tsodyks, Mitkov, &Sompolinsky, 1993) reveal that in the absence of external sources of strong stochastic noise, they tend to settle into temporally ordered states of tonic firing or oscillations. Recent extensive numerical study (Bush & Douglas, 1991; Hansel & Sompolinsky, 1992, 1996) of a model of local circuits in visual cortex with realistic conductance-based dynamics has shown the existence of parameter regimes in which these networks exhibit strongly irregularstates, denoted as synchronized chaotic states. These chaotic states are generated by the emergence of strong synchrony in the fluctuating activity of different neurons, which consistently generates a strong fluctuating feedback to each cell. Thus, this is a network realization of the scenario of correlated synaptic inputs. The resulting patterns of activity show strongly synchronized bursting patterns, tightly timed by the common inhibitory feedback. Although bursting patterns are sometimes observed in cortical networks, these synchronized chaotic states are hard to reconcile with the Poisson-like weakly correlated firing patterns, commonly observed in cortex.

In this work, we explore the possibility that local networks with intrinsic dynamics evolve toward states that are characterized by strong chaos in conjunction with weak cross-correlations, through the mechanism of balancing between excitation and inhibition. This possibility raises several questions:

• What are the conditions under which a network evolves to a state in which the excitatory and inhibitory inputs are balanced?
• What are the characteristics of this balanced state? Does the balanced state represent a cooperative state that is qualitatively distinct from the synchronized chaotic state?
• What are the functional advantages of the balanced state?

We study these questions using a network model with the simplified dynamics of binary elements. The architecture consists of excitatory and inhibitory populations connected by sparse random connections.

Indeed, our analytical solution of the model in the limit of large network size shows that in a broad range of parameters, the network settles into a stable balanced state. An interesting feature of our theory is that it goes far beyond calculating the properties of the macroscopic order parameters. The theory yields a complete statistical characterization of the balanced state. It shows that the balanced state is associated with a strong Poisson-like firing pattern and also with abroad inhomogeneity in the average rates of individual neurons. Finally, we address the possible functional implications of the balanced state by showing that the network is capable of fast tracking of temporal changes in the external input to the network.

In section 2, we present the model's dynamics and architecture. Section 3 presents the mean-field dynamic equations of the evolution in time of the two macroscopic order parameters, which are the rates of activity of the two subpopulations.

Section 6 addresses the stability of the balanced state. It shows that there is a comfortable parameter regime where the balanced state is stable. We also discuss what happens to the network when the balanced fixed point is unstable. Section 7 considers the effect of inhomogeneity in the local thresholds. We show that in the presence of inhomogeneity, the distribution of rates acquires a characteristic skewed shape with a long tail, qualitatively similar to the observed distribution of rates in populations of neurons in the cortex of behaving monkeys. In section 8, we evaluate the sensitivity of the temporal fluctuations in the local instantaneous activities to a small change in the initial condition. We conclude that a small change in the initial condition leads rapidly to a complete loss of memory of the unperturbed initial conditions. Thus, our network shows the main characteristics of chaotic systems. Section 9 studies the dynamical response of the system to dynamic changes in the external input and shows the fast tracking capabilities of the network. In section 10 we discuss the results and some open issues. Details of the theory are outlined in appendixes A and B .

A preliminary report on this work was published in van Vreeswijk and Sompolinsky (1996).

Figure 1: A schematic representation of the network architecture. Excitatory connections are shown as open circles; inhibitory ones as filled circles.

2 The Model

We consider a network of $\mathrm{N}_{1}$ excitatory and $\mathrm{N}_{2}$ inhibitory neurons. The network also receives input from excitatory neurons outside it(see Figure 1). We will use either the subscript 1 or E to denote the excitatory population and 2 or I for the inhibitory one. The pattern of connections is random but fixed in time.

We will call K the connectivity index. The state of each neuron is described by a binary variable $\sigma$. The value $\sigma=0(\sigma=1)$ corresponds to a quiescent (active) state. The network has an asynchronous dynamics where only one neuron updates its state at any given time. The updated state of the updating neuron at time $t$ is

$$\sigma_{\mathrm{k}}^{\mathrm{i}}(\mathrm{t})=\Theta\left(\mathrm{u}_{\mathrm{k}}^{\mathrm{i}}(\mathrm{t})\right),$$

where $\Theta(\mathrm{x})$ is the Heaviside function, $\Theta(\mathrm{x})=0$ for $\mathrm{x} \leq 0$ and $\Theta(\mathrm{x})=1$ for $\mathrm{x}>0$. The total synaptic input, $\mathrm{u}_{\mathrm{k}}^{\mathrm{i}}$ to the neuron, relative to the threshold, $\theta_{\mathrm{k}}$, at time t is

$$\mathrm{u}_{\mathrm{k}}^{\mathrm{i}}(\mathrm{t})=\sum_{\mathrm{l}=1}^{2} \sum_{\mathrm{j}=1}^{\mathrm{N}_{\mathrm{l}}} \mathrm{~J}_{\mathrm{kl}}^{\mathrm{ij}} \sigma_{\mathrm{l}}^{\mathrm{j}}(\mathrm{t})+\mathrm{u}_{\mathrm{k}}^{0}-\theta_{\mathrm{k}},$$

where $\mathrm{u}_{\mathrm{k}}^{0}$ denotes the constant external input to the kth population. As explained in appendix B , the precise definition of the order of updates is not essential. One model is a stochastic model in which each neuron updates its state at time intervals that have Poisson statistics. This model is the simplest to analyze. However, it has the drawback that it introduces a stochastic element (the random choice of the updating neuron). An alternative model is a fully deterministic one in which each neuron updates its state at equally spaced times where the time between updates is different for each neuron. We show in appendix B that the two models have the same mean-field equations. In both cases, the mean interval between consecutive updates of a neuron of the kth population is $\tau_{\mathrm{k}}$. We will use time units such that $\tau_{\mathrm{E}}=1$ so that the only independent time parameter is $\tau \equiv \tau_{\mathrm{I}}$.

To correspond to point processes, we define a spike as the transition from thepassive ( 0 ) to the active ( 1 ) state. Note that the firing rate, $\mathrm{r}_{\mathrm{k}^{\mathrm{k}}}^{\mathrm{i}}$, of neuron in population k is different from the average value, $\mathrm{m}_{\mathrm{k}}^{\mathrm{i}}(\mathrm{t})$, of $\sigma_{\mathrm{k}}^{\mathrm{i}}$ because before the cell can spike, it has to update to the passive state. However, if neuron i of the kth population updates to the active state in two consecutive updates, the synapses projecting from this cell stay active after the second update, even though no new spike is emitted. However, if $\mathrm{m}_{\mathrm{k}}^{\mathrm{i}}$, which we will call the activity rate, is small, the probability of two consecutive updates to the active state is low, and thus for small $\mathrm{m}_{\mathrm{k}}^{\mathrm{i}}$, the activity rate and the firing rate are nearly equal. Indeed, if we assume that at each update the probability of being in the active state is $\mathrm{m}_{\mathrm{k}}^{\mathrm{i}}$ (very nearly true in this model for low rates, as shown in section 5.3), the firing rate is given by $\mathrm{r}_{\mathrm{k}}^{\mathrm{i}}=\mathrm{m}_{\mathrm{k}}^{\mathrm{i}}\left(1-\mathrm{m}_{\mathrm{k}}^{\mathrm{i}}\right) / \tau_{\mathrm{k}}$.

A central ingredient of our model is the assumption that the total excitatory feedback current and the total inhibitory current into a cell are large compared to the neuronal threshold. We model this by choosing thresholds $\theta_{\mathrm{k}}$ that are of order 1 and by assuming that the strength of individual synapses is of order $1 / \sqrt{\mathrm{K}}$, that is, the coefficients $\mathrm{J}_{\mathrm{kl}}$ are of order unity. Furthermore, as will be seen later, it is crucial that the excitatory inputs from the external sources too are large compared to the threshold. This is modeled by denoting these inputs as

$$\mathrm{u}_{\mathrm{k}}^{0}=\mathrm{E}_{\mathrm{k}} \mathrm{~m}_{0} \sqrt{\mathrm{~K}} \quad \mathrm{k}=1,2,$$

where $\mathrm{E}_{\mathrm{k}}$ is of order unity and $0<\mathrm{m}_{0}<1$ represents the mean activity of the external neurons. We will use the notation

$$\mathrm{E}_{1}=\mathrm{E}$$

for the external input to the excitatory population and

$$E_{2}=I$$

for the external input to the inhibitory neurons. We assume that the external input is temporally regular.

Since the model neurons are threshold elements, the absolute scale of $u_{i}^{k}$ is irrelevant. We therefore set

$$\mathrm{J}_{\mathrm{EE}}=\mathrm{J}_{\mathrm{IE}}=1,$$

so that the only connection parameters from the network are the inhibitory and external ones. We will denote the former as

$$\mathrm{J}_{\mathrm{E}} \equiv-\mathrm{J}_{\mathrm{EI}} ; \mathrm{JI}_{\mathrm{I}} \equiv-\mathrm{J}_{\mathrm{II}},$$

where $\mathrm{J}_{\mathrm{I}}, \mathrm{J}_{\mathrm{E}}>0$.

3 Mean-Field Equations for Population Rates

The dynamics of our model can be described by mean-field theory, which is exact in the limit of large $\mathrm{N}_{\mathrm{k}}$. To define this limit, we assume that $\mathrm{N}_{\mathrm{I}} / \mathrm{N}_{\mathrm{E}}$ is held fixed as the network size $\mathrm{N}=\mathrm{N}_{\mathrm{E}}+\mathrm{N}_{\mathrm{I}}$ grows. The nature of the meanfield theory depends on the assumed relationship between the network size and the connectivity index. Conventional mean-field theory assumes that the networks are fully connected, defined here to mean that $\mathrm{K} / \mathrm{N}$ is fixed as $\mathrm{N} \rightarrow \infty$. Here we assume sparse connectivity defined by assuming that K is fixed as N grows. We are primarily interested in temporal variability that is present in highly connected networks, which are either fully connected or sparsely connected with large connectivity index K. Therefore we will focus on the case of large K . Technically, we will first take the limit $\mathrm{N} \rightarrow \infty$ and then the limit $\mathrm{K} \rightarrow \infty$. In reality, networks have a large fixed size and connectivity, so that the distinction between full and sparse connectivity may be problematic. Nevertheless, the sparse limit is appropriate as long as

$$1 \ll \mathrm{~K} \ll \mathrm{~N}_{\mathrm{k}}, \quad \mathrm{k}=1,2 .$$

The mean-field theory of our model for arbitrary fixed K is presented in appendix A. Taking the large K limit provides a substantial simplification of the mean-field equations. In this limit, most of the properties of the system can be expressed in terms of the first and second moments of the neuronal activity levels, as will be shown here and in the following sections. We first consider the population-averaged firing rates of the excitatory and inhibitory cells as

$$\mathrm{m}_{\mathrm{k}}(\mathrm{t})=\left[\sigma_{\mathrm{k}}^{\mathrm{i}}(\mathrm{t})\right]=\frac{1}{\mathrm{~N}_{\mathrm{k}}} \sum_{\mathrm{i}=1}^{\mathrm{N}_{\mathrm{k}}} \sigma_{\mathrm{k}}^{\mathrm{i}}(\mathrm{t}), \quad \mathrm{k}=1,2,$$

where [. . .] denotes a population average. In appendix A we show that the average activities satisfy in the large K limit

$$\tau_{\mathrm{k}} \frac{\mathrm{~d}}{\mathrm{dt}} \mathrm{~m}_{\mathrm{k}}(\mathrm{t})=-\mathrm{m}_{\mathrm{k}}(\mathrm{t})+\mathrm{H}\left(\frac{-\mathrm{u}_{\mathrm{k}}}{\sqrt{\alpha_{\mathrm{k}}}}\right) .$$

Here H is the complementary error function,

$$H(z) \equiv \int_{z}^{\infty} \frac{d x}{\sqrt{2 \pi}} e^{-x^{2} / 2},$$

Figure 2: Complementary error function $\mathrm{H}(\mathrm{x})$. The error function varies sigmoidally from 1 for $\mathrm{x} \rightarrow-\infty$ to 0 at $\mathrm{x} \rightarrow \infty$.

shown in Figure 2. The quantities $u_{k}(t)$ and $\alpha_{k}(t)$ are

$$\mathrm{u}_{\mathrm{k}}(\mathrm{t})=\sqrt{\mathrm{K}}\left(\sum_{\mathrm{l}=1}^{2} \mathrm{~J}_{\mathrm{kl}} \mathrm{~m}_{\mathrm{l}}(\mathrm{t})+\mathrm{E}_{\mathrm{k}} \mathrm{~m}_{0}\right)-\theta_{\mathrm{k}}$$

and

$$\alpha_{\mathrm{k}}(\mathrm{t})=\sum_{\mathrm{l}=1}^{2}\left(\mathrm{~J}_{\mathrm{kl}}\right)^{2} \mathrm{~m}_{\mathrm{l}}(\mathrm{t})$$

respectively. Equation 3.5 denotes the population average of the total input to a neuron in the kth population, relative to threshold. Equation 3.6 de notes the variance of this input. Note that the external population does not contribute to variance because we assumed that the input is the same for all the neurons in a population.

In the case of a constant external input, $\mathrm{m}_{0}(\mathrm{t})=\mathrm{m}_{0}$, the network settles into a state in which the average activities are constant, $\mathrm{m}_{\mathrm{k}}(\mathrm{t})=\mathrm{m}_{\mathrm{k}}$, given by the stable fixed points of equation 3.3,

$$\mathrm{m}_{\mathrm{k}}=\mathrm{H}\left(\frac{-\mathrm{u}_{\mathrm{k}}}{\sqrt{\alpha_{\mathrm{k}}}}\right)$$

where the mean inputs are

$$\begin{aligned} & \mathrm{u}_{\mathrm{E}}=\left(\mathrm{Em}_{0}+\mathrm{m}_{\mathrm{E}}-\mathrm{J}_{\mathrm{E}} \mathrm{~m}_{\mathrm{I}}\right) \sqrt{\mathrm{K}}-\theta_{\mathrm{E}}, \\ & \mathrm{u}_{\mathrm{I}}=\left(\mathrm{Im}_{0}+\mathrm{m}_{\mathrm{E}}-\mathrm{J}_{\mathrm{I}} \mathrm{~m}_{\mathrm{I}}\right) \sqrt{\mathrm{K}}-\theta_{\mathrm{I}} . \end{aligned}$$

The variance of the inputs is

$$\alpha_{\mathrm{k}}=\mathrm{m}_{\mathrm{E}}+\mathrm{J}_{\mathrm{k}}^{2} \mathrm{~m}_{\mathrm{I}} .$$

Equation 3.3 reflects the fact that the instantaneous input to each neuron $u_{i}^{k}(t)$ fluctuates across the population of neurons, and these fluctuations obey a gaussian statistics in the large K limit. The expressions for the mean and variance of the input to a cell can be derived in the large K limit by the following arguments. The population average inputs is

$$\mathrm{u}_{\mathrm{k}}(\mathrm{t})=\left[\mathrm{u}_{\mathrm{k}}^{\mathrm{i}}(\mathrm{t})\right]=\sum_{\mathrm{l}=1}^{2} \sum_{\mathrm{j}=1}^{\mathrm{N}_{\mathrm{l}}}\left[\mathrm{~J}_{\mathrm{kl}}^{\mathrm{ij}}\right]\left[\sigma_{\mathrm{l}}^{\mathrm{j}}(\mathrm{t})\right]+\mathrm{u}_{\mathrm{k}}^{0}-\theta_{\mathrm{k}} .$$

The population average $\left[\mathrm{J}_{\mathrm{kl}}^{\mathrm{ij}}\right]$ is equivalent to a quenched average over the random connectivity and is therefore equal to $\mathrm{J}_{\mathrm{kl}} \sqrt{\mathrm{K}} / \mathrm{N}$, yielding equation 3.5. Note that on the right-hand side (r.h.s.) of equation 3.11 we have neglected the correlations between the random fluctuations in the activity of a cell and the particular realization of its output connectivity. This is justified since such correlations are weak in the large N limit. Similarly, the variance $\alpha_{\mathrm{k}}$ of the input is

$$\alpha_{\mathrm{k}}(\mathrm{t})=\left[\left(\delta \mathrm{u}_{\mathrm{k}}^{\mathrm{i}}(\mathrm{t})\right)^{2}\right]=\sum_{1, \mathrm{l}^{\prime}=1}^{2} \sum_{\mathrm{j}, \mathrm{j}^{\prime}=1}^{\mathrm{N}_{\mathrm{l}}}\left[\left(\delta\left(\mathrm{~J}_{\mathrm{kl}}^{\mathrm{ij}} \sigma_{\mathrm{l}}^{\mathrm{j}}(\mathrm{t})\right)\right)^{2}\right],$$

where $\delta \mathrm{X} \equiv \mathrm{X}-[\mathrm{X}]$. Observing that $\left[\left(\mathrm{J}_{\mathrm{kl}}^{\mathrm{ij}} \sigma_{\mathrm{l}}^{\mathrm{j}}(\mathrm{t})\right)^{2}\right]=\mathrm{J}_{\mathrm{kl}}^{2} \mathrm{~m}_{\mathrm{l}} / \mathrm{N}$, whereas $\left[\left(\mathrm{J}_{\mathrm{kl}}^{\mathrm{ij}} \sigma_{\mathrm{l}}^{\mathrm{j}}(\mathrm{t})\right)\right]^{2}=\mathrm{J}_{\mathrm{kl}}^{2} \mathrm{~m}_{\mathrm{l}}^{2} \mathrm{~K} / \mathrm{N}^{2}$, which is negligible, one obtains equation 3.6.

4 Population Rates in the Balanced State

In a balanced state, the temporal fluctuations in the inputs are of the same order as the distance between the mean input relative to threshold (even when K is large). To show this, we need to probe the network temporal properties. Here we study the necessary consequences of the balanced state on the behavior of the population rates. A necessary condition for abalanced state is that both the excitatory and the inhibitory populations do not fire at their maximum rate, or are completely silent, when we take the limit $\mathrm{K} \rightarrow \infty$. In other words we look for solutions with $0<\mathrm{m}_{\mathrm{k}}<1$ in the large K limit.

To have equilibrium rates with $\mathrm{m}_{\mathrm{k}} \neq 0,1$ in the large K limit, both $\mathrm{u}_{\mathrm{E}}$ and $u_{\mathrm{I}}$ have to be finite in this limit. This means that the r.h.s. of equations 3.8 through 3.9 vanish to leading order. This leads to the following equations:

$$\begin{aligned} & \mathrm{Em}_{0}+\mathrm{m}_{\mathrm{E}}-\mathrm{J}_{\mathrm{E}} \mathrm{~m}_{\mathrm{I}}=\mathrm{O}(1 / \sqrt{\mathrm{K}}) \\ & \mathrm{Im}_{0}+\mathrm{m}_{\mathrm{E}}-\mathrm{J}_{\mathrm{I}} \mathrm{~m}_{\mathrm{I}}=\mathrm{O}(1 / \sqrt{\mathrm{K}}) \end{aligned}$$

Thus, in the large K limit we obtain

$$\begin{aligned} & \mathrm{m}_{\mathrm{E}}=\frac{\mathrm{J}_{\mathrm{I}} \mathrm{E}-\mathrm{J}_{\mathrm{E}} \mathrm{I}}{\mathrm{~J}_{\mathrm{E}}-\mathrm{J}_{\mathrm{I}}} \mathrm{~m}_{0} \equiv \mathrm{~A}_{\mathrm{E}} \mathrm{~m}_{0} \\ & \mathrm{~m}_{\mathrm{I}}=\frac{\mathrm{E}-\mathrm{I}}{\mathrm{~J}_{\mathrm{E}}-\mathrm{J}_{\mathrm{I}}} \mathrm{~m}_{0} \equiv \mathrm{~A}_{\mathrm{I}} \mathrm{~m}_{0} \end{aligned}$$

Since both $\mathrm{A}_{\mathrm{E}}$ and $\mathrm{A}_{\mathrm{I}}$ have to be positive, the coupling strengths have to satisfy

$$\frac{\mathrm{E}}{\mathrm{I}}>\frac{\mathrm{J}_{\mathrm{E}}}{\mathrm{~J}_{\mathrm{I}}}>1$$

or

$$\frac{\mathrm{E}}{\mathrm{I}}<\frac{\mathrm{J}_{\mathrm{E}}}{\mathrm{~J}_{\mathrm{I}}}<1 .$$

Besides this balanced solution, we should also examine the possibility of unbalanced solutions in which either $\mathrm{m}_{\mathrm{k}}=0$ and $\mathrm{u}_{\mathrm{k}}$ is of order $\sqrt{\mathrm{K}}$ and negative, or $\mathrm{m}_{\mathrm{k}}=1$ and $\mathrm{u}_{\mathrm{k}}$ is of order $\sqrt{\mathrm{K}}$ and positive. Equation 4.6 admits an unbalanced solution in which $\mathrm{m}_{\mathrm{E}}=0$. In this solution, $\mathrm{m}_{\mathrm{I}}$ is to leading order in k given by $\mathrm{m}_{\mathrm{I}}=\mathrm{Im}_{0} / \mathrm{J}_{\mathrm{I}}$ (since the leading order in $\mathrm{u}_{\mathrm{I}}$ should vanish) so that

$$\mathrm{u}_{\mathrm{E}}=\sqrt{\mathrm{K}}\left(\mathrm{E}-\mathrm{J}_{\mathrm{E}} \mathrm{I} / \mathrm{J}_{\mathrm{I}}\right) \mathrm{m}_{0}<0$$

Furthermore if $\mathrm{J}_{\mathrm{E}}<1$ and $\mathrm{J}_{\mathrm{I}}<1$, there exists a solution with $\mathrm{m}_{\mathrm{E}}=\mathrm{m}_{\mathrm{I}}=1$ even for $\mathrm{m}_{0}=0$. In this solution, $\mathrm{u}_{\mathrm{k}}$ satisfies to leading order

$$\mathrm{u}_{\mathrm{k}}=\sqrt{\mathrm{K}}\left(1-\mathrm{J}_{\mathrm{k}}\right) \quad(\mathrm{k}=\mathrm{E}, \mathrm{I})$$

so $\mathrm{u}_{\mathrm{k}}$ is of order $\sqrt{\mathrm{K}}$ and positive.

Thus if we require that there be no stationary solutions with $\mathrm{m}_{\mathrm{E}}=0,1$ or $\mathrm{m}_{\mathrm{I}}=0,1$ for small $\mathrm{m}_{0}$, the following constraints have to be satisfied:

$$\begin{aligned} & \frac{\mathrm{E}}{\mathrm{I}}>\frac{\mathrm{J}_{\mathrm{E}}}{\mathrm{~J}_{\mathrm{I}}}>1 . \\ & \mathrm{J}_{\mathrm{E}}>1 . \end{aligned}$$

It is straightforward to show that these constraints eliminate all possible unbalanced states.

Throughout the article, we will assume that equations 4.9 and 4.10 are satisfied and that $\mathrm{m}_{0}$ is small enough so that $\mathrm{A}_{\mathrm{k}} \mathrm{m}_{0}<1$. Equations 4.3 and 4.4 imply that the network activity rates grow linearly with the external rate, $\mathrm{m}_{\mathrm{k}}=\mathrm{A}_{\mathrm{k}} \mathrm{m}_{0}$, even though microscopic dynamics is highly nonlinear. This is because the network dynamically finds an operating point at which the net input in both populations is balanced. Thus, the linearity in the network rates reflects the linearity of the synaptic summation underlying our model. 4.1 The Net Input. Equations 4.3 and 4.4 determine the average rates of the populations, but they must be consistent also with the general equilibrium results of equation 3.7. According to equations 4.3 and 4.4, the leading $\mathrm{O}(\sqrt{\mathrm{K}})$ contributions to $\mathrm{u}_{\mathrm{k}}$ cancel each other. Thus, the net value of $\mathrm{u}_{\mathrm{k}}$ is determined by subleading contributions, such as corrections of order $1 / \sqrt{\mathrm{K}}$ to $\mathrm{m}_{\mathrm{k}}$. In fact, equations 3.7 should be viewed as equations that determine the net synaptic inputs $u_{k}$ given the mean activity rates $m_{k}$, equations 4.3 and 4.4. It is useful to denote by $h(m)$ the scaled input of $m$, defined as the solution of the equation

$$\mathrm{m}=\mathrm{H}(-\mathrm{h}) .$$

Thus, equation 3.7 reduces to

$$\mathrm{u}_{\mathrm{k}}=\sqrt{\alpha_{\mathrm{k}}} \mathrm{~h}\left(\mathrm{~m}_{\mathrm{k}}\right)$$

The activity of neurons in cortex is usually much less than the saturation rate. It is therefore useful to consider the limit where $\mathrm{m}_{0} \ll 1$. In this regime, $\mathrm{m}_{\mathrm{k}} \ll 1$, and we can use the approximation

$$\mathrm{H}(\mathrm{x}) \approx \frac{\exp \left(-\mathrm{x}^{2} / 2\right)}{\sqrt{2 \pi}|\mathrm{x}|}$$

to obtain

$$h(m) \approx-\sqrt{2|\log m|}$$

Substituting this result in equation 4.12 yields

$$\mathrm{u}_{\mathrm{k}} \approx-\sqrt{2 \alpha_{\mathrm{k}}\left|\log \mathrm{~m}_{\mathrm{k}}\right|}$$

This relation between $\mathrm{m}_{\mathrm{k}}$ and $\mathrm{u}_{\mathrm{k}}$ will be needed to calculate the rate distribution (see section 5.1). 4.2 Finite K Corrections. For finite K the residuals of order $1 / \sqrt{\mathrm{K}}$ in the rates are not negligible, so that equations 4.3 and 4.4 no longer hold exactly. For finite K , the equilibrium activities satisfy $\mathrm{m}_{\mathrm{k}}=\mathrm{F}_{\mathrm{k}}\left(\mathrm{m}_{\mathrm{E}}, \mathrm{m}_{\mathrm{I}}\right)$, with $\mathrm{F}_{\mathrm{k}}$ given by equation A.5. However, as long as

$$\mathrm{m}_{\mathrm{k}} \gg \mathrm{~K}^{-1}$$

the gaussian assumption of the input statistics is a good approximation; hence, equations 3.7 still hold. Thus, the leading finite K corrections can be incorporated by resorting to the full mean-field equations: equations 3.3 through 3.10. In particular, the finite K equations for the fixed point are

$$\begin{aligned} & \mathrm{Em}_{0}+\mathrm{m}_{\mathrm{E}}-\mathrm{J}_{\mathrm{E}} \mathrm{~m}_{\mathrm{I}}=\left(\theta_{\mathrm{E}}+\sqrt{\alpha_{\mathrm{E}}} \mathrm{~h}\left(\mathrm{~m}_{\mathrm{E}}\right)\right) / \sqrt{\mathrm{K}} \\ & \mathrm{Im}_{0}+\mathrm{m}_{\mathrm{E}}-\mathrm{J}_{\mathrm{I}} \mathrm{~m}_{\mathrm{I}}=\left(\theta_{\mathrm{I}}+\sqrt{\alpha_{\mathrm{I}}} \mathrm{~h}\left(\mathrm{~m}_{\mathrm{I}}\right)\right) / \sqrt{\mathrm{K}} \end{aligned}$$

As long as $m_{0}$ is not small, the r.h.s. of these equations are small for large K ; hence the corrections to the linear solution, equations 4.3 and 4.4, are small. When $\mathrm{m}_{0}$ becomes sufficiently small (i.e., of order $1 / \sqrt{\mathrm{K}}$ or less), the strong nonlinearity in the single neuron dynamics reveals itself in a strong nonlinearity in the population response. In particular, the effect of the single neuron threshold $\theta_{\mathrm{k}}$ becomes important. This is seen in Figure 3, where the population rates are evaluated by the finite K equations-equations 4.17 and 4.18 -with $\mathrm{K}=1000$. For comparison, we also show the straight lines predicted by the large K limit. Since the steady-state rates in cortical networks are usually low, it is sometimes useful to incorporate the leading finite K corrections. Whenever we refer in subsequent figures to explicit values for K , we use equations 4.17 and 4.18 , unless otherwise stated. Except for thresholding the population rates, the finite K corrections affect only the quantitative results, not the qualitative predictions of the simple large K theory.

5 Spatial and Temporal Variability

So far we have been concerned only with the population average rates $\mathrm{m}_{\mathrm{k}}$. However, the fact that the population averages are not saturated does not necessarily imply that the system's state exhibits strong temporal variations.

Figure 3: The mean activity of the excitatory population (thick solid line) and the inhibitory population (thick dashed line) as a function of the input activity. For a network in which cells receive input from, on average, 1000 cells in each population. Forcomparison, the activities in the largeK limitarealso shown (thin solid line for the excitatory and thin dashed line for the inhibitory population, respectively). Parameter values: $\mathrm{E}=1, \mathrm{I}=0.8, \mathrm{~J}_{\mathrm{E}}=2, \mathrm{~J}_{\mathrm{I}}=1.8, \theta_{\mathrm{E}}=1$, and $\theta_{\mathrm{I}}=0.7$.

Specifically, a population average excitatory rate, $\mathrm{m}_{\mathrm{k}}$, may be the outcome of a fluctuating state where all the cells in the kth population fire a fraction $\mathrm{m}_{\mathrm{k}}$ of the time. However, it can also be achieved by a frozen state in which a fraction $\mathrm{m}_{\mathrm{k}}$ of the cells fire every time these cells are updated, while all other cells never fire. In other words, the population average does not distinguish between temporal and spatial fluctuations of activity levels.

Fortunately, the mean-field theory fully characterizes the statistics of both the spatial and the temporal fluctuations in the activities in the balanced state. The statistics can be expressed by writing the instantaneous activity of a cell as a threshold function of two random variables $x_{i}$ and $y_{i}(t)$,

$$\sigma_{\mathrm{k}}^{\mathrm{i}}(\mathrm{t})=\Theta\left(-\mathrm{u}_{\mathrm{k}}+\sqrt{\beta_{\mathrm{k}}} \mathrm{x}_{\mathrm{i}}+\sqrt{\alpha_{\mathrm{k}}-\beta_{\mathrm{k}}} \mathrm{y}_{\mathrm{i}}(\mathrm{t})\right) .$$

The means $\mathrm{u}_{\mathrm{k}}$ are given by equations 3.8 and 3.9. The parameter $\beta_{\mathrm{k}}$ is given by

$$\beta_{\mathrm{k}}=\mathrm{q}_{\mathrm{E}}+\mathrm{J}_{\mathrm{k}}^{2} \mathrm{q}_{\mathrm{I}}$$

The order parameter $q_{\mathrm{k}}$ is defined as

$$\mathrm{q}_{\mathrm{k}}=\frac{1}{\mathrm{~N}_{\mathrm{k}}} \sum_{\mathrm{i}=1}^{\mathrm{N}_{\mathrm{k}}}\left(\mathrm{~m}_{\mathrm{k}}^{\mathrm{i}}\right)^{2}$$

where $\mathrm{m}_{\mathrm{k}}^{\mathrm{i}}$ is the time-averaged activity rate of the ith cell,

$$\mathrm{m}_{\mathrm{k}}^{\mathrm{i}} \equiv\left\langle\sigma_{\mathrm{k}}^{\mathrm{i}}(\mathrm{t})\right\rangle$$

The symbol $\langle\ldots\rangle$ denotes average over long time. Both $\mathrm{x}_{\mathrm{i}}$ and $\mathrm{y}_{\mathrm{i}}(\mathrm{t})$ are independent gaussian variables with zero mean and unit variance.

Quenched fluctuations of synaptic inputs. The term proportional to $\mathrm{x}_{\mathrm{i}}$ represents a quenched random component of the synaptic input received by different cells and thus represents a spatial inhomogeneity in the rates. The origin of this inhomogeneity is twofold. Since the connectivity is random in our model, cells may differ in the number of synaptic inputs they have. This component is given by

$$\delta_{1}\left\langle\mathrm{u}_{\mathrm{k}}^{\mathrm{i}}\right\rangle=\sum_{\mathrm{l}=1}^{2} \sum_{\mathrm{j}=1}^{\mathrm{N}_{\mathrm{l}}} \delta \mathrm{~J}_{\mathrm{kl}}^{\mathrm{ij}}\left[\mathrm{~m}_{\mathrm{l}}^{\mathrm{j}}\right]$$

Here, $\delta \mathrm{J}_{\mathrm{kl}}^{\mathrm{ij}}=\mathrm{J}_{\mathrm{kl}}^{\mathrm{ij}}-\left[\mathrm{J}_{\mathrm{kl}}^{\mathrm{ij}}\right]$. In addition, different neurons are connected to different cells so that even if all the cells would have received the same number of inputs, the system would evolve into a state with a self-consistently developed spatial inhomogeneity. The second component can be written as

$$\delta_{2}\left\langle\mathrm{u}_{\mathrm{k}}^{\mathrm{i}}\right\rangle=\sum_{\mathrm{l}=1}^{2} \sum_{\mathrm{j}=1}^{\mathrm{N}_{\mathrm{l}}} \mathrm{~J}_{\mathrm{kl}}^{\mathrm{ij}} \delta \mathrm{~m}_{\mathrm{l}}^{\mathrm{j}}$$

where $\delta \mathrm{m}_{\mathrm{l}}^{\mathrm{i}} \equiv \mathrm{m}_{\mathrm{l}}^{\mathrm{i}}-\mathrm{m}_{\mathrm{l}}$. Adding the two contributions yields

$$\left[\left(\delta\left\langle\mathrm{u}_{\mathrm{k}}^{\mathrm{i}}(\mathrm{t})\right\rangle\right)^{2}\right]=\sum_{\mathrm{l}=1}^{2} \mathrm{~J}_{\mathrm{kl}}^{2} q_{\mathrm{I}}=\mathrm{q}_{\mathrm{E}}+\mathrm{J}_{\mathrm{k}}^{2} \mathrm{q}_{\mathrm{I}}=\beta_{\mathrm{k}}$$

Thus, this variance represents the fluctuation in both the number and the identity of input cells to the different cells.

Temporal fluctuations of synaptic inputs. The term in equation 5.1 that is proportional to $y_{i}(t)$ represents the stochastic component of the inputs to a cell-a temporally fluctuating component with a short-time correlation. This can be written as

$$\mathrm{u}_{\mathrm{k}}^{\mathrm{i}}(\mathrm{t})-\left\langle\mathrm{u}_{\mathrm{k}}^{\mathrm{i}}\right\rangle=\sum_{\mathrm{l}=1}^{2} \sum_{\mathrm{j}=1}^{\mathrm{N}_{\mathrm{l}}} \mathrm{~J}_{\mathrm{kl}}^{\mathrm{ij}}\left(\sigma_{\mathrm{k}}^{\mathrm{j}}(\mathrm{t})-\mathrm{m}_{\mathrm{k}}^{\mathrm{j}}\right),$$

from which one obtains

$$\begin{aligned} {\left[\left(u_{\mathrm{k}}^{\mathrm{i}}(\mathrm{t})-\left\langle u_{\mathrm{k}}^{\mathrm{i}}\right\rangle^{2}\right]\right.} & =\sum_{\mathrm{l}=1}^{2} \mathrm{~J}_{\mathrm{kl}}^{2}\left(\mathrm{~m}_{\mathrm{l}}-\mathrm{q}_{\mathrm{l}}\right)=\mathrm{m}_{\mathrm{E}}-\mathrm{q}_{\mathrm{E}}+\mathrm{J}_{\mathrm{k}}^{2}\left(\mathrm{~m}_{\mathrm{I}}-\mathrm{q}_{\mathrm{I}}\right) \\ & =\alpha_{\mathrm{k}}-\beta_{\mathrm{k}} \end{aligned}$$

Note that the variance of the temporal fluctuations in the inputs depends on $m_{k}-q_{k}$, which measures the temporal variability of the state. 5.1 Distribution of Time-Averaged Rates. The distribution of rates in the kth population is defined as

$$\rho_{\mathrm{k}}(\mathrm{~m}) \equiv \mathrm{N}_{\mathrm{k}}^{-1} \sum_{\mathrm{i}=1}^{\mathrm{N}_{\mathrm{k}}} \delta\left(\mathrm{~m}-\mathrm{m}_{\mathrm{i}}^{\mathrm{k}}\right)$$

The statistics of the time-averaged local rates can be derived by averaging equation 5.1 over $y_{i}(t)$ (which is equivalent to average over time),

$$\mathrm{m}_{\mathrm{k}}^{\mathrm{i}}=\mathrm{m}_{\mathrm{k}}\left(\mathrm{x}_{\mathrm{i}}\right)=\mathrm{H}\left(\frac{-\mathrm{u}_{\mathrm{k}}+\sqrt{\beta_{\mathrm{k}}} \mathrm{x}_{\mathrm{i}}}{\sqrt{\alpha_{\mathrm{k}}-\beta_{\mathrm{k}}}}\right)$$

Thus, the distribution of $\mathrm{m}_{\mathrm{k}}^{\mathrm{i}}$ is fully determined by its first two moments. Averaging this equation over $\mathrm{x}_{\mathrm{i}}$ yields equation 3.7. Similarly, squaring equation 5.11 and averaging over $x_{i}$ yields

$$\mathrm{q}_{\mathrm{k}}=\int \mathrm{Dx}\left[\mathrm{H}\left(\frac{-\mathrm{u}_{\mathrm{k}}+\sqrt{\beta_{\mathrm{k}}} \mathrm{x}}{\sqrt{\alpha_{\mathrm{k}}-\beta_{\mathrm{k}}}}\right)\right]^{2}$$

Here we have used the gaussian measure, $\mathrm{Dx} \equiv \mathrm{dx} \exp \left(-\mathrm{x}^{2} / 2\right) / \sqrt{2 \pi}$. In general, $q_{k}$ satisfies $\left(m_{k}\right)^{2} \leq q_{k} \leq m_{k}$. The smaller $q_{k}$ is, the more homogeneous is the rate distribution. In a frozen state in which a fraction $\mathrm{m}_{\mathrm{k}}$ of the cells are active every time they are updated, while all other cells are always quiescent, $\mathrm{q}_{\mathrm{k}}$ is given by $\mathrm{q}_{\mathrm{k}}=\mathrm{m}_{\mathrm{k}}$. On the other hand, if all cells in the population have a probability $\mathrm{m}_{\mathrm{k}}$ of being active each time they are updated, $\mathrm{m}_{\mathrm{i}}^{\mathrm{k}}=\mathrm{m}_{\mathrm{k}}, \mathrm{q}_{\mathrm{k}}=\left(\mathrm{m}_{\mathrm{k}}\right)^{2}$. Equations 5.12 have two solutions: an unstable solution with $q_{k}=m_{k}$, corresponding to a frozen state, and a stable solution, $\left(\mathrm{m}_{\mathrm{k}}\right)^{2}<\mathrm{q}_{\mathrm{k}}<\mathrm{m}_{\mathrm{k}}$, which corresponds to a temporally fluctuating state. Although the frozen solution is unstable, its existence highlights the fact that the temporal variability in our system is purely of deterministic origin and is not induced by external stochastic sources.

Generalizing equation 5.12, we can write,

$$\rho_{\mathrm{k}}(\mathrm{~m})=\int \mathrm{D} \mathrm{x} \delta\left(\mathrm{~m}-\mathrm{m}_{\mathrm{k}}(\mathrm{x})\right)$$

In section A. 1 we analyze the properties of this distribution. A numerical evaluation of $\rho_{\mathrm{k}}(\mathrm{m})$ is shown in Figure4, which displays the rate distribution of the excitatory activity for different values of $\mathrm{m}_{\mathrm{E}}$. The distribution is plotted against $\mathrm{m} / \mathrm{m}_{\mathrm{E}}$. The synaptic couplings were kept constant, while the mean rates were varied by adjusting the external rate $\mathrm{m}_{0}$. For high mean activity levels, the distribution has a pronounced skewed shape. Note, however, that according to equation 5.11 the distribution of the time-averaged inputs $u_{i}^{k}$ to the cells is gaussian for all values of $m_{k}$.

In the low rate limit, $\mathrm{m}_{0} \ll 1$, equation 5.12 can be solved using equations 4.13 through 4.15 , yielding to leading order,

$$\mathrm{q}_{\mathrm{k}}=\mathrm{m}_{\mathrm{k}}^{2}+0\left(\mathrm{~m}_{\mathrm{k}}^{3}\left|\log \mathrm{~m}_{\mathrm{k}}\right|\right)$$

Thus, if the network evolves to a state with low average activity levels, $\mathrm{m}_{\mathrm{k}} \ll 1$, $\mathrm{q}_{\mathrm{k}}$ is slightly larger than $\mathrm{m}_{\mathrm{k}}^{2}$. The fact that $\mathrm{q}_{\mathrm{k}} \ll \mathrm{m}_{\mathrm{k}}$ implies that the balanced state is characterized by strong temporal fluctuations in the activity of the individual cells. On the other hand, the fact that $\mathrm{q}_{\mathrm{k}}$ is not exactly equal to $\mathrm{m}_{\mathrm{k}}^{2}$ reflects the spatial inhomogeneity in the time-averaged rates within a population. Equation 5.14 implies that when the mean activity $\mathrm{m}_{\mathrm{k}}$ decreases, the width of the distribution is proportional to $\left(\mathrm{m}_{\mathrm{k}}\right)^{3 / 2}$; it decreases faster than the mean $\mathrm{m}_{\mathrm{k}}$. Thus, for low mean activity, $\rho_{\mathrm{k}}(\mathrm{m})$ becomes narrowly peaked at $\mathrm{m}=\mathrm{m}_{\mathrm{k}}$, as shown in Figure 4. The reason for the narrow peak is that in our model, the fluctuations in the input are related to the fluctuations in the feedback from the network, hence their variance becomes small as the activity in the network decreases (see equation 5.2.) 5.2 Time-D elayed Autocorrelations. In order to complete the statistical characterization of the balanced state, we have to determine the statistics of the temporal fluctuations in the activities of single cells or, equivalently, the temporal fluctuations in their input. We have already stated that the temporal fluctuations in $\mathrm{u}_{\mathrm{i}}^{\mathrm{k}}(\mathrm{t})$ obey gaussian statistics, with variance given by $\alpha_{\mathrm{k}}-\beta_{\mathrm{k}}$. To characterize its statistics fully, we have to evaluate its autocorrelations. Using arguments similar to those already outlined, it is straightforward to show that the autocorrelation of the input is linearly related to

Figure 4: Distribution of the activity rates of the excitatory population for two different values of the average rate in the large K limit: $\mathrm{m}_{\mathrm{E}}=0.01$ (solid line) and $\mathrm{m}_{\mathrm{E}}=0.1$ (dashed line). The distributions are shown as a function of the local rate divided by the mean rate. Parameter values as in Figure 3.

the autocorrelations in the local activities,

$$\beta_{\mathrm{k}}(\tau)=\left[\left\langle\delta \mathrm{u}_{\mathrm{i}}^{\mathrm{k}}(\mathrm{t}) \delta \mathrm{u}_{\mathrm{i}}^{\mathrm{k}}(\mathrm{t}+\tau)\right\rangle\right]=\mathrm{q}_{\mathrm{E}}(\tau)+\mathrm{J}_{\mathrm{k}}^{2} \mathrm{q}_{\mathrm{I}}(\tau),$$

where $\mathrm{q}_{\mathrm{k}}(\tau)$ is the time-delayed autocorrelations of the local activities,

$$\mathrm{q}_{\mathrm{k}}(\tau)=\mathrm{N}_{\mathrm{k}}^{-1} \sum_{\mathrm{i}=1}^{\mathrm{N}_{\mathrm{k}}}\left\langle\sigma_{\mathrm{i}}^{\mathrm{k}}(\mathrm{t}) \sigma_{\mathrm{i}}^{\mathrm{k}}(\mathrm{t}+\tau)\right\rangle,$$

and as before $\langle\ldots\rangle$ denotes average over t . Note that $\mathrm{q}_{\mathrm{k}}(0)=\mathrm{m}_{\mathrm{k}}$, whereas $\mathrm{q}_{\mathrm{k}}(\tau \rightarrow \infty)=\mathrm{q}_{\mathrm{k}}$. Likewise, $\beta_{\mathrm{k}}(0)=\alpha_{\mathrm{k}}$, whereas $\beta_{\mathrm{k}}(\tau \rightarrow \infty)=\beta_{\mathrm{k}}$. Using this relation, the following self-consistent equation for $\mathrm{q}_{\mathrm{k}}(\tau)$ (with $\tau \geq 0$ ) is obtained:

$$\begin{aligned} \tau_{\mathrm{k}} \frac{\mathrm{dq}_{\mathrm{k}}(\tau)}{\mathrm{d} \tau}= & -\mathrm{q}_{\mathrm{k}}(\tau)+\int_{0}^{\infty} \frac{\mathrm{dt}}{\tau_{\mathrm{k}}} \\ & \times \exp \left(-\mathrm{t} / \tau_{\mathrm{k}}\right) \int \mathrm{Dx}\left[\mathrm{H}\left(\frac{-\mathrm{u}_{\mathrm{k}}-\sqrt{\beta_{\mathrm{k}}(\mathrm{t}+\tau)} \mathrm{x}}{\sqrt{\alpha_{\mathrm{k}}-\beta_{\mathrm{k}}(\mathrm{t}+\tau)}}\right)\right]^{2} . \end{aligned}$$

Figure 5: Population-averaged autocorrelation for the excitatory population in the large K limit (solid line). The dashed line shows the autocorrelation for a population of cells with the same rate distribution but Poissonian updating. Parameter values: $\tau=0.9, \mathrm{~m}_{0}=0.1$, and other parameters as in Figure 3.

Note that the integral overt in equation 5.17 results from averaging over the distribution of update time intervals. The solution of this integral equation yields a function $q_{k}(t)$, which decays to its equilibrium value with a time constant of the order of $\tau_{\mathrm{k}}$. A numerical solution of equation 5.17 for $q_{\mathrm{k}}(\tau)$ is shown in Figure 5. As can be seen, the autocorrelations are larger than those predicted by Poisson statistics. This enhancement of short-time correlations reflect the refractoriness in the activities of the cells that project the cell. 5.3 Numerical Realization of Synaptic Inputs to a Cell. In order to demonstrate the nature of the fluctuating synaptic inputs to a single excitatory cell in the balanced state, we have numerically generated samples of stochastic gaussian processes, which simulate the fluctuations of the synaptic inputs to a single excitatory cell. In order to show explicitly the effect of balancing, we have simulated separately the total excitatory and inhibitory components of $\mathrm{u}_{\mathrm{E}}^{\mathrm{i}}(\mathrm{t})$. The time average of the total excitatory (inhibitory) component is itself sampled from a gaussian distribution with a mean $\sqrt{\mathrm{K}}\left(\mathrm{m}_{\mathrm{E}}+\mathrm{Em}_{0}\right)\left(\sqrt{\mathrm{K}} \mathrm{J}_{\mathrm{E}} \mathrm{m}_{\mathrm{I}}\right)$ and a variance $\mathrm{q}_{\mathrm{E}}\left(\mathrm{J}_{\mathrm{E}} \mathrm{q}_{\mathrm{I}}\right)$. The time-dependent fluctuations of the total excitatory (inhibitory) inputhave atime-delayed au-

Figure 6: Temporal structure of the input to an excitatory cell. The upper panel shows the total excitatory input, consisting of the external input and the excitatory feedback (upper trace), the total inhibitory input (lower trace), as well as the net input (middle trace). They are calculated by sampling from the timecorrelated gaussian statistics predicted by the theory. The times when the cell switched to the active states are indicated. Parameter values: $\mathrm{m}_{0}=0.04$ and other parameters as in Figure 5. $\mathrm{K}=1000$ was used to calculate the average input.

tocorrelation equal to $\mathrm{q}_{\mathrm{E}}(\tau)-\mathrm{q}_{\mathrm{E}}\left(\mathrm{J}_{\mathrm{E}}^{2}\left(\mathrm{q}_{\mathrm{I}}(\tau)-\mathrm{q}_{\mathrm{I}}\right)\right)$. The results are shown in Figure 6, where we have used $\mathrm{K}=1000$. They demonstrate that the total excitatory and inhibitory inputs are large compared to the threshold and have fluctuations that are small compared to their mean. Because the network is in the balanced state, the net input is of the same order as the threshold, and the fluctuations bring the input above threshold at irregular intervals. In the lower part of the figure, we show the output state of the cell. This is evaluated by generating the sequence of update times and thresholding the net input at these times. Because of the update rule, the cell does not switch from passive to active every time the net input crosses the threshold.

In Figure 7 we present the ISI histogram of the cell. Because the interval between spikes is a convolution of two random events-first, a transition from 1 to 0, and then a transition from 0 to 1-the ISI vanishes at small intervals. Thus, the above definition of a spike to some extent captures the refractoriness of real spikes. In fact, if we ignore the short time correlations in the activities, the ISI of the ith (say, excitatory) cell with an average rate

Figure 7: Interspike interval distribution for the cell shown in Figure 6 (solid line). The distribution was determined by measuring the time between consecutive switches from the inactive to the active state until $5 \times 10^{5}$ intervals had been accumulated. The dashed line shows the interspike interval distribution for Poissonian updating with $\mathrm{m}_{\mathrm{i}}=0.06$.

$\mathrm{m}_{\mathrm{i}}$ can be shown to be simply

$$\mathrm{I}_{\mathrm{i}}(\mathrm{t})=\frac{\mathrm{m}_{\mathrm{i}}\left(1-\mathrm{m}_{\mathrm{i}}\right)}{\tau_{\mathrm{E}}\left(1-2 \mathrm{~m}_{\mathrm{i}}\right)}\left(\exp \left(-\mathrm{m}_{\mathrm{i}} \mathrm{t} / \tau_{\mathrm{E}}\right)-\exp \left(-\left(1-\mathrm{m}_{\mathrm{i}}\right) \mathrm{t} / \tau_{\mathrm{E}}\right)\right), \mathrm{t} \geq 0$$

This function rises linearly from zero and peaks at $\mathrm{t} \propto \tau_{\mathrm{E}}$. For intervals of the order of $\tau_{\mathrm{E}} / \mathrm{m}_{\mathrm{i}}$ or longer, $\mathrm{I}(\mathrm{t})$ decays purely exponentially with a decay constant $\mathrm{m}_{\mathrm{i}} / \tau_{\mathrm{E}}$ as in the ISI of a single Poisson process with a rate $\mathrm{m}_{\mathrm{i}} / \tau_{\mathrm{E}}$. Comparison with Figure 7 shows that this is indeed a very good approximation of the ISI of our model.

Finally, it should be noted that because of the sparsity of the connectivity, different cells receive input from different subpopulations, so that the fluctuations in their input will be only very weakly correlated. As a result, the correlations in their activity will be very small.

6 Stability of the Balanced State

To determine the stability of the balanced state, we have to study the response of the system to small perturbations in the population activity rates. However, because of the nature of the balanced state, we have to distinguish two scales of perturbations: local perturbations, in which the deviations in the rates are small compared to $1 / \sqrt{\mathrm{K}}$, and global perturbations, in which these deviations are large compared to $1 / \sqrt{\mathrm{K}}$. 6.1 Local Stability. Local stability of the balanced state requires that a sufficiently small perturbation in the populations rates will decay to zero. In our case, a sufficiently small perturbation means that it initially causes only a small disruption of the balanced state. This means that the perturbations are small not only compared to $\mathrm{m}_{\mathrm{k}}$ but also compared to $1 / \sqrt{\mathrm{K}}$, so that the perturbations of the inputs to the cells are initially small. We therefore considera solution of equations 3.3 with an initial condition $\mathrm{m}_{\mathrm{k}}(0)=\mathrm{m}_{\mathrm{k}}+\delta \mathrm{m}_{\mathrm{k}}(0)$ with a small $\delta \mathrm{m}_{\mathrm{k}}(0)$ where $\left|\delta \mathrm{m}_{\mathrm{k}}(0)\right| \ll 1 / \sqrt{\mathrm{K}}$. In this case, the perturbation of the total mean input $\mathrm{u}_{\mathrm{k}}$ is also small; hence we can linearize the dynamic equations 3.3 around their fixed point. Thus, $\delta \mathrm{m}_{\mathrm{k}}(\mathrm{t})=\mathrm{m}_{\mathrm{k}}(\mathrm{t})-\mathrm{m}_{\mathrm{k}}$ satisfy a linear equation of the form

$$\tau_{\mathrm{k}} \frac{\mathrm{~d}}{\mathrm{dt}} \delta \mathrm{~m}_{\mathrm{k}}(\mathrm{t})=-\delta \mathrm{m}_{\mathrm{k}}(\mathrm{t})+\sqrt{\mathrm{K}} \sum_{\mathrm{l}=1,2} \mathrm{f}_{\mathrm{kl}} \delta \mathrm{~m}_{\mathrm{l}}(\mathrm{t}) .$$

Calculating $\mathrm{f}_{\mathrm{kl}}$ by partial differentiation of the r.h.s. of equation 3.3 yields

$$\mathrm{f}_{\mathrm{kl}}=\frac{\exp \left(-\mathrm{u}_{\mathrm{k}}^{2} / 2 \alpha_{\mathrm{k}}\right) \mathrm{J}_{\mathrm{kl}}}{\sqrt{2 \pi \alpha_{\mathrm{k}}}} .$$

Solving equations 6.1 one obtains $\delta \mathrm{m}_{\mathrm{k}}(\mathrm{t})=\delta \mathrm{m}_{\mathrm{k}, 1} \exp \left(\lambda_{1} \mathrm{t}\right)+\delta \mathrm{m}_{\mathrm{k}, 2} \exp \left(\lambda_{2} \mathrm{t}\right)$ where the eigenvalues $\lambda_{1}$ and $\lambda_{2}$ of the 2 by 2 equations (see equations 6.1) are both of order $\sqrt{\mathrm{K}}$. Requiring that their real part be negative yields a condition on $\tau$ of the form

$$\tau<\tau_{\mathrm{L}},$$

where $\tau_{\mathrm{L}}$ is of order 1; its precise value depends on the system parameters. Since both $\lambda_{1}$ and $\lambda_{2}$ are of order $\sqrt{\mathrm{K}}$, if $\tau<\tau_{\mathrm{L}}$, small perturbations will decay with an extremely short time constant of order $1 / \sqrt{\mathrm{K}}$. This is due to the strong negative feedback, of order $\sqrt{\mathrm{K}}$, generated by the strong synaptic couplings. 6.2 Global Stability. The local stability condition in equation 6.3 guarantees that a perturbation smaller than $\mathrm{O}(1 / \sqrt{\mathrm{K}})$ will die out. It is therefore important to ask whether the balanced state is stable also to perturbations that are large compared to this order. However, such perturbations will generate a large disruption in the inputs $\mathrm{u}_{\mathrm{k}}$, of order $\sqrt{\mathrm{K}}$; hence, linearization of the dynamic equations is inadequate. We therefore have to consider the nonlinear evolution of perturbations in the rates under equations 3.3. In fact, since the perturbation destroys the balance between excitation and inhibition, $\mathrm{H}\left(-\mathrm{u}_{\mathrm{k}} / \sqrt{\alpha_{\mathrm{k}}}\right)$ of equation 3.3 can be approximated by $\Theta\left(\mathrm{u}_{\mathrm{k}}\right)$; hence the evolution of the perturbations is described by

$$\tau_{\mathrm{k}} \frac{\mathrm{~d}}{\mathrm{dt}} \delta \mathrm{~m}_{\mathrm{k}}(\mathrm{t})=-\delta \mathrm{m}_{\mathrm{k}}(\mathrm{t})+\Theta\left(\delta \mathrm{m}_{\mathrm{E}}-\mathrm{J}_{\mathrm{k}} \delta \mathrm{~m}_{\mathrm{I}}\right)-\mathrm{m}_{\mathrm{k}}$$

These equations are piecewise linear and therefore can be solved explicitly. One finds that the solution of these equations decays to zero provided that the inhibitory time constant satisfies

$$\tau<\tau_{\mathrm{G}},$$

where

$$\tau_{\mathrm{G}}=\mathrm{J}_{\mathrm{E}} \min \left\{\sqrt{\frac{\mathrm{~J}_{\mathrm{I}} \mathrm{~m}_{\mathrm{I}}\left(1-\mathrm{m}_{\mathrm{I}}\right)}{\mathrm{J}_{\mathrm{E}} \mathrm{~m}_{\mathrm{E}}\left(1-\mathrm{m}_{\mathrm{E}}\right)}}, \frac{\mathrm{m}_{\mathrm{I}}}{\mathrm{~m}_{\mathrm{E}}}, \frac{1-\mathrm{m}_{\mathrm{I}}}{1-\mathrm{m}_{\mathrm{E}}}\right\} .$$

In conclusion, the global stability condition guarantees that starting from arbitrary initial values $\mathrm{m}_{\mathrm{k}}(0)$, the population rates eventually will approach the balanced regime characterized by local fields $\mathrm{u}_{\mathrm{k}}$, which are of order 1 and not $\sqrt{\mathrm{K}}$. In other words, the rates will deviate from the values of the balanced fixed point by at most $\mathrm{O}(1 / \sqrt{\mathrm{K}})$ quantities. Whether they will actually approach this fixed point or will converge to a limit cycle around it depends on the local stability condition-equation 6.3. Depending on the system parameters $\tau_{\mathrm{G}}$ may be greater or smaller than $\tau_{\mathrm{L}}$. Figure 8 shows the evolution of $\mathrm{m}_{\mathrm{E}}$ and $\mathrm{m}_{\mathrm{I}}$ in a network with $\mathrm{K}=1000$, for $\tau=1.3$, when the network starts far away from the balanced state. The initial evolution is similar to the global dynamics. It converges to the neighborhood of the balanced fixed point in an oscillatory manner characteristic of the dynamics of equation 6.4. In Figure 8B, we show the late portion of the dynamics, which corresponds to the local dynamics (see equations 6.1). For the parameters used in this figure, the large K critical inhibitory time constants are $\tau_{\mathrm{L}}=1.61$ and $\tau_{\mathrm{G}}=1.50$.

To illustrate the region of stability of the balanced state, we have calculated the phase diagram of the network in terms of two parameters: (1) the inhibitory time constant $\tau$ and (2) the ratio between the external input into the inhibitory population and the external input into the excitatory population. We have chosen to scale $\mathrm{m}_{0}$ so that the excitatory population rate is held fixed. The results areshown in Figure9, where both the local and global

Figure 8: Evolution to the stable fixed point. The average inhibitory rate is plotted against the average excitatory rate. (A) The evolution of the rates when the rates are initialized far from their steady-state values. (B) A close-up view of the approach to the fixed point. Parameters: $\tau=1.3$. The other parameter values as in Figure 3.

stability lines are presented. For these parameters, $\tau_{\mathrm{L}}$ is always smaller than $\tau_{\mathrm{G}}$. 6.3 Regimes of Instability. Stability of the balanced state requires that $\tau$ be smaller than both $\tau_{\mathrm{L}}$ and $\tau_{\mathrm{G}}$. It is of interest to consider what happens if this condition is not fulfilled.

1. Unbalanced limit cycle. $\tau>\max \left\{\tau_{\mathrm{L}}, \tau_{\mathrm{G}}\right\}$. In this case equations 3.3 possess a stable unbalanced limit cycle-a stable oscillatory solution with $\mathrm{u}_{\mathrm{k}}(\mathrm{t})$ of order $\sqrt{\mathrm{K}}$. This is shown in Figure 10A.
2. Balanced limit cycle. $\tau_{\mathrm{L}}<\tau<\tau_{\mathrm{G}}$. In this case, perturbations that are of order 1 will decrease until they are of order $1 / \sqrt{\mathrm{K}}$, while perturbations that are small compared to $1 / \sqrt{\mathrm{K}}$ will increase until they are of order $1 / \sqrt{\mathrm{K}}$. Since there are no fixed points with $\delta \mathrm{m}_{\mathrm{k}}$ of order $1 / \sqrt{\mathrm{K}}$, this means that there has to be a stable limit cycle with an amplitude of order $1 / \sqrt{\mathrm{K}}$. Thus, in this regime, the system converges to a limit cycle that maintains the approximate balance between excitation and inhibition. This is described schematically in Figure 10B.
3. Balanced fixed point with shrinking basin: $\tau_{\mathrm{G}}<\tau<\tau_{\mathrm{L}}$. In this case, perturbations of order 1 go to a global limit cycle, while perturbations much smaller than $1 / \sqrt{\mathrm{K}}$ evolve to the fixed point. There must be an unstable limit cycle with amplitude of order $1 / \sqrt{\mathrm{K}}$ that separates perturbations that go to the global limit cycle and perturbations that go to the fixed point.

Figure 9: Critical time constants $\tau_{\mathrm{G}}$ (solid line) and $\tau_{\mathrm{L}}$ (dashed line) as a function of I/E. The external rate was adjusted to keep the excitatory activity level constant at $\mathrm{m}_{\mathrm{E}}=0.1 \mathrm{I} / \mathrm{E}$ was varied from 0 to $\left(\mathrm{J}_{\mathrm{I}}-\mathrm{m}_{\mathrm{E}}\right) /\left(\mathrm{J}_{\mathrm{E}}-\mathrm{m}_{\mathrm{E}}\right)$. For this range of $\mathrm{I} / \mathrm{E}, \mathrm{m}_{\mathrm{I}}$ varies from $\mathrm{m}_{\mathrm{E}} / \mathrm{J}_{\mathrm{I}}$ to 1 . Parameters: $\mathrm{E}=1, \mathrm{~J}_{\mathrm{E}}=2.0$, and $\mathrm{J}_{\mathrm{I}}=1.8$.

7 Inhomogeneous Thresholds

So far we have considered networks of identical neurons, except for their connectivity. Real neuronal systems exhibit a substantial inhomogeneity in single neuron properties. It is therefore important to consider how such inhomogeneities affect the behavior of our system. We will model the inhomogeneity by a variability in the thresholds of the neurons. Inhomogeneities in the local thresholds may have a particularly strong effect in a balanced state with low mean activity. The reason is that the intrinsic fluctuations are all generated by feedback from the network activity. Hence, they decrease in amplitude as the mean activity in the network drops. In particular, under these conditions, the intrinsic temporal fluctuations may not be of sufficiently large amplitude to overcome the quenched dispersion of local thresholds. Therefore, the important issue we address here is whether the balanced state remains temporally fluctuating in the limit of low mean activity in the presence of inhomogeneous thresholds, or whether it becomes a frozen state. We will show that the answer to these questions depends not only on the width of the threshold distribution but also on the form of its

Figure 10: Different scenarios when the fixed point is unstable. (A)A case where $\tau$ is larger than both $\tau_{\mathrm{G}}$ and $\tau_{\mathrm{L}}$. Here the rates evolve to the global limit cycle where the amplitude of the oscillations is of order 1. The solid line shows the evolution when the network is initiated outside the limit cycle; the dashed line corresponds to the trajectory for initial rates inside the limit cycle. Parameters as in Figure 8, except $\tau=1.8$. (B) Shows schematically a case where $\tau_{\mathrm{L}}<\tau<\tau_{\mathrm{G}}$. The network evolves to a limit cycle with the amplitude of order $\mathrm{K}^{-1 / 2}$. The figure shows the evolution of the rates with initial conditions far from the fixed point. The insert shows an expanded view of the area around the fixed point, with the trajectory of rates starting outside the limit cycle (solid line) and the trajectory of a network that was initiated close to the fixed point.

tail. We denote the local threshold of a neuron by $\theta_{\mathrm{i}}^{\mathrm{k}}+\theta_{\mathrm{k}}$, where $\theta_{\mathrm{k}}$ is the population-averaged threshold and $\theta_{\mathrm{i}}^{\mathrm{k}}$ is a quenched random variable with zero mean. We will call $\theta_{\mathrm{i}}^{\mathrm{k}}$ the local threshold. The mean activity rate of neurons in the kth population that have a local threshold $\theta$ is

$$\mathrm{m}_{\mathrm{k}}(\theta)=\mathrm{H}\left(\frac{\theta-\mathrm{u}_{\mathrm{k}}}{\sqrt{\alpha_{\mathrm{k}}}}\right),$$

and hence the population-averaged rate is

$$\mathrm{m}_{\mathrm{k}}=\int \mathrm{d} \theta \mathrm{P}(\theta) \mathrm{H}\left(\frac{\theta-\mathrm{u}_{\mathrm{k}}}{\sqrt{\alpha_{\mathrm{k}}}}\right),$$

where $\mathrm{P}(\theta)$ denotes the quenched distribution of $\theta$, and $\mathrm{u}_{\mathrm{k}}$ and $\alpha_{\mathrm{k}}$ are given as before by equations 3.8, 3.9, and 3.10. Note that we have absorbed the mean threshold, $\theta_{\mathrm{k}}$ in the definition of $\mathrm{u}_{\mathrm{k}}$ (see equations 3.8 and 3.9). 7.1 Distribution of Thresholds with Long Tails. We first consider a distribution with a long tail of low thresholds. A concrete example is

$$P(\theta)=\frac{1}{\Delta \sqrt{2 \pi}} \exp \left(-\frac{1}{2}(\theta / \Delta)^{2}\right) .$$

In this case, the spatial fluctuations in the inputs (relative to thresholds) consist of two gaussian terms. One is induced by the random connectivity and has a variance $\alpha_{\mathrm{k}}$, and the other is induced by the thresholds and has a variance $\Delta$. The balance conditions that determine the population rates (equations 4.3 and 4.4) still hold. In addition,

$$\mathrm{m}_{\mathrm{k}}=\mathrm{H}\left(\frac{-\mathrm{u}_{\mathrm{k}}}{\sqrt{\alpha_{\mathrm{k}}+\Delta^{2}}}\right)$$

which determines $\mathrm{u}_{\mathrm{k}}$, and

$$\mathrm{q}_{\mathrm{k}}=\int \mathrm{D} \mathrm{x}\left[\mathrm{H}\left(\frac{-\mathrm{u}_{\mathrm{k}}-\sqrt{\Delta+\beta_{\mathrm{k}}} \mathrm{x}}{\sqrt{\alpha_{\mathrm{k}}-\beta_{\mathrm{k}}}}\right)\right]^{2}$$

Now let us consider the limit of low mean rates, which is achieved by assuming that $\mathrm{m}_{0}$ is small. For fixed $\Delta$, if the mean rates become sufficiently low so that $\mathrm{m}_{\mathrm{k}} \ll \Delta$, the intrinsic variances $\alpha_{\mathrm{k}}$ and $\beta_{\mathrm{k}}$ can be neglected compared with $\Delta$; hence one obtains

$$\mathrm{m}_{\mathrm{k}} \approx \mathrm{q}_{\mathrm{k}} \approx \mathrm{H}\left(\frac{-\mathrm{u}_{\mathrm{k}}}{\Delta}\right)$$

The fact that $q_{k} \approx m_{k}$ implies that the state is essentially frozen, namely,

$$\mathrm{m}_{\mathrm{k}}(\theta) \approx \Theta\left(\mathrm{u}_{\mathrm{k}}-\theta\right)$$

and, consequently, the distribution of mean rates has a distinct bimodal shape,

$$\rho_{\mathrm{k}}(\mathrm{~m}) \approx\left(1-\mathrm{m}_{\mathrm{k}}\right) \delta(\mathrm{m})+\mathrm{m}_{\mathrm{k}} \delta(\mathrm{~m}-1), \quad \mathrm{m}_{\mathrm{k}} \ll 1,$$

as shown in Figure 11A. Thus, an unbounded threshold distribution has a relatively strong qualitative effect on the balance state, in the limit of low mean rate. 7.2 Bounded D istribution. We next consider the case of a bounded distribution of thresholds. As an example, we take a distribution of $\theta$ that is uniform between $-\Delta / 2$ and $+\Delta / 2$, and zero otherwise. In this case, equation 7.2 yields

$$\mathrm{m}_{\mathrm{k}}=\frac{1}{\Delta} \int_{-\Delta / 2}^{\Delta / 2} \mathrm{~d} \theta \mathrm{H}\left(\frac{-\mathrm{u}_{\mathrm{k}}+\theta}{\sqrt{\alpha_{\mathrm{k}}}}\right)$$

To assess the effect of $\Delta$, we analyze equation 7.9 in the low $\mathrm{m}_{\mathrm{k}}$ limit. In this case, the solution for $\mathrm{u}_{\mathrm{k}}$ is

$$\mathrm{u}_{\mathrm{k}}+\Delta / 2=\mathrm{O}\left(\sqrt{\mathrm{~m}_{\mathrm{k}}}\right) .$$

Figure 11: Distribution of the activities of the cells in the excitatory population in the large K limit. (A) Distribution for a network of neurons with a gaussian distribution of thresholds. The distribution is shown for population-averaged rates $\mathrm{m}_{\mathrm{E}}=0.01$ (solid line) and $\mathrm{m}_{\mathrm{E}}=0.1$ (dashed line). The insert shows the divergence $a t \mathrm{~m}=1$ of the distribution for $\mathrm{m}_{\mathrm{E}}=0.01$ with the density in arbitrary units. Parameter values $\Delta=0.2$ and other values as in Figure 3. (B) Distribution of activity levels of the cells in the excitatory population in the large K limit for a network of neurons with a bounded distribution of thresholds. The distribution is shown for mean rates $\mathrm{m}_{\mathrm{E}}=0.01$ (solid line) and $\mathrm{m}_{\mathrm{E}}=0.1$ (dashed line). Parameter values as in A. (C) Firing rate distribution for neurons in the right prefrontal cortex of a monkey attending to a complex stimulus (light source and sound) and executing a reaching movement. The rates were averaged over the duration of events that showed a significant response. The average rate was 15.8 Hz .

Thus, the population rates adjust themselves so that synapticinputisslightly below the smallest threshold in the population, $\theta_{\mathrm{k}}-\mathrm{D} / 2$; see equation 3.8. The small gap between the mean synaptic input and the minimal threshold is such that the temporal fluctuations of the network, with the low variance $\alpha_{\mathrm{k}}$, are sufficient to bring the neurons to threshold levels. Indeed, analyzing the rate distribution for this case, we find that it is unimodal with width $\sqrt{\mathrm{q}_{\mathrm{k}}}$, where

$$\mathrm{q}_{\mathrm{k}} \propto \Delta \alpha_{\mathrm{k}}^{3 / 2}$$

This means that the rate distribution is extremely broad and skewed. The full shape of the rate distribution is given by

$$\rho_{\mathrm{k}}(\mathrm{~m}) \approx \frac{\sqrt{\alpha_{\mathrm{k}} / 2 \Delta^{2}}}{\mathrm{~m} \sqrt{|\log \mathrm{~m}|}}, \mathrm{m}_{-}<\mathrm{m}<\mathrm{m}_{+}$$

and zero otherwise. The bounds of $m$ are:

$$m_{-} \propto \exp \left(-\Delta^{2} /\left(2 \alpha_{k}\right)\right)$$ $$\mathrm{m}_{+} \propto \Delta \sqrt{\alpha_{\mathrm{k}} /\left|\log \left(\alpha_{\mathrm{k}}\right)\right|} \gg \mathrm{m}_{\mathrm{k}}$$

The results (see equations 7.9-7.14) show that in the case of a bounded threshold distribution, the temporal variability remains strong even in the limit of low mean rates. However, the inhomogeneity strongly affects the shape of the rate distribution, making it more skewed and broader. Figure 11 B shows the results of numerical calculation of the rate distribution for the excitatory population, with a uniform distribution of thresholds between $-\Delta / 2$ and $\Delta / 2$, for different values of mean rates. Comparing Figures $4,11 \mathrm{~A}$, and 11 B , we see that for moderate mean rates $\mathrm{m}_{\mathrm{k}}=0.1, \Delta$ does not have a big effect on the shape of the distribution. However, when the network mean activity is lowered, the distribution peak shifts to values that are much smaller than the mean, while its tail extends to rates of the order of $\sqrt{\mathrm{m}_{\mathrm{E}}}$. In contrast, in the case of a homogeneous threshold, lowering the mean rates shifts the peak toward the mean and decreases the width of $\rho(\mathrm{m})$ (see Figure 4). In the case of a gaussian distribution, lowering the mean rates creates a pronounced bimodal distribution, characteristic of a frozen state, as seen in Figure 11A.

In general, for small $\mathrm{m}_{\mathrm{k}}$, a threshold distribution $\mathrm{P}(\theta)$ will yield a rate distribution $\rho_{\mathrm{k}}$ for population k that is given by

$$\rho_{\mathrm{k}}(\mathrm{~m})=\sqrt{2 \pi} \mathrm{P}\left(-\sqrt{\alpha_{\mathrm{k}}}\left(\mathrm{~h}(\mathrm{~m})+\tilde{\mathrm{h}}_{\mathrm{k}}\right)\right) \mathrm{e}^{\mathrm{h}^{2}(\mathrm{~m}) / 2}$$

where $\tilde{\mathrm{h}}_{\mathrm{k}} \equiv \mathrm{h}\left(\mathrm{m}_{\mathrm{k}}\right)$ is determined by

$$\int \mathrm{dm} \mathrm{~m} \rho_{\mathrm{k}}(\mathrm{~m})=\mathrm{m}_{\mathrm{k}}$$

If $\mathrm{P}(\theta)$ has tails that fall off as slow as or slower than a gaussian, $\rho_{\mathrm{k}}$ will diverge for $\mathrm{m}=0$ and $\mathrm{m}=1$; if $\mathrm{P}(\theta)$ falls off faster than a gaussian, $\rho_{\mathrm{k}}$ will be negligible for $\mathrm{m}<\mathrm{m}_{-}$and for $\mathrm{m}>\mathrm{m}_{+}$with, for small $\mathrm{m}_{\mathrm{k}}, \mathrm{m}_{-} \ll \mathrm{m}_{\mathrm{k}}$ and $\mathrm{m}_{\mathrm{k}} \ll \mathrm{m}_{+} \ll 1$. In this case, $\rho_{\mathrm{k}}$ can be approximated by

$$\rho_{\mathrm{k}}(\mathrm{~m}) \propto \frac{\mathrm{P}\left(-\sqrt{\alpha_{\mathrm{k}}}\left(\sqrt{2 \log (\mathrm{~m})}-\tilde{\mathrm{h}}_{\mathrm{k}}\right)\right)}{\mathrm{m} \sqrt{\log (\mathrm{~m})}}$$

form $-\mathrm{m}<\mathrm{m}_{+}$. Furthermore $\mathrm{P}\left(-\sqrt{\alpha_{\mathrm{k}}}\left(\sqrt{2 \log (\mathrm{~m})}-\tilde{\mathrm{h}}_{\mathrm{k}}\right)\right)$ varies only slowly with m for these rates.

Thus for a threshold distribution with a tail that falls off faster than a gaussian, the distribution of the rates goes to 0 for $\mathrm{m}=0$ and $\mathrm{m}=1$ and has a long power-law tail that extends up to a rate $\mathrm{m}_{+}$that is much larger than the average rate. In contrast, if the tails of the distribution fall off as slow as or slower than a gaussian, the rate distribution will peak at $\mathrm{m}=0$ and $\mathrm{m}=1$ if the average rate is sufficiently low. 7.3 Experimental Rate Distribution. The above results make a clear prediction about the shape of the rate distribution in a local population of neurons with low mean rates. It seems reasonable to compare these predictions with the distribution of rates in cortical neuronal pools of behaving animals. Figure 11C presents an experimentally determined rate histogram of neurons in the right prefrontal cortex of a monkey (Abeles, Bergman, & Vaadia, 1988). The data was taken from time intervals when the monkey was attending to a variety of stimuli (light sources and sound) or executing simple reaching movements. The average rate (of the neurons that showed any activity during the time of measurement) was 15.8 Hz . The observed histogram has a distinct unimodal skewed shape with a tail extending up to 80 Hz . These results are consistent with the theoretical predictions of Figure 11B.

8 Chaotic Nature of the Balanced State

The strong temporal fluctuations of the neuronal activity in our model and the resultant fast decay of temporal correlations strongly suggest that the balanced state corresponds to a chaotic attractor. However, to justify characterizing this state as chaotic, we need to study the sensitivity of the dynamic trajectory to small perturbations in the initial conditions. If the network evolves to a chaotic attractor, small perturbations in the state of the network should grow at least exponentially. After some time, the state of the network is far from the state the network would have been in had it not been perturbed. This definition of chaos is technically inapplicable to a system with discrete degrees of freedom such as ours, since in this case the size of a perturbation of the system state is bounded by the discreteness of the system's state. In our case, the minimum perturbation is changing the state of a single neuron. Nevertheless, in the limit of large network size, we can consider such a perturbation as infinitesimal, as described below.

We consider two copies of the network. In one copy, the states of the neurons are given by $\sigma_{1, k}^{i}(t)$; in the other, they are given by $\sigma_{2, k}^{i}(t)$. Both networks have the same connection matrices $\mathrm{J}_{\mathrm{kl}}^{\mathrm{ij}}$ and the same update schedule. The networks get the same constant input $\mathrm{m}_{0}(\mathrm{t})=\mathrm{m}_{0}$ and are assumed to have reached a balanced state with the same population rates,

$$\frac{1}{\mathrm{~N}_{\mathrm{k}}} \sum_{\mathrm{i}=1}^{\mathrm{N}_{\mathrm{k}}}\left\langle\sigma_{\mathrm{p}, \mathrm{k}}^{\mathrm{i}}(\mathrm{t})\right\rangle=\mathrm{m}_{\mathrm{k}} \quad \text { for } \mathrm{p}=1,2$$

The distance between the network states at time $t$ is defined as

$$\begin{aligned} \mathrm{D}_{\mathrm{k}}(\mathrm{t}) & =\frac{1}{\mathrm{~N}_{\mathrm{k}}} \sum_{\mathrm{i}=1}^{\mathrm{N}_{\mathrm{k}}}\left\langle\left(\sigma_{1, \mathrm{k}}^{\mathrm{i}}(\mathrm{t})-\sigma_{2, \mathrm{k}}^{\mathrm{i}}(\mathrm{t})\right)^{2}\right\rangle \\ & =\frac{1}{\mathrm{~N}_{\mathrm{k}}} \sum_{\mathrm{i}=1}^{\mathrm{N}_{\mathrm{k}}}\left\{\left\langle\sigma_{1, \mathrm{k}}^{\mathrm{i}}(\mathrm{t})\right\rangle+\left\langle\sigma_{2, \mathrm{k}}^{\mathrm{i}}(\mathrm{t})\right\rangle-2\left\langle\sigma_{1, \mathrm{k}}^{\mathrm{i}}(\mathrm{t}) \sigma_{2, \mathrm{k}}^{\mathrm{i}}(\mathrm{t})\right\rangle\right\} . \end{aligned}$$

Here the angular brackets do not mean average over time but average over all initial conditions of the two networks subject to the constraints that each individual network is at equilibrium (e.g., its $\mathrm{m}_{\mathrm{k}}$ and $\mathrm{q}_{\mathrm{k}}$ have the equilibrium values), and that the distance between the initial states of the two networks equals a given $\mathrm{D}_{\mathrm{k}}(0)$. If the network is in a chaotic state, the distance $\mathrm{D}_{\mathrm{k}}(\mathrm{t})$ of the cells in population k , defined by equation 8.2, should grow at least exponentially for small $\mathrm{D}_{\mathrm{k}}$, The maximum Lyapunov exponent $\lambda_{\mathrm{L}}$, defined by

$$\lambda_{\mathrm{L}} \equiv \lim _{\mathrm{D}_{\mathrm{k}} \rightarrow 0} \mathrm{D}_{\mathrm{k}}^{-1} \frac{\mathrm{dD}_{\mathrm{k}}}{\mathrm{dt}}$$

should be positive. Note that in calculating $\lambda_{\mathrm{L}}$, we will first take the large N limit of $\mathrm{D}_{\mathrm{k}}$ and then $\mathrm{D}_{\mathrm{k}} \rightarrow 0$ limit. To write the dynamics of $\mathrm{D}_{\mathrm{k}}$, it is useful to write $\mathrm{D}_{\mathrm{k}}(\mathrm{t})$ as

$$\mathrm{D}_{\mathrm{k}}(\mathrm{t})=2\left(\mathrm{~m}_{\mathrm{k}}-\mathrm{Q}_{\mathrm{k}}(\mathrm{t})\right)$$

where $\mathrm{Q}_{\mathrm{k}}(\mathrm{t})$ denotes the overlap of the two trajectories. In appendix A , we show that $\mathrm{Q}_{\mathrm{k}}(\mathrm{t})$ satisfies an equation similar to that of $\mathrm{q}_{\mathrm{k}}(\tau)$,

$$\begin{aligned} \tau_{\mathrm{k}} \frac{\mathrm{dQ}_{\mathrm{k}}}{\mathrm{dt}}= & -\mathrm{Q}_{\mathrm{k}} \\ & +\int \mathrm{Dx}\left[\mathrm{H}\left(\frac{-\mathrm{u}_{\mathrm{k}}+\sqrt{\gamma_{\mathrm{k}}(\mathrm{t})} \mathrm{x}}{\sqrt{\alpha_{\mathrm{k}}-\gamma_{\mathrm{k}}(\mathrm{t})}}\right)\right]^{2} \end{aligned}$$

with $\mathrm{u}_{\mathrm{k}}$ and $\alpha_{\mathrm{k}}$ are as above, and $\gamma_{\mathrm{k}}(\mathrm{t})$ given by

$$\gamma_{\mathrm{k}}(\mathrm{t})=\sum_{\mathrm{l}=1}^{2}\left(\mathrm{~J}_{\mathrm{kl}}\right)^{2} \mathrm{Q}_{\mathrm{l}}(\mathrm{t})=\mathrm{Q}_{\mathrm{E}}(\mathrm{t})+\mathrm{J}_{\mathrm{k}}^{2} \mathrm{Q}_{\mathrm{I}}(\mathrm{t}) .$$

This equation has two stationary solutions. One is

$$\mathrm{Q}_{\mathrm{k}}=\mathrm{m}_{\mathrm{k}},$$

which corresponds to a fully locked trajectories. This solution is unstable, as will be shown below. The stable fixed point is

$$\mathrm{Q}_{\mathrm{k}}=\mathrm{q}_{\mathrm{k}},$$

which corresponds to a fully desynchronized trajectory so that at long times, the correlations between the two trajectories at the same time are those induced by the time-independent average activities. Starting from any nonidentical states, the two trajectories eventually will desynchronize them completely. To find the initial rate of divergence, we expand equation 8.5 for small $\mathrm{D}_{\mathrm{k}}$ and find that to leading order, the distances satisfy

$$\tau_{\mathrm{k}} \frac{\mathrm{dD}_{\mathrm{k}}}{\mathrm{dt}}=\frac{2}{\pi} \frac{\mathrm{e}^{-\mathrm{u}_{\mathrm{k}}^{2} / 2 \alpha_{\mathrm{k}}}}{\sqrt{\alpha_{\mathrm{k}}}} \sqrt{\alpha_{\mathrm{k}}-\gamma_{\mathrm{k}}} .$$

Since $\alpha_{\mathrm{k}}-\gamma_{\mathrm{k}} \propto \mathrm{D}_{\mathrm{k}}$, equation 8.9 has a growing solution even if $\mathrm{D}_{\mathrm{k}}(0)=$ 0 . This implies that the Lyapunov exponent $\lambda_{\mathrm{L}}$ is infinitely large in the balanced state. Figure 12 shows the evolution of $\mathrm{D}_{\mathrm{E}} . \mathrm{D}_{\mathrm{E}}$ increases rapidly to the equilibrium value $\mathrm{D}_{\mathrm{E}}=2\left(\mathrm{~m}_{\mathrm{E}}-\mathrm{q}_{\mathrm{E}}\right)$, for arbitrarily small initial positive value. This should be contrasted with systems with finite positive Lyapunov exponents, where the initial rate of growth depends on the magnitude of the initial perturbation of the initial conditions. The divergence of $\lambda_{\mathrm{L}}$ in our system is related to the discreteness of the degrees of freedom, which implies an infinitely high microscopic gain: a small change in the inputs to a cell can cause a finite change in its state.

9 Tracking of Time-D ependent Input

We have shown that for a large range of parameters, anetwork with synaptic strengths of order $1 / \sqrt{\mathrm{K}}$ will evolve to a balanced state, and we investigated some of the characteristics of this state. But so far we have not addressed the question of potential functional advantages of this state. Why should a network generate an excitatory input that is much larger than the threshold input and then counterbalance this with a nearly equally large inhibitory input? If we consider the metabolic costs of such large currents, it seems

Figure 12: Evolution of the distance $\mathrm{D}_{\mathrm{E}}$ starting from a small initial distance in the large K limit. Parameters as in Figure 3 and $\tau=0.9$.

clear that a biological system would not choose such a mechanism unless it has some advantages over other mechanisms.

In this section we present one possible advantage of the balanced network. We have already shown that perturbations in the network rates, which are small compared to $1 / \sqrt{\mathrm{K}}$, die out in a time on the order of $1 / \sqrt{\mathrm{K}}$. Therefore, the network is very stable against small fluctuations in the rates. We now consider the consequences of this for the response of the system to time-dependent change in the external driving force $\mathrm{m}_{0}$.

If the external activity $\mathrm{m}_{0}$ changes suddenly by a small amount, on the order of $1 / \sqrt{\mathrm{K}}$, the equilibrium rates will change by an amount that is of the same order. So just after the change in external rate, the network rates differ slightly from the equilibrium rate. They will approach the new equilibrium rate on a time scale that is of the order $1 / \sqrt{\mathrm{K}}$, so the network rates adapt very fast to a sudden change in $\mathrm{m}_{0}$. This means that if $\mathrm{m}_{0}$ changes continuously with time, the network rates will track $\mathrm{m}_{0}$ very fast, provided that $\mathrm{m}_{0}$ does not change too rapidly. To quantify the speed of the tracking of a balanced network, we compare the network rates with the rates of a hypothetical network that tracks changes in the external rates instantaneously. In such a network, the rates $\mathrm{m}_{\mathrm{k}}^{\infty}$ satisfy $\mathrm{m}_{\mathrm{k}}^{\infty}(\mathrm{t})=\mathrm{m}_{\mathrm{k}}\left(\mathrm{m}_{0}(\mathrm{t})\right)$, where $\mathrm{m}_{\mathrm{k}}\left(\mathrm{m}_{0}\right)$ is the equilibrium rate for $m_{0}(t)=m_{0}$, which are given by

$$\mathrm{m}_{\mathrm{k}}^{\infty}(\mathrm{t})=\mathrm{H}\left(-\frac{\mathrm{u}_{\mathrm{k}}^{\infty}(\mathrm{t})}{\sqrt{\alpha_{\mathrm{k}}^{\infty}(\mathrm{t})}}\right),$$

with

$$\mathrm{u}_{\mathrm{k}}^{\infty}(\mathrm{t})=\sqrt{\mathrm{K}}\left(\mathrm{~J}_{\mathrm{k} 0} \mathrm{~m}_{0}(\mathrm{t})+\sum_{\mathrm{l}=1}^{2} \mathrm{~J}_{\mathrm{kl}} \mathrm{~m}_{\mathrm{l}}^{\infty}(\mathrm{t})\right)-\theta_{\mathrm{k}}$$

and

$$\alpha_{\mathrm{k}}^{\infty}(\mathrm{t})=\sum_{\mathrm{l}=1,2}\left(\mathrm{~J}_{\mathrm{kl}}\right)^{2} \mathrm{~m}_{\mathrm{l}}^{\infty}(\mathrm{t})$$

Note that to leading order in $\mathrm{K} \mathrm{m}_{\mathrm{k}}^{\infty}$ satisfies the balance condition

$$\mathrm{m}_{\mathrm{k}}^{\infty}(\mathrm{t})=\mathrm{A}_{\mathrm{k}} \mathrm{~m}_{0}(\mathrm{t})$$

However, equations 9.1 through 9.3 take into account also the $1 / \sqrt{\mathrm{K}}$ corrections in $\mathrm{m}_{\mathrm{k}}^{\infty}(\mathrm{t})$.

We now assume that

$$\mathrm{m}_{\mathrm{k}}(\mathrm{t})=\mathrm{m}_{\mathrm{k}}^{\infty}(\mathrm{t})+\mathrm{m}_{\mathrm{k}}^{1}(\mathrm{t}) / \sqrt{\mathrm{K}}$$

namely, that the deviation from perfect tracking of the instantaneous is only of order $1 / \sqrt{\mathrm{K}}$. The rates $\mathrm{m}_{\mathrm{k}}$ satisfy equation 3.3. To leading order in K this is

$$\tau_{\mathrm{k}} \frac{\mathrm{dm}_{\mathrm{k}}^{\infty}(\mathrm{t})}{\mathrm{dt}}=-\mathrm{m}_{\mathrm{k}}^{\infty}(\mathrm{t})+\mathrm{H}\left(-\frac{\mathrm{u}_{\mathrm{k}}^{\infty}(\mathrm{t})+\sum_{\mathrm{l}} \mathrm{~J}_{\mathrm{kl}} \mathrm{~m}_{\mathrm{l}}^{1}(\mathrm{t})}{\sqrt{\alpha_{\mathrm{k}}^{\infty}(\mathrm{t})}}\right), \mathrm{k}=1.2 .$$

Using equation 9.4 we obtain,

$$A_{k}\left(\tau_{k} \frac{\mathrm{dm}_{0}(\mathrm{t})}{\mathrm{dt}}+\mathrm{m}_{0}(\mathrm{t})\right)=\mathrm{H}\left(-\frac{\mathrm{u}_{\mathrm{k}}^{\infty}(\mathrm{t})+\sum_{\mathrm{l}} \mathrm{~J}_{\mathrm{kl}} \mathrm{~m}_{\mathrm{l}}^{1}(\mathrm{t})}{\sqrt{\alpha_{\mathrm{k}}^{\infty}(\mathrm{t})}}\right), \mathrm{k}=1,2$$

which determines the small deviations $\mathrm{m}_{\mathrm{k}}^{1}(\mathrm{t}) / \sqrt{\mathrm{K}}$ as functions of the timedependent drive $\mathrm{m}_{0}(\mathrm{t})$. Since $\mathrm{H}(\mathrm{x})$ is between 0 and 1 , equations 9.7 have a solution only for

$$0<\mathrm{m}_{0}+\tau_{\mathrm{k}} \frac{\mathrm{dm}_{0}}{\mathrm{dt}}<1 / \mathrm{A}_{\mathrm{k}}$$

This implies that the almost perfect tracking occurs for rates of change of the external input, which obeys the following bounds:

$$\max _{\mathrm{k}=1,2}-\frac{\mathrm{m}_{0}}{\tau_{\mathrm{k}}}<\frac{\mathrm{dm}_{0}}{\mathrm{dt}}<\min _{\mathrm{k}=1,2} \frac{1}{\tau_{\mathrm{k}}}\left(\frac{1}{\mathrm{~A}_{\mathrm{k}}}-\mathrm{m}_{0}\right) .$$

To understand these results qualitatively, let us consider a system in the balanced state with a fixed $\mathrm{m}_{0}$, where at timet $=\mathrm{t}_{0}, \mathrm{~m}_{0}$ is suddenly changed to $\mathrm{m}_{0}+\delta \mathrm{m}_{0}$. We assume that $\delta \mathrm{m}_{0}$ is much smaller than $\mathrm{m}_{0}$ but $\delta \mathrm{m}_{0} \sqrt{\mathrm{~K}}$ is of order 1. This is shown in Figure 13, where $\mathrm{m}_{0}$ is increased by a series of small steps. Because the input is $\sqrt{\mathrm{K}} \mathrm{m}_{0}(\mathrm{t})$, the small change in $\mathrm{m}_{0}$ initially causes a change of order 1 in the total input. Hence, the probability $\mathrm{P}_{\mathrm{k}}$ that the cells in the kth population, which are updated at time $t_{0}$, will go to the active state is initially increased by a large amount. This is denoted as $\Delta \mathrm{P}_{\mathrm{k}}$; $\Delta \mathrm{P}_{\mathrm{k}}$ is of order 1, as shown by the dashed curve in the figure. In fact, this probability is given by the r.h.s. of equation 9.4, which differs substantially from the previous equilibrium probability, $\mathrm{A}_{\mathrm{k}} \mathrm{m}_{0}$. This initial increase in the number of active cells causes a large inhibitory feedback, which causes $\mathrm{P}_{\mathrm{k}}$ to decrease quickly to its new equilibrium value, which is only slightly increased from its original equilibrium value, as seen in Figure 13. Thus, the initial response is highly nonlinear due to the initial disruption of the balance in the inputs and the highly nonlinear dynamics of single cells. This initial large response causes a fast rate of increase in the population rates, since $\delta \mathrm{m}_{\mathrm{k}} \approx \tau_{\mathrm{k}}^{-1} \mathrm{dt} \Delta \mathrm{P}_{\mathrm{k}}$, implying that $\delta \mathrm{m}_{\mathrm{k}}$ reaches the value $\mathrm{A}_{\mathrm{k}} \delta \mathrm{m}_{0}$ in time of order $\tau_{\mathrm{k}} / \delta \mathrm{m}_{0} \approx \tau_{\mathrm{k}} / \sqrt{\mathrm{K}}$; see the dotted line in Figure 13. The final change in the population rates linearly follows the change in the external input as required to maintain the balance between excitation and inhibition.

The limitation on the change in the external rate is readily explained by the maximum increase (decrease) in the network rate that the microscopic dynamics allows. The fastest the network rates can increase (decrease) is by putting all newly updated cells in the active (passive) state, that is, $\mathrm{P}_{\mathrm{k}}=1$ ( $\mathrm{P}_{\mathrm{k}}=0$ ), so that the change in the network rates is bounded by

$$-\mathrm{m}_{\mathrm{k}}<\tau_{\mathrm{k}} \frac{\mathrm{dm}_{\mathrm{k}}}{\mathrm{dt}}<1-\mathrm{m}_{\mathrm{k}}$$

If the external rates increase (decrease) faster than the bound, equations 9.8 the network will not stay in the balanced state during the rate change, so $\mathrm{u}_{\mathrm{k}}$ is of order $\sqrt{\mathrm{K}}$. Consequently, the input is above (below) the threshold for all cells of the kth population that are updated, and all updated cells are in the active (passive) state.

To compare the tracking capabilities of balanced networks with those of an unbalanced network, we consider a network of threshold linear neurons with synapses of strength $\mathrm{J}_{\mathrm{kl}} / \mathrm{K}$ for internetwork connections and $\tilde{\mathrm{J}}_{\mathrm{k} 0} / \mathrm{K}$ for the strengths of the synapses projecting from the external population and

Figure 13: Reaction of the excitatory population to input that is increased by small steps. The solid line shows the activity of a network that responds instantaneously. The dashed line shows the probability $\mathrm{P}_{\mathrm{E}}$ of updating to the active state for neurons that happen to update. The dotted line represents the populationaveraged activity of the population in a balanced network. Parameter values: $\mathrm{K}=1000, \tau=0.9$. Other parameters as in Figure 3.

the thresholds $\mathrm{T}_{\mathrm{k}}$ chosen so that the equilibrium rates of this network are the same as those for the balanced network. We choose the same neuronal time constants as in the balanced network. In this network, the rates satisfy

$$\tau_{\mathrm{k}} \frac{\mathrm{dm}_{\mathrm{k}}}{\mathrm{dt}}=-\mathrm{m}_{\mathrm{k}}+\left(\tilde{\mathrm{J}}_{\mathrm{k} 0} \mathrm{~m}_{0}+\mathrm{J}_{\mathrm{kE}} \mathrm{~m}_{\mathrm{E}}+\mathrm{J}_{\mathrm{kI}} \mathrm{~m}_{\mathrm{I}}-\mathrm{T}_{\mathrm{k}}\right)_{+},$$

with $(\mathrm{x})_{+}=(\mathrm{x}+|\mathrm{x}|) / 2$. If we set $\mathrm{m}_{\mathrm{k}}(\mathrm{t})=\mathrm{m}_{\mathrm{k}}^{\infty}\left(\mathrm{m}_{0}(\mathrm{t})\right)+\mathrm{m}_{\mathrm{k}}^{1}(\mathrm{t})$, the difference between the network rates and the rates of a perfectly tracking network, $\mathrm{m}_{\mathrm{k}}^{1}(\mathrm{t})$, satisfies

$$\begin{aligned} \tau_{\mathrm{E}} \frac{\mathrm{dm}_{\mathrm{E}}^{1}}{\mathrm{dt}} & =\left(\mathrm{J}_{\mathrm{EE}}-1\right) \mathrm{m}_{\mathrm{E}}^{1}+\mathrm{JeI}_{\mathrm{EI}} \mathrm{~m}_{\mathrm{I}}^{1}-\tau_{\mathrm{E}} \mathrm{~A}_{\mathrm{E}} \frac{\mathrm{dm}_{0}}{\mathrm{dt}} \\ \tau_{\mathrm{I}} \frac{\mathrm{dm}_{\mathrm{I}}^{1}}{\mathrm{dt}} & =\mathrm{J}_{\mathrm{IE}} \mathrm{~m}_{\mathrm{E}}^{1}+\left(\mathrm{JII}_{\mathrm{II}}-1\right) \mathrm{m}_{\mathrm{I}}^{1}-\tau_{\mathrm{I}} \mathrm{~A}_{\mathrm{I}} \frac{\mathrm{dm}_{0}}{\mathrm{dt}} \end{aligned}$$

In other words, $\mathrm{m}_{\mathrm{k}}^{1}$ will be of order 1 .

Figure 14: Population-averaged activity of the excitatory cells for an input that varies with time. The input is constant from time $t=0$ to $t=1$. Between $t=1$ and $\mathrm{t}=2$, the input increases linearly. Aftert $=2$, the input is again constant. The solid line shows the excitatory rate of a network that responds infinitely quickly; the dashed line shows the response of the balanced network. Also shown is the response of an unbalanced network of threshold linear neurons (dotted line). Parameters for the balanced network as in Figure 13. For the unbalanced network, see the text.

Thus, in the unbalanced network, the difference between the network rates and the rates in a perfectly tracking network will be of order $\sqrt{\mathrm{K}}$ times larger than in a balanced network.

Figures 14 and 15show a comparison of the tracking capabilities of a balanced network with $\mathrm{K}=1000$ and an unbalanced network with threshold linear units. Between $t=0$ and $t=1$, the networks are at equilibrium. In Figure 14 the external activity is ramped between $t=1$ and $t=2$,

$$\mathrm{m}_{0}(\mathrm{t})=\mathrm{m}_{0}+\mathrm{v}_{0} \mathrm{t},$$

and after $\mathrm{t}=2 \mathrm{~m}_{0}$ is kept constant again. The graph shows $\mathrm{m}_{\mathrm{E}}^{\infty}$ and $\mathrm{m}_{\mathrm{E}}$ for both networks plotted against time. Clearly the balanced network is much better than the unbalanced network in tracking the change in external rate. Similar results are seen in the case of a sinusoidal external input (see Figure 15). Finally, in Figure 16, we plot the rate of change of $\mathrm{m}_{\mathrm{E}}$ versus $v_{0}$ for the ramped input case. The results of this, as well as Figures 14 and

Figure 15: Average rate of the excitatory population for a sinusoidally varying input. The rates of the excitatory population in an instantaneously responding network (solid line), a balanced network (dashed line), and an unbalanced network (dotted line) are shown. Parameters as in Figure 13.

15 are based on a full finite K solution of the dynamics. We also show in Figure 16 the large K predictions according to which there is a sharp upper bound for fast tracking at a value of $\mathrm{v}_{0}$ given in equation 9.8.

10 Discussion

10.1 Asynchronous and Synchronized Chaos. The purpose of our the ory is to identify the different mechanisms by which the deterministic dynamics generates strongly irregularstates in large neural networks, in which each cell receives input from many other cells. To understand these mechanisms from a theoretical point of view, it is important to study the network behavior in the limits of large system size N and large connectivity index K . In a finite network with fluctuating dynamics, there will always be some degree of synchrony and some compensation between inhibition and excitation. It is thus impossible to single out balancing between excitation and inhibition as a mechanism for variability separate from synchronized chaos (Bush & Douglas 1991; Hansel & Sompolinsky, 1992, 1996). It is only in the limit of large N , where states with synchrony that does not vanish

Figure 16: Rate $\mathrm{v}_{\mathrm{E}}$ with which the average excitatory rate changes as a function of the rate of change $\mathrm{v}_{0}$ of the external input. The solid line shows $\mathrm{v}_{\mathrm{E}}$ for a network with $\mathrm{K}=1000$. The dashed line shows the same for the large K limit. Parameters as in Figure 13.

in this limit can be distinguished from states where the synchrony does vanish, that the different mechanisms become clearly separate. Similarly, the importance of the limit of large K is that for fixed finite K , network parameters may be tuned so that fluctuations in individual synaptic inputs generate fluctuations in the membrane potential of the postsynaptic cells. These fluctuations can be due to stochastic synaptic failures or variability in the presynaptic cells from within or outside the network. In other words, for a finite fixed K , the issue of balancing between excitation and inhibition is a quantitative one. Only in the large K limit can the distance between the net input and the threshold clearly be separated, in the balanced state, from the corresponding distance for the excitatory and inhibitory components.

The outcome of the present theory combined with our previous studies is that chaotic states in large, highly connected networks can be classified as synchronized chaos and asynchronous chaos. Synchronized chaos is likely to occur in fully connected networks, where K is proportional to N , yielding astrong overlap between inputs to different neurons. In this case, the chaotic state is characterized by cross-correlations between neuronal pairs whose amplitude is of order 1 even in the limit of $\mathrm{N} \rightarrow \infty$, thereby creating strong fluctuations of the common feedback. Thus, synchronized chaos can be viewed as resulting from an instability in the dynamics of the macroscopic degrees of freedom that comprise the common fluctuating mean field.

Asynchronous chaoticstates are distinguished by weak cross-correlations. In the present case, this is due to the sparseness of the connections. More specifically, in our networks, the amplitude of the cross-correlations has a broad distribution in the network due to inhomogeneity in the connectivity. Most of the cross-correlations are of the order $1 / \mathrm{N}$, where N is the network size. The maximal value of the cross-correlations occurs for pairs that are directly connected, and this cross-correlation is of the order of the strength of the synapse, $O(1 / \sqrt{K})$. Thus, chaos in this state is the result of instability in local degrees of freedom, similar to chaos in asymmetric spin glasses and neural networks. 10.2 Balanced State with Strong or Weak Synapses. The scaling of connection strength in our theory of the balanced state is different from conventional mean-field theories of highly connected networks. Most mean-field theories of large, highly connected neural networks assume that each connection is scaled as the inverse of the mean number of inputs to a neuron, K . In contrast, we scale the connections as $1 / \sqrt{\mathrm{K}}$. This aspect, together with the relative sparseness of the connections and the asynchrony of the dynamics, yields a highly irregular dynamical state, despite the fact that the singleneuron dynamics in our model is the simple threshold updates of binary units. The presence of relatively large connections is again analogous to the scaling of connections in highly connected spin glasses and random neural networks, where the interactions have to scale as the inverse square root of the connectivity index (Derrida, Gardner, & Zippelius, 1987; Sompolinsky, Crisanti, & Sommers, 1988). In Sompolinsky et al. (1988) the network is a fully connected asymmetric analog circuit with connections that are independent random variables with zero mean. The connections possess a square-root scaling with the number of inputs, as is natural for mean-field spin glasses (Mezard, Parisi, & Virasoro, 1987). In Derrida et al. (1987), the connectivity is randomly sparse, as in our model. The connections store random memories so that in the limit of a large K (and correspondingly large number of stored patterns), they are effectively random in sign and exhibit chaotic dynamics similar to the asymmetric spin glass. In contrast, in our case the connections are not random in sign but are organized in an excitatory-inhibitory two-population architecture. Consequently, the balance between excitation and inhibition that gives rise to the temporally disordered state is entirely a dynamic effect. Our results should be contrasted with a conventional, fully connected network with the same simple twopopulation architecture with more conventional 1/K scaling of connections. Such networks converge to either static states or globally coherent limit cycles (Abbott & van Vreeswijk 1993; Gerstner & van Hemmen, 1993; Grannan et al., 1992; Hansel et al., 1995; van Vreeswijk, 1996; Wilson & Cowan, 1972).

An important consequence of our assumption of relatively strong synaptic connections concerns the size of the external input to the local network. According to our theory, the balanced state is robust only when the DC external input to the local network is large-of the same order as the local excitatory and inhibitory feedback, and much larger than the net synaptic input to a cell. In the notation of our model, the external input to an excitatory cell is $\mathrm{Em}_{0} \sqrt{\mathrm{~K}}$, whereas the net input to this cell, $\mathrm{u}_{\mathrm{E}}$, is smaller than the external input by a factor of the order of $1 /\left(\mathrm{m}_{0} \sqrt{\mathrm{~K}}\right)$, where $\mathrm{m}_{0}$ is the rate of an input cell and is assumed to be much larger than $1 / \sqrt{\mathrm{K}}$. Figure3 shows that lowering the strength of the external input-that is, reducing $\mathrm{m}_{0}$-will turn off the activity of the network. In fact, equation 3.3 implies that to maintain the balanced activity in the case of an external input that is only of order 1 requires the vanishing of the denominators in equations 4.3 and 4.4, which means that the interaction strengths have to be fine-tuned to a very narrow range.

Because of the importance of the scaling of the synapses in our theory, it is very informative to consider the behavior of our model if we use a conventional scaling of synapses: that each synapse scales as $1 / \mathrm{K}$-the weak synapses scenario. In this scenario, each component of the synaptic inputs, including the total external input to a cell, is of order 1. The solution of this model (van Vreeswijk & Sompolinsky, 1997) shows that when K is not large, the network settles in a strongly disordered state. This is not surprising given that the connectivity is randomly asymmetric and there is no danger of averaging of the fluctuating inputs to a cell. However, the fate of the temporal variability as K is increased is highly sensitive to the presence of local inhomogeneity. If the neurons have the same threshold, the chaotic state is maintained as K increases. In this case, the population rates adjust their value so that the net input is close to the threshold level within a distance of the order of $1 / \sqrt{\mathrm{K}}$. This is shown in Figure 17. This figure displays the time course of the various synaptic inputs to a cell, evaluated by simulating a sample from the statistics predicted by the mean-field solution with the weak-synapses scaling. The results of this figure should be contrasted with the behavior of the strong synapses scenario (see Figure 6). In Figure 17 the variability is caused by the fact that the cell is always hovering close to its threshold. In contrast, in the case of Figure 6, the distance between the net input and the threshold is not small compared to the distance between threshold and rest. In this case, the variability is caused by the presence of excitatory and inhibitory inputs, each of which is much larger than the threshold.

Despite the difference in behavior between the two scenarios, the dynamic mechanism for these balanced chaotic states is the same. In both cases, the distance between the net input and the threshold is smaller (by

Figure 17: Inputs to an excitatory cell in a network with synaptic strengths $\mathrm{J}_{\mathrm{kl}} / \mathrm{K}$. The total excitatory input (upper trace), the total inhibitory input (lower trace), and the net input (middle trace) are shown in the upper panel. At the bottom, the times when the cell switches from the passive to the active states are indicated. Parameters: $\mathrm{K}=1000, \mathrm{E}=2.0, \mathrm{I}=1.6, \mathrm{~J}_{\mathrm{E}}=2.0, \mathrm{~J}_{\mathrm{I}}=1,8, \theta_{\mathrm{E}}=1.0, \theta_{\mathrm{I}}=0.8$, $\tau=0.9$, and $\mathrm{m}_{0}=0.503$.

a factor of $1 / \sqrt{\mathrm{K}}$ ) from the distance to threshold of the excitatory and inhibitory components. Thus, it would seem that choosing between these scenarios is largely a matter of biological interpretation. However, there are some qualitative differences between the two scenarios. Since the synaptic inputs are all of the same order as the threshold, it is harder to obtain states with low mean rates in both the excitatory and inhibitory populations. To achieve low rates, the ratio between the external inputs to the two populations (I/E in the notation of equations 2.4 and 2.5) has to be close to the ratio of their thresholds $\theta_{\mathrm{I}} / \theta_{\mathrm{E}}$. More important, the weak-synapses scenario of Figure 17 breaks down in the presence of inhomogeneity in the local thresholds. In this case, the population rates are incapable of accommodating the different thresholds. As a result, in the case of inhomogeneous thresholds, when K increases, the network state becomes increasingly frozen; neurons with high thresholds become inactive, whereas neurons with low ones fire close to saturation. This freezing occurs as soon as the width of the inhomogeneity in threshold is larger than $1 / \sqrt{\mathrm{K}}$. In contrast, in the scenario of strong synapses, the state becomes frozen only when the inhomogeneity is large compared to 1. Note that in the case of Figure 17, the external input is of the same order as the net input to the cell. Equally important is the fact that in the weak scenario case and homogeneous networks, the collective time constants are of the same order as the single cell time constant so that the network will not exhibit the phenomenon of fast tracking predicted in our theory.

Finally, the model can be generalized to a model with synaptic strengths that scale as $\mathrm{K}^{-\alpha}$, with $0<\alpha<1$. Of course, these models can be distinguished from the present model only in the large K limit. In this limit, the net average inputs into the populations scale as $\mathrm{K}^{1-\alpha}$, while the quenched and temporal fluctuations in the inputs scale as $\mathrm{K}^{1 / 2-\alpha}$. Therefore the leading order in the inputs has to cancel, leading to the balance condition. For any $\alpha$, this leads to asynchronous chaotic activity in a homogeneous network, similar to the case $\alpha=1 / 2$. However, if we introduce a distribution of the threshold with width of order 1 , we have to distinguish two regimes, apart from $\alpha=1 / 2$ of our model. If $\alpha>1 / 2$, the fluctuations in the input decrease with K so that the network goes to the frozen state in the large K limit. On the other hand, if $\alpha<1 / 2$, the fluctuations grow with K , and therefore the inhomogeneity in the threshold becomes negligible in the large K limit. Thus, a network with inhomogeneous thresholds will act in the same way as a network with homogeneous thresholds. Specifically for low rates, the rate distribution will become narrow. Thus, only for a network with synaptic strengths of order $1 / \sqrt{\mathrm{K}}$ is there a nontrivial interaction between the fluctuation in the input and the threshold inhomogeneities. 10.3 Comparison with Other Netw ork Models. Some of our results are consistent with those of the integrate-and-fire network models of Tsodyks and Sejnowski (1995) and Amit and Brunel (1997a, 1997b). Although constructing an exact mean-field theory for the integrate-and-fire dynamics similar to the one presented here for binary units is much more difficult, we believe that most of the predictions of our mean-field theory are applicable to the integrate-and-fire dynamics as well, provided that the same connectivity architecture and scaling of parameters with N and K are used. However, a direct comparison between our theory and the results of Tsodyks and Sejnowski (1995) and Amit and Brunel (1997a, 1997b) is difficult because of their introduction of stochasticity in the network, the combination of mechanisms such as resetting potential close to threshold, and the lack of full, explicit specification of scaling of parameters with N and K . Tsodyks and Sejnowski (1995) show numerically that their model is capable of "fast switching" in response to a fast change in the external stimulus. This may be related to the fast tracking predicted in our model. The fact that our model does not respond quickly to a sudden switching of the stimulus (see Figure 13) is probably a result of the dynamics of binary neurons. However, the switching time constants observed in Tsodyks and Sejnowski (1995) is of the same order as the single-cell integration time constant, while the fast tracking should occur on a much shorter time constant. In recent numerical simulations of integrate-and-fire networks, Amit and Brunel (1997b) showed that the strength of the average cross-correlations decreases as N increases (keeping the connectivity index constant). However, they do not show whether as N increases, the variability in the single cell remains the same. If this would be the case, their results are consistent with our predictions regarding asynchronous chaotic state. 10.4 Biological Implications. With regard to the biological systems, we should reemphasize that most likely temporal irregularity is a result of several mechanisms, including those mentioned in section 1. Our discussion makes it clear that even with regard to deterministic network mechanism in a finite system, the temporally irregular state is likely to be at best intermediate between the synchronized and the balanced chaotic states. An important question is whether external input is large relative to net input to a cortical cell. Recent experimental findings of Ferster, Chung, and Wheat (1996) in cat primary visual cortex suggest that the input from the lateral geniculate nucleus (LGN) to layer 4 cortical cells are in fact a fraction of the net input. Stratford, Tarczy-Hornoch, Martin, Bannister, and Jack (1996) show that the total strength of the LGN synapses is about 2.5 to 3 times smaller than the total strength of the excitatory feedback synapses from layer 4 cells; however, this study does not measure the strength of the feedback from the inhibitory interneurons, so it does not allow for the estimation of the net feedback. Further experimental clarification of this issue is called for. Measurements of the distribution of time-averaged rates within a local population of neurons and the change in its shape when the overall level of response increases, similar to those of Figure 11B, would be an interesting test of the underlying statistical characteristics of the network spatiotemporal fluctuations. 10.5 Future Work. On theoretical grounds, our work raises several interesting issues worth pursuing. First, it would be important to know whether the theory of the balanced state applies also to networks with more interesting connectivity architecture. Thus, it would be interesting to extend our theory to networks that model associative memory or hypercolumns in visual cortex. It is important to study the consequences of nonlinearities of synaptic summations, for example, by treating synaptic inputs as conductance changes instead of currents.

In considering the functional implications of our theory, it is important to distinguish between the sensitivity of a chaotic autonomous system to changes in its initial condition and its ability to lock to a changing external drive. The analysis of tracking capabilities of our network shows that the macroscopic state of the network responds fast to changing input. In the case of a homogeneous input, it can be shown that the microscopic state is not tightly locked to the changing stimulus. On the other hand, preliminary analysis (van Vreeswijk & Sompolinsky, 1998) shows that in the case of spatially inhomogeneous input fluctuations, the microscopic state of the network will tightly lock to the stimulus temporal variations. These findings are consistent with recent findings that cortical cells respond highly reliably to the fluctuations in the stimulus (Bair & Koch, 1996; Britten, Shadlen, Newsom, & Movshon, 1992). Elucidation of the computational aspects of balanced states in neuronal networks is a challenging issue.

Recently Markram and Tsodyks (1996) have shown that the synapses between cortical pyramidal cells show a marked degree of depression. It should be investigated how such dynamical synapses affect the balanced state. If one assumes synaptic depression between the excitatory-to-excitatory synapses only, and facilitation between the synapses from the inhibitory to the excitatory and from the excitatory to inhibitory populations (Thomson, West, & Deuchars, 1995; Thomson, West, Hahn, & Deuchars, 1996), the equilibrium rates in the network decrease, relative to those in a network without facilitation. This synaptic depression and facilitation also has the effect that the constraints on the synaptic strengths (see equations 4.9 and 4.10) can be relaxed. Because the synaptic depression and facilitation become effective only on a time scale that is as slow as or slower than the membrane time constant, the response of such a network to an external input that changes with time is more complicated than in the model studied here. If the input is suddenly increased by a small amount, the network rates increase to the rate the network would have in equilibrium if the synaptic strengths were not changed in a time of order $1 / \sqrt{\mathrm{K}}$ and then, on a much slower time scale, the rates decrease due to the change in the synaptic strengths.

Since in the balanced state the finite K corrections of the rates are determined by both the first and the second moment of input, the change in rate due to synaptic depression or facilitation depends not only on the average change in the synaptic strength but also on its fluctuation. Thus, synaptic depression due to a change in the height of the excitatory postsynaptive potentials (EPSPs), but without a change in the probability of release, will affect the rates differently from synaptic depression that leaves the height of the EPSPs unaffected but decreases the probability of release, even if both mechanisms result in the same average depression. Another effect of synaptic depression to take into account is that the effect of a spike is decreased if it follows shortly after the preceding spike. This will decrease the fluctuations in the input, relative to the fluctuations in the activity of the presynaptic cells. Facilitation will have the opposite effect, since it increases the effect of a spike if it closely follows the previous one. These issues warrant further study.

A ppendix A: Derivation of the Mean-Field Theory

A. 1 Population Rates. We first consider the population-averaged activities $\mathrm{m}_{\mathrm{E}}(\mathrm{t})$ and $\mathrm{m}_{\mathrm{I}}(\mathrm{t})$ in the limit of large $\mathrm{N}_{\mathrm{E}}$ and $\mathrm{N}_{\mathrm{I}}$ and finite K . We first assume that each cell in the kth population is updated stochastically at a rate $\tau_{\mathrm{k}}$. When a cell is updated, it moves to the active state if its total input is above threshold. Otherwise its updated state is 0 . It is convenient to define a time-dependent local rate variable,

$$\mathrm{m}_{\mathrm{k}}^{\mathrm{i}}(\mathrm{t})=\left\langle\sigma_{\mathrm{k}}^{\mathrm{i}}(\mathrm{t})\right\rangle .$$

Here, the symbol $\langle\ldots\rangle$ does not mean average over time, as in equation 5.4 and thereafter. Instead, it means average over all initial conditions that are consistent with given values for $\mathrm{m}_{\mathrm{k}}(0)$ and also over the random sequence of update times. It is well known that the rate of a binary variable that obeys the update rule satisfies the following continuous time dynamics (Ginzburg & Sompolinsky, 1994; Glauber, 1963),

$$\tau_{\mathrm{k}} \frac{\mathrm{~d}}{\mathrm{dt}} \mathrm{~m}_{\mathrm{k}}^{\mathrm{i}}(\mathrm{t})=-\mathrm{m}_{\mathrm{k}}^{\mathrm{i}}(\mathrm{t})+\Theta\left(\mathrm{u}_{\mathrm{k}}^{\mathrm{i}}(\mathrm{t})\right),$$

where $\mathrm{u}_{\mathrm{k}}^{\mathrm{i}}(\mathrm{t})$ is the total synaptic input into a cell i in the kth population relative to its threshold and is given in our case by equation 2.2. If a cell receives $\mathrm{n}_{\mathrm{E}}(\mathrm{t})$ and $\mathrm{n}_{\mathrm{I}}(\mathrm{t})$ excitatory and inhibitory feedback inputs, respectively, then its input is

$$\mathrm{u}_{\mathrm{k}}^{\mathrm{i}}(\mathrm{t})=\sqrt{\mathrm{K}} \mathrm{~J}_{\mathrm{k} 0} \mathrm{~m}_{0}+\frac{\mathrm{J}_{\mathrm{kE}}}{\sqrt{\mathrm{~K}}} \mathrm{n}_{\mathrm{E}}(\mathrm{t})+\frac{\mathrm{J}_{\mathrm{kI}}}{\sqrt{\mathrm{~K}}} \mathrm{n}_{\mathrm{I}}(\mathrm{t})-\theta_{\mathrm{k}} .$$

The main assumption underlying the mean-field theory is that the activities of the different input cells to a given cell are uncorrelated. Technically, this holds rigorously provided that $\mathrm{K} \ll \log \mathrm{N}_{\mathrm{k}}$ (Derrida et al., 1987). Using this assumption, the population average of equation A. 2 yields the following mean-field equations for the population activities,

$$\tau_{\mathrm{k}} \frac{\mathrm{~d}}{\mathrm{dt}} \mathrm{~m}_{\mathrm{k}}(\mathrm{t})=-\mathrm{m}_{\mathrm{k}}(\mathrm{t})+\mathrm{F}_{\mathrm{k}}\left(\mathrm{~m}_{\mathrm{E}}(\mathrm{t}), \mathrm{m}_{\mathrm{I}}(\mathrm{t})\right),$$

where $\mathrm{F}_{\mathrm{k}}$ denotes the probability that the updating cell at time t will be in an updated active state. It is given by

$$\mathrm{F}_{\mathrm{k}}\left(\mathrm{~m}_{\mathrm{E}}, \mathrm{~m}_{\mathrm{I}}\right)=\sum_{\mathrm{n}_{1}, \mathrm{n}_{2}=0}^{\infty} \mathrm{p}_{1}\left(\mathrm{n}_{1}\right) \mathrm{p}_{2}\left(\mathrm{n}_{2}\right) \Theta\left(\sqrt{\mathrm{K}} \mathrm{~J}_{\mathrm{k} 0} \mathrm{~m}_{0}+\sum_{\mathrm{l}} \frac{\mathrm{~J}_{\mathrm{kl}}}{\sqrt{\mathrm{~K}}} \mathrm{n}_{\mathrm{l}}-\theta_{\mathrm{k}}\right),$$

where $p_{l}(n)$ is the probability that a cell receives $n$ active inputs from the lth population.

For $\mathrm{N}_{\mathrm{E}}, \mathrm{N}_{\mathrm{I}} \rightarrow \infty$ the probability of s synapses of population l projecting to a cell is $\mathrm{K}^{\mathrm{s}} \mathrm{e}^{-\mathrm{K}} / \mathrm{s}$ !. On average each of these synapses has a probability $\mathrm{m}_{1}$ to be active, hence,

$$\begin{aligned} p_{l}(n) & =\sum_{s=n}^{\infty} \frac{K^{s}}{s!} e^{-K}\binom{s}{n} m_{l}^{n}\left(1-m_{l}\right)^{s-n} \\ & =\frac{\left(m_{l} K\right)^{n}}{n!} e^{-m_{l} K} \end{aligned}$$

According to equation A.6, the means and variances of this distribution are $\left[\mathrm{n}_{\mathrm{k}}\right]=\left[\left(\delta \mathrm{n}_{\mathrm{k}}\right)^{2}\right]=\mathrm{Km}_{\mathrm{k}}$. Therefore, in the limit $\mathrm{K} \rightarrow \infty, \mathrm{F}_{\mathrm{k}}\left(\mathrm{m}_{\mathrm{E}}, \mathrm{m}_{\mathrm{I}}\right)$ is given by

$$\mathrm{F}_{\mathrm{k}}\left(\mathrm{~m}_{\mathrm{E}}, \mathrm{~m}_{\mathrm{I}}\right)=\int \mathrm{D} \mathrm{x} \Theta\left(\mathrm{u}_{\mathrm{k}}+\sqrt{\alpha_{\mathrm{k}}} \mathrm{x}\right)=\mathrm{H}\left(\frac{-\mathrm{u}_{\mathrm{k}}}{\sqrt{\alpha_{\mathrm{k}}}}\right)$$

where $\mathrm{Dx}=\mathrm{dx} \exp \left(-\mathrm{x}^{2} / 2\right) / \sqrt{2 \pi}$. From the above statistics of $\mathrm{n}_{1}$, one obtains that the average input, relative to threshold, $\mathrm{u}_{\mathrm{k}}$ into a cell of population k given by

$$\mathrm{u}_{\mathrm{k}}=\left(\mathrm{J}_{\mathrm{k} 0} \mathrm{~m}_{0}+\mathrm{J}_{\mathrm{kE}} \mathrm{~m}_{\mathrm{E}}+\mathrm{J}_{\mathrm{kl}} \mathrm{~m}_{\mathrm{I}}\right) \sqrt{\mathrm{K}}-\theta_{\mathrm{k}}$$

and standard deviation of the input $\alpha_{\mathrm{k}}$,

$$\alpha_{\mathrm{k}}=\left(\mathrm{J}_{\mathrm{kE}}\right)^{2} \mathrm{~m}_{\mathrm{E}}+\left(\mathrm{J}_{\mathrm{kI}}\right)^{2} \mathrm{~m}_{\mathrm{I}}$$

from which equations 3.3 through 3.6 follow. A. 2 Autocorrelations. We now extend the analysis to evaluate the dynamics of the autocorrelation function $\mathrm{q}_{\mathrm{k}}(\tau)$, (see equations 5.16). Using similar arguments as for equations A.2, $\mathrm{q}_{\mathrm{k}}(\tau)$ satisfies an equation of the following form,

$$\tau_{\mathrm{k}} \frac{\mathrm{dq}_{\mathrm{k}}}{\mathrm{~d} \tau}=-\mathrm{q}_{\mathrm{k}}(\tau)+\int_{0}^{\infty} \frac{\mathrm{dt}^{\prime}}{\tau_{\mathrm{k}}} \exp \left(-\mathrm{t}^{\prime} / \tau_{\mathrm{k}}\right) \mathrm{F}_{\mathrm{k}}\left(\left\{\mathrm{~m}_{\mathrm{l}}\right\} ;\left\{\mathrm{q}_{\mathrm{l}}\left(\mathrm{t}^{\prime}+\tau\right)\right\}\right),$$

where

$$\mathrm{F}_{\mathrm{k}}=\left[\left\langle\Theta\left(\mathrm{u}_{\mathrm{k}}^{\mathrm{i}}(\mathrm{t})\right) \Theta\left(\mathrm{u}_{\mathrm{k}}^{\mathrm{i}}\left(\mathrm{t}+\mathrm{t}^{\prime}+\tau\right)\right)\right\rangle\right] .$$

Here the averaging is also over absolute timet. The integral over time in the r.h.s. of equation A. 11 takes into account the correlation between the inputs to a cell that updates its state at time $\mathrm{t}+\tau$ and its inputs at the last update before time t . Thus, the time integral is an integral over the exponential distribution of update interval of the last update before time t. Separating the total number of active inputs into those that come from sources in the kth population that are active in both times ( $\mathrm{n}_{1, \mathrm{k}}$ ) and those that are active only in one of the times ( $\mathrm{n}_{2, \mathrm{k}}$ and $\mathrm{n}_{3, \mathrm{k}}$, respectively), one can write

$$\begin{aligned} \mathrm{F}_{\mathrm{k}}= & \prod_{\mathrm{l}=1,2} \sum_{\mathrm{n}_{\mathrm{kl}}} \mathrm{p}_{\mathrm{l}}\left(\mathrm{n}_{1 \mathrm{l}}, \mathrm{n}_{2 \mathrm{l}}, \mathrm{n}_{3 \mathrm{l}}\right) \Theta\left(\sqrt{\mathrm{K}} \mathrm{~J}_{\mathrm{k} 0} \mathrm{~m}_{0}+\sum_{\mathrm{l}} \frac{\mathrm{~J}_{\mathrm{kl}}}{\sqrt{\mathrm{~K}}}\left(\mathrm{n}_{1 \mathrm{l}}+\mathrm{n}_{2 \mathrm{l}}\right)-\theta_{\mathrm{k}}\right) \\ & \times \Theta\left(\sqrt{\mathrm{K}} \mathrm{~J}_{\mathrm{k} 0} \mathrm{~m}_{0}+\sum_{\mathrm{l}} \frac{\mathrm{~J}_{\mathrm{kl}}}{\sqrt{\mathrm{~K}}}\left(\mathrm{n}_{1 \mathrm{l}}+\mathrm{n}_{3 \mathrm{l}}\right)-\theta_{\mathrm{k}}\right) \end{aligned}$$

where

$$\mathrm{p}_{\mathrm{l}}\left(\mathrm{n}_{1}, \mathrm{n}_{2}, \mathrm{n}_{3}\right)=\frac{\left(\mathrm{q}_{1} \mathrm{~K}\right)^{\mathrm{n}_{1}}}{\mathrm{n}_{1}!} \frac{\left(\left(\mathrm{m}_{1}-\mathrm{q}_{1}\right) \mathrm{K}\right)^{\mathrm{n}_{2}+\mathrm{n}_{3}}}{\mathrm{n}_{2}!\mathrm{n}_{3}!} \mathrm{e}^{-\left(2 \mathrm{~m}_{1}-\mathrm{q}_{1}\right) \mathrm{K}}$$

In the large K limit, this can be written as

$$\begin{aligned} \mathrm{F}_{\mathrm{k}}= & \int \mathrm{D} \mathrm{x}_{1} \int \mathrm{D} \mathrm{x}_{2} \int \mathrm{D} \mathrm{x}_{3} \Theta\left(\mathrm{u}_{\mathrm{k}}+\sqrt{\beta_{\mathrm{k}}} \mathrm{x}_{1}+\sqrt{\alpha_{\mathrm{k}}-\beta_{\mathrm{k}}} \mathrm{x}_{2}-\theta_{\mathrm{k}}\right) \\ & \times \Theta\left(\mathrm{u}_{\mathrm{k}}+\sqrt{\beta_{\mathrm{k}}} \mathrm{x}_{1}+\sqrt{\alpha_{\mathrm{k}}-\beta_{\mathrm{k}}} \mathrm{x}_{3}-\theta_{\mathrm{k}}\right) \\ = & \int \mathrm{Dx}\left[\mathrm{H}\left(\frac{\theta_{\mathrm{k}}-\mathrm{u}_{\mathrm{k}}-\sqrt{\beta_{\mathrm{k}}} \mathrm{x}}{\sqrt{\alpha_{\mathrm{k}}-\beta_{\mathrm{k}}}}\right)\right]^{2} \end{aligned}$$

with $\mathrm{u}_{\mathrm{k}}$ and $\alpha_{\mathrm{k}}$ as above, and $\beta_{\mathrm{k}}$ given by

$$\beta_{\mathrm{k}}(\tau)=\sum_{\mathrm{l}=1,2}\left(\mathrm{~J}_{\mathrm{kl}}\right)^{2} \mathrm{q}_{\mathrm{l}}(\tau)$$

So $q_{k}$ satisfies equation 5.17. A. 3 Sensitivity to Initial Conditions. The derivation of equation 8.5 for the overlaps $\mathrm{Q}_{\mathrm{k}}(\mathrm{t})=\left(\mathrm{m}_{\mathrm{k}}+\mathrm{D}_{\mathrm{k}}\right) / 2$ of two trajectories, equation 8.2 which start with slightly different initial conditions, is similar to that of $\mathrm{q}_{\mathrm{k}}$. Here the inputs $\mathrm{n}_{1, \mathrm{k}}$ are the sources to a given cell that are active at timet in both trajectories. The only difference between the equation for the delayed-time autocorrelations and the equal-time overlap between two trajectories is the integral over the previous update times which appears in equations A. 10 and 5.17. This results from the fact that in the latter case, the update sequence is identical in the two trajectories.

Appendix B: Determinstic Update Rules

The general form of equation A. 4 is usually derived for a binary variable that is updated stochastically at a rate $\tau_{\mathrm{k}}$. One might therefore argue that the irregular firing in our model is due to the stochasticity of the update times of the model neurons. To show that this is not the case, we define here a completely deterministic dynamic model and show that it leads to exactly the same equations for the mean rates of activity as those given above.

Consider the same network model, except that a neuron i of population k is updated at times $\mathrm{t}=\left(\mathrm{n}+\delta_{\mathrm{k}}^{\mathrm{i}}\right) \tau_{\mathrm{k}}$ with $\mathrm{n}=0,1,2, \ldots$ and $\delta_{\mathrm{k}}^{\mathrm{i}}$ is randomly chosen between 0 and 1 . Let $\mathrm{m}_{\mathrm{k}}^{+}(\mathrm{t})$ be the probability that the neuron of population k, updated at timet, goes into (or stays in) the active state. Since all neurons of population k are updated exactly once between times $\mathrm{t}-\tau_{\mathrm{k}}$ and $\mathrm{t}, \mathrm{m}_{\mathrm{k}}(\mathrm{t})$ is given by

$$\mathrm{m}_{\mathrm{k}}(\mathrm{t})=\frac{1}{\tau_{\mathrm{k}}} \int_{0}^{\tau_{\mathrm{k}}} \mathrm{dt}^{\prime} \mathrm{m}_{\mathrm{k}}^{+}\left(\mathrm{t}-\mathrm{t}^{\prime}\right)$$

Going through arguments similar to those shown above, one can show that $\mathrm{m}_{\mathrm{k}}$ satisfies

$$\mathrm{m}_{\mathrm{k}}(\mathrm{t})=\frac{1}{\tau_{\mathrm{k}}} \int_{0}^{\tau_{\mathrm{k}}} \mathrm{dt}^{\prime} \mathrm{F}_{\mathrm{k}}\left(\mathrm{~m}_{\mathrm{E}}\left(\mathrm{t}-\mathrm{t}^{\prime}\right), \mathrm{m}_{\mathrm{I}}\left(\mathrm{t}-\mathrm{t}^{\prime}\right)\right)$$

with $\mathrm{F}_{\mathrm{k}}$ given by equation A.5. If we introduce inhomogeneities in the rate with which the cells are updated so that cell i of population is updated at times $\mathrm{t}=\left(\mathrm{n}+\delta_{\mathrm{k}}^{\mathrm{i}}\right) \tau_{\mathrm{i}}^{\mathrm{k}}$, where $\tau_{\mathrm{i}}^{\mathrm{k}}$ has a probability $\mathrm{R}_{\mathrm{k}}(\tau) \mathrm{d} \tau$ of being between $\tau$ and $\tau+\mathrm{d} \tau$, we find that $\mathrm{m}_{\mathrm{k}}$ evolves as

$$\mathrm{m}_{\mathrm{k}}(\mathrm{t})=\int_{0}^{\infty} \mathrm{d} \tau \frac{\mathrm{R}_{\mathrm{k}}(\tau)}{\tau} \int_{0}^{\tau} \mathrm{dt}^{\prime} \mathrm{F}_{\mathrm{k}}\left(\mathrm{~m}_{\mathrm{E}}\left(\mathrm{t}-\mathrm{t}^{\prime}\right), \mathrm{m}_{\mathrm{I}}\left(\mathrm{t}-\mathrm{t}^{\prime}\right)\right)$$

For $\mathrm{R}_{\mathrm{k}}(\mathrm{t})=\mathrm{te}^{-\mathrm{t} / \tau_{\mathrm{k}}} / \tau_{\mathrm{k}}^{2}$, this can be written as

$$\mathrm{m}_{\mathrm{k}}(\mathrm{t})=\frac{1}{\tau_{\mathrm{k}}} \int_{0}^{\infty} \mathrm{dt}^{\prime} \mathrm{e}^{-\mathrm{t} / \tau_{\mathrm{k}}} \mathrm{~F}_{\mathrm{k}}\left(\mathrm{~m}_{\mathrm{E}}\left(\mathrm{t}-\mathrm{t}^{\prime}\right), \mathrm{m}_{\mathrm{I}}\left(\mathrm{t}-\mathrm{t}^{\prime}\right)\right)$$

and this is equivalent to equation A.4. Thus, in this completely deterministic model, the mean rates $\mathrm{m}_{\mathrm{k}}$ satisfy exactly the same equations as the model with stochastic updating. This also holds true for the other mean-field equations of the model.

Acknowledgments

We thank D. J. Amit, D. Hansel, T. Sejnowski, and M. Tsodyks for very helpful discussions. We are grateful to M. Abeles, H. Bergman, and E. Vaadia for permission to present their data. This work is partially supported by the Fund for Basic Research of the Israeli Academy of Science.

References

Abbott, L. F., & van Vreeswijk, C. (1993). Asynchronous states in networks of pulse-coupled oscillators. Phys. Rev. E, 48, 1483-1490. Abeles, M. (1991). Corticonics: Neural circuits of the cerebral cortex. Cambridge: Cambridge University Press. Abeles, M., Bergman, H., & Vaadia, E. (1988). Unpublished data. Amit, D. J., &Brunel, N. (1997a). Model of global spontaneous activity and local structured activity during delay periods in the cerebral cortex Cereb. Cortex, 7, 237-252. Amit, D. J., & Brunel, N. (1997b). Dynamics of a recurrent network of spiking neurons before and following learning. Network, 8, 373-404. Bair, W., & Koch, C. (1996). Temporal precision of spike trains in extrastriate cortex of the behaving Macaque monkey. Neural Comput., 8, 1185-1202. Bair, W., Koch, C., Newsome, W., &Britten, K. (1994). Power spectrum analysis of bursting cells in area MT in the behaving monkey. J. Neurosci., 14, 2870-2892. Bell, A., Mainen, Z. F., Tsodyks, M., & Sejnowski, T. (1994). Why do cortical neurons fire irregularly? Soc. Neurosci. Abstr., 20, 1527. Britten, K. H., Shadlen, M. J., Newsome, W. T., & Movshon, J. A. (1992). The analysis of visual motion: A comparison of neuronal and psychophysical performance. J. Neurosci., 12, 4745-4765. Burns, B. D., & Webb, A. C. (1976). The spontaneous activity of neurons in the cat's visual cortex. Proc. R. Soc. Lond. B, 194, 211-223. Bush, P. C., & Douglas, R. J. (1991). Synchronization of bursting action potential discharge in a model network of neocortical neurons. Neural Comp., 3, 19-30. Derrida, B., Gardner, E., & Zippelius, A. (1987). An exactly soluble asymmetric neural network model. Europhys. Lett., 4, 167-173. Douglas, R.J., &Martin, K.A. C. (1991). Opening the grey box. Trends in Neurosci., 14, 286-293. Douglas, R.J., Martin, K.A. C., &Whitteridge, D. (1991). An intracellularanalysis of the visual response of neurons in cat visual cortex. J. Physiol., 440, 659-696. Ermentrout, J. B., & Gutkin, B. (in press). Dynamics of spike generation determines in vivo spike train statistics. Ferster, D., Chung, S., & Wheat, H. (1996). Orientation selectivity of thalamic input to simple cells of cat visual cortex. Nature, 380, 249-252. Ferster, D., & Jagadeesh, B. (1992). EPSP-IPSP interactions in cat visual cortex studied with in vivo whole-cell patch recording. J. Neurosci., 12, 1262-1274. Gerstein, G. L., & Mandelbrot, B. (1964). Random walk models for the spike activity of a single cell. Biophys. J., 4, 41-68. Gerstner, W., & van Hemmen, J. L. (1993). Associative memory in a network of "spiking" neurons. Network, 3, 139-164. Ginzburg, I., & Sompolinsky, H. (1994). Theory of correlations in stochastic neural networks. Phys. Rev. E, 50, 3171-3190.

Glauber, R.J. (1963). Time-dependent statistics of the Ising model. J. Math. Phys., 4, 294-307. Grannan, E., Kleinfeld, D., & Sompolinsky, H. (1992). Stimulus-dependent synchronization of neuronal assemblies. Neural Comp., 4, 550-569. Gray, C. M., &Singer, W. (1989). Oscillatory responses in cat visual cortex exhibit inter-columnar synchronization which reflects global stimulus properties. Nature, 338, 334-337. Hansel, D., Mato, G., & Meunier, C. (1995). Synchronization in excitatory neural networks. Neural Comp., 3, 307-338. Hansel, D., & Sompolinsky, H. (1992). Synchronization and computation in a chaotic neural network. Phys. Rev. Lett., 68, 718-721. Hansel, D., & Sompolinsky, H. (1996). Chaos and synchrony in a model of a hypercolumn in visual cortex. J. Comput. Neurosci., 3, 7-34. Holt, G. R., Softky, W. R., Koch, C., & Douglas, R. J. (1996). A comparison of discharge variability in vitro and in vivo in cat visual cortex neurons. J. Neurophysiol., 75, 1806-1814. Mainen, Z. J., & Sejnowski, T. (1995). Reliability of spike timing in neocortical neurons. Science, 268, 1503-1506. Markram, H., & Tsodyks, M. (1996). Redistribution of synaptic efficacy between neocortical pyramidal neurons. Nature, 382, 807-810. Mézard, M., Parisi, G., & Virasoro, M. A. (1987). Spin glass theory and beyond. Singapore: World Scientific. Perkel, D. H., Gerstein, G. L., & Moore, G. P. (1967a). Neuronal spike trains and stochastic pointprocesses. I. The single spike train. Biophys. J., 7, 391-418. Perkel, D. H., Gerstein, G. L., & Moore, G. P. (1967b). Neuronal spike trains and stochastic point processes. II. Simultaneous spike trains. Biophys. J., 7, 419-440. Shadlen, M. N., & Newsome, W. T. (1994). Noise, neural codes and cortical organization. Curr. Opin. Neurobiol., 4, 569-579. Shadlen, M. N., & Newsome, W. T. (1995). Is there a signal in the noise? Curr. Opin. Neurobiol., 5, 248-250. Softky, W. R. (1995). Simple codes versus efficient codes. Curr. Opin. Neurobiol., 5, 239-247. Softky, W. R., Koch, C. (1993). The highly irregular firing of cortical cells is inconsistent with temporal integration of random EPSPs. J. Neurosc., 13, 334350.

Sompolinsky, H., Crisanti, A., &Sommers, H.J. (1988). Chaos in neural networks Phys. Rev. Lett., 61, 259-262. Stratford, K. J., Tarczy-Hornoch, K., Martin, K. A. C., Bannister, N. J., & Jack, J. J. B. (1996). Excitatory synaptic inputs to spiny stellate cells in cat visual cortex. Nature, 382, 258-261. Thomson, A. M., West, D. C., & Deuchars, J. (1995). Properties of single axon excitatory postsynaptic potentials elicited in spiny interneurons by action potentials in pyramidal neurons in slices of rat neocortex. Neurosci., 69, 727738.

Thomson, A. M., West, D. C., Hahn, J., & Deuchars, J. (1996). Single axon IPSPs elicited in pyramidal cells by three classes of interneurones in slice of rat neocortex. J. Physiol., 496, 81-102. Troyer, T. W., &Miller, K. D. (1997). Physiological gainleads to high ISI variability in a simple model of a cortical regular firing cell. Neural Comput., 9, 733-745. Tsodyks, M., Mitkov, I., & Sompolinsky, H. (1993). Pattern of synchrony in inhomogeneous networks of oscillators with pulse interaction. Phys. Rev. Lett., 71, 1282-1285. Tsodyks, M., & Sejnowski, T. (1995). Rapid state switching in balanced cortical network models. Network, 6, 111-124. Vaadia, E., Haalman, I., Abeles, M., Prut, Y., Slovin, H., & Aertsen, A. (1995). Dynamics of neural interaction in monkey cortex in relation to behavioral events. Nature, 373, 515-518 van Vreeswijk, C. (1996). Partial synchrony in populations of pulse-coupled oscillators. Phys. Rev. E, 54, 5522-5537. van Vreeswijk, C., & Sompolinsky, H. (1996). Chaos in neuronal networks with balanced excitatory and inhibitory activity. Science, 274, 1724-1726. van Vreeswijk, C., & Sompolinsky, H. (1997). Irregular firing in sparse networks of weakly coupled neurons. Unpublished manuscript. van Vreeswijk, C., & Sompolinsky, H. (1998). Locking to fluctuating input in balanced networks. Unpublished manuscript. Wilson, H. R., & Cowan, J. D. (1972). Excitatory and inhibitory interactions of localized populations of model neurons. Biophys. J., 12, 1-24.

Received April 9, 1997; accepted December 10, 1997.