We model neurons as integrate-and-fire (I&F) units (Dayan and Abbott, 2001; Newhall et al., 2010; Zhou et al., 2010). The governing equation for the membrane potential vik of the ith neuron in the kth population is

where g_L denotes the leakage conductance, ϵRk is the resting voltage of the kth population and Iik(t) is the corresponding input current (k = E, I). The voltage vik evolves according to Equation (1) when vik≤ϵTk, where ϵTk is the threshold of the kth population. When vik crosses ϵTk, the neuron spikes, and then vik is reset to ϵRk. Upon resetting, vik is immediately governed by Equation (1) again. In simulation, gL=50 s-1 corresponds to the membrane time constant 20 ms. The dimensionless values used in simulations are ϵRE=ϵRI=0.0, ϵTE=1.0 and ϵTI=0.7, which correspond to the parameters in Vreeswijk and Sompolinsky (1998).

The instantaneous current projecting into the ith neuron in the kth population,

consists of two terms, where IikE(t)=fk∑sδ(t-ζisk)+JkE∑j=1NECijkE∑sδ(t-τjsE) is the excitatory input, IikI(t)=-JkI∑j=1NICijkI∑sδ(t-τjsI) is the inhibitory input, δ(·) is the Dirac delta function, f^k is the strength of the external input, and J^kl is the coupling strength from the lth population to the kth population, k, l = E, I. The coupling constant Cijkl=0 or 1 is an element of the adjacency matrix of the network. It describes the connection from the jth neuron in the lth population to the ith neuron in the kth population. The first term in IikE(t) corresponds to the current arriving from the external input. The spike-time sequence, {ζisk,s=1,2,…}, corresponds to the external input into the ith neuron in the kth population. At the arrival time ζisk of the sth spike, the voltage of the ith neuron in the kth population jumps by the amount of f^k. In the simulations below, we use Poisson trains for the external inputs. The second term in IikE(t) and the term in IikI(t) correspond to the presynaptic input from neurons of the excitatory and inhibitory populations in the network, respectively, where τjsk is the arrival time of the sth spike from the jth neuron in the kth population for k = E, I.

In our simulation, each neuron in the network has, on average, K excitatory and K inhibitory presynaptic neurons. Because each neuron may receive a large number of synapses in the cortex (Peters, 1987; Braitenberg and Schuz, 1998), and the connection between cortical neurons is often sparse with low connection probability (Holmgren et al., 2003), we choose K to be sufficiently large but much smaller than the total number of neurons in the network. Experimentally, for example, it is observed that cells in the primary visual cortex of adult cats fire much more irregularly than cells in vitro when they are both stimulated by injecting direct current through the electrode. Therefore, there is a substantial influence of fluctuations of synaptic inputs on the irregular activity (Holt et al., 1996). To capture the effect of fluctuations as observed in the experiment, we follow the balanced network theory to set the scaling of the coupling strength J^kl to be of order 1/K, thus leading to the scaling of fluctuations in the total synaptic inputs as of order 1. As a consequence, the fluctuations persist in the large-K limit (van Vreeswijk and Sompolinsky, 1996; Vreeswijk and Sompolinsky, 1998; Vogels and Abbott, 2005).

The connection from the jth neuron in the lth population to the ith neuron in the kth population Cijkl in our simulation follows a Bernoulli distribution, i.e., the probability P(Cijkl=1)=K/Nl and P(Cijkl=0)=1-K/Nl, where N^l is the total number of the lth population, k, l = E, I. The parameter values in our simulations are as follows: N^E = 32000, N^I = 8000, K = 400, JEE=JIE=1.0/K, JII=1.8/K, JEI=2.0/K, fE=1.0/K, fI=0.8/K, and the external inputs are Poisson processes with rate ν^E = ν^I = ν⁰K, where ν⁰ controls the magnitude of Poisson rate.

Simulations of the neuronal network model are carried out to the machine accuracy using the event-driven algorithm (Brette et al., 2007). The event-driven algorithm for a general point process as external inputs proceeds by generating the time of the next external spike. Some neurons in the network receive this spike. Upon receiving an external spike, these neurons' voltages increase by the amount of f^E for the excitatory population and f^I for the inhibitory population. If some of them have their voltages exceed the threshold, they fire. Their voltages are held at the reset voltage, and the voltages of their postsynaptic neurons are instantaneously increased (for excitatory inputs) or decreased (for inhibitory inputs). It is possible that the voltages of these postsynaptic neurons may now be above threshold as well, then, these neurons also fire. Their voltages are held at the reset voltage, while their postsynaptic neurons' voltages are changed. This process repeats until no new neurons spike. We emphasize that in our dynamics we hold the voltage of the neurons that just fired at the resting potential in order to prevent any of these neurons from firing more than once at any given time. After all the neurons are updated at this time, we release these neurons from the reset voltage to follow the dynamics governed by Equation (1) until the next external spike.

As is well known, for the balanced state in binary neuronal networks, the population-averaged firing rate is a linear function of the external input. We now address the question of whether the balanced state in the I&F network also possesses the linear response property of the population-averaged firing rate to the external drive.

We consider the mean firing rate m_E and m_I in the large-K limit. In a balanced network, the firing events of different neurons are nearly independent of one another (Vreeswijk and Sompolinsky, 1998). This is illustrated in Figure 1 in which the cross-correlation between firing events of pairs of neurons is narrowly distributed around zero. According to Equations (1) and (2), it is obvious that, under Poisson-train external inputs, spiking events of a neuron in the network, in general, are not a Poisson train, i.e., {τjsE} and {τjsI} do not follow Poisson point process for a fixed j in Equation (2). However, the input to the ith neuron is a spike train summed over outputs from many neurons in the network. Since the firing events of neurons are statistically independent from one another, the summed spike train from a large number of output spike trains of neurons in the network asymptotically approaches a Poisson spike process (Cinlar, 1972). The recurrent excitatory (inhibitory) input of each neuron can be treated as a Poisson train with rate Km^E (Km^I), where m^E (m^I) is the mean firing rate per neuron averaged over the excitatory (inhibitory) population. Each neuron receives three Poisson spike trains in Equation (2). By homogeneity of our networks, the subscript i in Equation (2) will be dropped for the remainder of this discussion.

Since the voltage of neuron in the kth population is reset to ϵRk after spiking, we consider Equation (1) with initial value vk(0)=ϵRk for k = E, I. We can readily obtain vk(t)=ϵRk+fkvk,ext(t)+JkEvk,E(t)-JkIvk,I(t) with vk,l(t)=∑s=1Mk,le-gL(t-Usk,l), where M^k,l is the total spike count of the lth input to the kth population, described by a Poisson distribution with the average number of successes ν^k,lt, k = E, I and l = ext, E, I. For each given M^k,l, the spike times of the lth input to the kth population, Usk,l, s = 1, 2, …, are uniformly distributed on the interval [0, t]. Clearly, the random variable Rsk,l(t)≡e-gL(t-Usk,l) takes value in the interval [e-gLt,1] with the following probability density PR(t)(r)=1gLt1r for r∈[e-gLt,1]. Since v^k,l(t) is given by a sum of independent identically distributed random variables for given M^k,l, the average neuronal voltage at time t can be simply expressed as

for k = E, I. As discussed above, ϵRk=0, f^k, J^kE, and J^kI are of order 1/K; ν^k is of order K; and m^E and m^I are of order 1. Then, the leading order of u^k(t) is K. Obviously, the membrane potential cannot become infinite as K → +∞. Therefore, the leading order of K should vanish and one can obtain

in the large-K limit. Equation (4) describes the linear dependence of the mean firing rate on the external input. This relation is a defining feature of a balanced state, similar to what is obtained for the binary neuronal system (Vreeswijk and Sompolinsky, 1998).

In Vreeswijk and Sompolinsky (1998), the properties of the balanced state are shown in detail with binary neuronal networks. We use numerical simulations to investigate whether the current-based I&F neuronal network coupled with delta-pulse interactions can exhibit the dynamical characteristics of a balanced state. Our results demonstrate that there indeed exists a balanced state in I&F neuronal networks. Figure 1 summarizes the defining characteristics of the balanced state as exhibited in I&F neuronal networks: balanced inputs (Figure 1A), irregular spiking activity (Figure 1B), heterogeneous firing rate (Figure 1C), weak correlation (Figure 1D), stationary asynchronous dynamics (Figure 1E), and linear response (Figure 1F). The results above confirm the important properties of balanced networks as discussed in previous theoretical work (van Vreeswijk and Sompolinsky, 1996; Vreeswijk and Sompolinsky, 1998; Renart et al., 2010; Litwin-Kumar and Doiron, 2012).

Next, we turn to the investigation of the persistence of the balanced state in the I&F system and derive conditions under which the balanced network state exists. Note that, in our simulation, we choose J^IE = J^EE. Therefore, from Equation (4), requiring firing rates to be non-negative implies

or

Note that Equation (5) is equivalent to the balance condition derived in Vreeswijk and Sompolinsky (1998) for a binary neuronal system. However, Equation (6) admits a solution for which m^E = 0 — that is, only inhibitory neurons fire in the system, the mean firing rate of the inhibitory population is then given by m^I = f^Iν⁰/J^II. Therefore, for k = E in Equation (3), we can find uE(t)≈ϵRE-CK, where C is a positive constant independent of K, which demonstrates that membrane potential is highly negative, and the excitatory neurons' firing activity is suppressed. In the simulation, we choose a range of parameters, in particular, f^E, f^I, J^EI, and J^II to examine the competition between the excitatory and inhibitory inputs to a neuron and verify whether the system maintains the characteristic balanced-state properties. If the neuronal network dynamics indeed exhibits the features displayed in Figure 1, that network is classified as balanced. Figure 2A summarizes our results for the parameters f^E/f^I and J^EI/J^II that we have scanned. Each dot in the figure can represent a set of parameters f^E, f^I, J^EI, and J^II with fixed ratios of f^E/f^I and J^EI/J^II. Each red dot in the parameter space indicates the fact that the system with the corresponding parameters can reach a balanced state; each blue dot represents the system with the corresponding parameters that exhibits synchronous dynamics; and each green dot indicates that at those parameter values only inhibitory neurons fire in the system but the inhibitory population still retains the characteristic balanced-state properties in Figure 1.

The balanced state with only the inhibitory population firing is established by the balance between strong excitatory external input to the inhibitory neurons and strong inhibitory recurrent current. As shown in Figure 2C, in this state, only the inhibitory population exhibits asynchronous firing activity while the excitatory population is silenced due to relatively strong recurrent inhibition as compared with the inhibitory population.

The parameter values for the results of the representative case as presented in Figure 2B are marked by the black asterisk in the red-dot area in Figure 2A. The parameters for the representative case shown in Figure 2C are marked by the black asterisk in the green-dot area in Figure 2A. The broad region covered by the red and green dots demonstrates the existence of the balanced state over a wide range of parameter space. The parameter values marked by the black asterisk in the blue-dot area in Figure 2A with strong f^E and J^II in comparison to f^I and J^EI, respectively, render the system robustly synchronized as shown in Figure 2D.

Note that the area represented by Equation (5) is smaller than the red-dot region that supports stable balanced states in our numerical results. We also point out that the area described by Equation (6) contains a small number of stable balanced states in which both excitatory and inhibitory neurons are spiking in addition to balanced states in which excitatory neurons are almost silent while only inhibitory neurons are spiking. As a matter of fact, upon inserting m^E = 0 and m^I = f^Iν⁰/J^II into Equation (3), and requiring that u^E becomes sufficiently negative as K becomes large so the excitatory population rarely fires, we obtain the condition

which contains Equation (6) and covers a greater area than the parameter regions in Figure 2A filled with green dots. These differences between the numerical results and the theoretical conditions (Equations 5–7) arise from the finite size effect of K.

As shown above, the I&F networks with delta-pulse coupling can persistently manifest the dynamics of a balanced state. We now address our central question of whether the irregular firing activity of neurons in the balanced network is a consequence of chaotic dynamics of the network. We can mathematically prove that the I&F networks with delta-pulse coupling cannot exhibit chaotic dynamics. Therefore, chaos may not underpin the irregular firing activity of neurons in a balanced state. Below is a proof that there is no chaos in the dynamics of the current-based I&F network coupled with delta-pulses and with a general point process (not limited to Poisson) as external inputs, the details of whose dynamics are described in the section Materials and Methods.

For each reference voltage trajectory, v(t)=(v1E(t),v2E(t),…,vNEE(t),v1I(t), v2I(t),…,vNII(t)), we consider the perturbed voltage trajectory ṽ(t)=(ṽ1E(t),ṽ2E(t),…,ṽNEE(t),ṽ1I(t),ṽ2I(t),… ,ṽNII(t)), with sufficiently small perturbation size ϵ at initial time t = t₀, i.e., ϵ = |ṽ(t₀) − v(t₀)| ≪ 1. The dynamics of the perturbed trajectory ṽ(t) is described by the equation

where τ~jsk is the sth spike of the jth neuron in the kth population along ṽ(t), and k = E, I. The external spike times {ζisk} are the same as those along the reference trajectory v(t). Then the largest Lyapunov exponent,

can be calculated for this I&F network dynamics (Zhou et al., 2009). As is well known, positive λ_max measures the average exponential spreading of nearby trajectories, while negative λ_max measures the exponential convergence of trajectories onto the attractor (Oseledec, 1968; Ott, 2002; Parker and Chua, 2012). Generically, an attractor is defined to be non-chaotic if λ_max is non-positive. In what follows, we show that the largest Lyapunov exponent λ_max of this I&F network is always negative and in fact approaches negative infinity for any spike train input.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.