Skip to main content

ORIGINAL RESEARCH article

Front. Comput. Neurosci., 28 June 2018

The Dynamics of Balanced Spiking Neuronal Networks Under Poisson Drive Is Not Chaotic

  • 1School of Mathematical Sciences, MOE-LSC, Institute of Natural Sciences, Shanghai Jiao Tong University, Shanghai, China
  • 2Mathematical Sciences Department, Rensselaer Polytechnic Institute, Troy, NY, United States
  • 3Courant Institute of Mathematical Sciences, Center for Neural Science, New York University, New York, NY, United States
  • 4NYUAD Institute, New York University Abu Dhabi, Abu Dhabi, United Arab Emirates

Some previous studies have shown that chaotic dynamics in the balanced state, i.e., one with balanced excitatory and inhibitory inputs into cortical neurons, is the underlying mechanism for the irregularity of neural activity. In this work, we focus on networks of current-based integrate-and-fire neurons with delta-pulse coupling. While we show that the balanced state robustly persists in this system within a broad range of parameters, we mathematically prove that the largest Lyapunov exponent of this type of neuronal networks is negative. Therefore, the irregular firing activity can exist in the system without the chaotic dynamics. That is the irregularity of balanced neuronal networks need not arise from chaos.

1. Introduction

Neural spiking activity in the brain is highly irregular (Britten et al., 1993; Shadlen and Newsome, 1998; Compte et al., 2003; London et al., 2010). It is believed that the irregularity of the spiking activity can reflect an underlying rich coding structure for information processing (Hertz and Prügel-Bennett, 1996; Gütig and Sompolinsky, 2006; Sussillo and Abbott, 2009; Monteforte and Wolf, 2012). This viewpoint naturally leads to the investigation of the origin of the irregular neuronal activity. A number of theoretical studies have postulated that a balance between excitatory and inhibitory inputs into an individual neuron can give rise to irregular activity (van Vreeswijk and Sompolinsky, 1996; Troyer and Miller, 1997; Vreeswijk and Sompolinsky, 1998; Vogels and Abbott, 2005; Miura et al., 2007). The idea behind the theory of balanced networks is that the excitatory and inhibitory components of inputs nearly cancel each other, and the neuronal firing activity is driven by strong fluctuations that intermittently interrupt this cancellation. Consistent with the hypothesized scenario, balanced synaptic inputs have been observed in slices of the ferret prefrontal and occipital cortex (Shu et al., 2003). Moreover, it has also been found that, in vivo studies, the balanced excitation and inhibition in ferrets' prefrontal cortex can substantially influence the neuronal activity (Haider et al., 2006).

It has been shown theoretically that small perturbations of the balanced state in a network with binary neurons grow exponentially, indicating the chaotic nature of the balanced activity (Vreeswijk and Sompolinsky, 1998). Some studies suggest that highly irregular activity in the balanced state originates from chaotic network dynamics (Vogels et al., 2005; Wallace et al., 2013; Ostojic, 2014). Albeit not specifically addressing the source of irregular activity in a balanced state, studies exist that have drawn the opposite conclusions in some special situations. For example, numerical simulations of neural networks consisting of pulse-coupled spiking neurons of only inhibitory type can display irregular activity in a dynamical state with a negative Lyapunov exponent (Zillmer et al., 2006; Jahnke et al., 2008; Monteforte and Wolf, 2012). Meanwhile, in the limit of fast synaptic response, any generic trajectory was shown to be asymptotically stable in inhibition-dominated networks (only a small fraction of connections can be excitatory) with the inhomogeneous delay distribution and strong coupling (Jahnke et al., 2009). Therefore, the question of whether the irregularity in a balanced state arises from chaos in a neuronal network with both excitatory and inhibitory neurons under more realistic Poisson drive, remains an important issue to be further clarified.

In this work, we first numerically show that, over a broad range of parameters, the balanced state can exist in current-based integrate-and-fire (I&F) neuronal networks consisting of both excitatory and inhibitory neurons with delta-pulse coupling currents and pulse-like external inputs. We then mathematically prove that, driven by any point process in time—not limited to Poisson point processes, the current-based I&F neuronal networks with delta-pulse interactions cannot exhibit chaotic dynamics. In fact, two nearby trajectories of such a network generically coalesce after a finite time regardless of whether the dynamics is in a balanced state or not. Our results demonstrate that in the delta-pulse coupled, current-based I&F system the irregular activity of the balanced state is not a consequence of a chaotic dynamical state. Our proof remains valid when the system possesses only excitatory or inhibitory population. This conclusion extends the previous results (Jin, 2002; Zillmer et al., 2006; Jahnke et al., 2009) that dynamics of the strongly inhibition-dominated networks are stable. By our analysis, stable dynamics with the irregular firing activity can occur in the delta-pulse coupled, current-based I&F neuronal network of any size.

2. Materials and Methods

2.1. I&F Model

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

dvikdt=-gL(vik-ϵRk)+Iik(t),    (1)

where gL 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,

Iik(t)=IikE(t)+IikI(t),    (2)

consists of two terms, where IikE(t)=fksδ(t-ζisk)+JkEj=1NECijkEsδ(t-τjsE) is the excitatory input, IikI(t)=-JkIj=1NICijkIsδ(t-τjsI) is the inhibitory input, δ(·) is the Dirac delta function, fk is the strength of the external input, and Jkl 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 fk. 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 Jkl 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 Nl is the total number of the lth population, k, l = E, I. The parameter values in our simulations are as follows: NE = 32000, NI = 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 = ν0K, where ν0 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 fE for the excitatory population and fI 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.

2.2. Analysis of the I&F Network

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 mE and mI 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 KmE (KmI), where mE (mI) 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.

FIGURE 1
www.frontiersin.org

Figure 1. Defining characteristics of a balanced network. (A) Balanced excitatory and inhibitory inputs into a sample I&F neuron. The amplitude of excitatory (red) and inhibitory (blue) inputs (normalized by gL) deviate significantly from the firing threshold (green solid line) whereas the summed total input (scaled by gL) (black) fluctuates around the threshold while keeping its mean (cyan) below the threshold; (B) The average fano-factor across the network as a function of the time bin size (black line). The values of fano-factor for different time bin sizes (larger than 100 ms) are above unity and nearly constant. The inset is the distribution of the fano-factor calculated from each neurons spike events over the entire network with the time bin size 400 ms. All the values are above unity, again exemplifying the irregularity of the spiking activity; (C) Distribution of neuronal firing rates normalized by the entire population-averaged firing rate. The distribution is broad and skewed as a consequence of the heterogeneity of the single neuron activity; (D) The distribution of cross-correlation coefficient of spike events for each pair of neurons in the network. It has a sharp peak around zero, therefore, the system is nearly uncorrelated. We choose the bin size Δt = 2ms to calculate the cross-correlation; (E) The upper panel is a raster plot of 500 excitatory (red dots) and 500 inhibitory (blue dots) neurons; the lower panel shows the number of the spiking neurons in the excitatory population (red) and inhibitory population (blue) in each time window. The raster plot shows there is no synchrony in the neuronal dynamics. The lower panel exhibits stationarity of the number of firing neurons in the network; (F) The mean firing rate of the excitatory and inhibitory population as a linear function of the external driving rate. For comparison, we vary the value of K here: K = 100 (dashed-dot), K = 400 (solid lines), and K = 3600 (dashed lines); red for excitatory population whereas blue for inhibitory population. The mean firing rates in the K → +∞ limit, mE=ν0 and mI=ν0, are also shown (black line, the theoretical gain curves for the excitatory and inhibitory populations overlap). Note that the minor disagreement between the theoretical prediction and the simulated gain curve arises from the finite size effect, and the gain curve converges to the infinite-K limit as K increases. The parameters are NE = 32000, NI = 8000, JEE=JIE=1.0/K, JII=1.8/K, JEI=2.0/K, fE=1.0/K, fI=0.8/K, νE = νI = ν0K and ν0 = 30Hz. In (A–E), K = 400.

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 Mk,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 Mk,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 vk,l(t) is given by a sum of independent identically distributed random variables for given Mk,l, the average neuronal voltage at time t can be simply expressed as

uk(t)=ϵRk+1-e-gLtgL(fkνk+KJkEmE-KJkImI)    (3)

for k = E, I. As discussed above, ϵRk=0, fk, JkE, and JkI are of order 1/K; νk is of order K; and mE and mI are of order 1. Then, the leading order of uk(t) is K. Obviously, the membrane potential cannot become infinite as K → +∞. Therefore, the leading order of K should vanish and one can obtain

mE=JIIfE-JEIfIJIEJEI-JEEJIIν0,mI=JIEfE-JEEfIJIEJEI-JEEJIIν0    (4)

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).

3. Results

3.1. Existence of Balanced States

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 JIE = JEE. Therefore, from Equation (4), requiring firing rates to be non-negative implies

fEfI>JEIJII>1    (5)

or

fEfI<JEIJII<1.    (6)

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 mE = 0 — that is, only inhibitory neurons fire in the system, the mean firing rate of the inhibitory population is then given by mI = fIν0/JII. 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, fE, fI, JEI, and JII 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 fE/fI and JEI/JII that we have scanned. Each dot in the figure can represent a set of parameters fE, fI, JEI, and JII with fixed ratios of fE/fI and JEI/JII. 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.

FIGURE 2
www.frontiersin.org

Figure 2. Dynamical network state as a function of fE/fI and JEI/JII. The range of parameters here is as follows: NE = 32000, NI = 8000, K = 400, 0 < fE < 0.2, 0.04 < fI < 1.0, νEfE = 10, νI = νE, JEE=JIE=1.0/K, 1.0/K<JEI<5.0/K, 0.0<JII<10.0/K. We choose values of the four parameters fE, fI, JEI and JII at random to verify whether the state of the network with these parameters is a balanced state. Note that each dot (C1, C2) in the panel (A) can represent a group of different parameter values as long as they satisfy fE/fI=C1 and JEI/JII=C2, where C1 and C2 are constant. A red dot represents a balanced state in which both excitatory and inhibitory neurons fire, with a sample raster in panel (B) (which has the parameters marked by the black asterisk in the area of red dots, and only 1000 excitatory and 1000 inhibitory neurons are shown here). A green dot represents a state in which only the inhibitory population is firing. This dynamical state also exhibits all the defining characteristics of a balanced state for the inhibitory population, while the excitatory population fires few spikes over a long time. A corresponding sample raster is shown in the panel (C) (which has the parameters marked by the black asterisk in the area of green dots, and only 1000 excitatory and 1000 inhibitory neurons are shown here). Finally, a blue dot represents a state of strong synchrony, for which a sample raster is shown in the panel (D) (which has the parameters marked by the black asterisk in the area of blue dots, and only 1000 excitatory and 1000 inhibitory neurons are shown here). The black dashed lines in the panel (A) are y = x, y = 1 and x = 1. In the raster plot, red dots denote the spike times of the excitatory population and blue dots denote the spike times of the inhibitory population.

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 fE and JII in comparison to fI and JEI, 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 mE = 0 and mI = fIν0/JII into Equation (3), and requiring that uE becomes sufficiently negative as K becomes large so the excitatory population rarely fires, we obtain the condition

fEfI<JEIJII,    (7)

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.

3.2. Absence of Chaos

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 = t0, i.e., ϵ = |(t0) − v(t0)| ≪ 1. The dynamics of the perturbed trajectory (t) is described by the equation

dikdt=-gL(ik-ϵRk)+fksδ(t-ζisk)        +JkEj=1NECijkEsδ(t-τ~jsE)-JkIj=1NICijkIsδ(t-τ~jsI),    (8)

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,

λmax=limTlimϵ01Tln (|(T)-v(T)|ϵ),    (9)

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.

3.2.1. Spike Train Sorting

For the preparation of the proof of the absence of chaos, we label all the neurons as {1, 2, …, NE + NI}, in which {1, 2, …, NE} labels the excitatory neurons and {NE + 1, NE + 2, …, NE + NI} labels the inhibitory neurons for easy description. That is τjs (τ~js) is the sth spike of the jth neuron, where j = 1, 2, …, NE stands for an excitatory neuron and j = NE + 1, NE + 2, …, NE + NI stands for an inhibitory neuron.

Then for any fixed finite time T, we sort the spike times of the two trajectories {τpq} and {τ~pq} into the increasing lists. Recall the fact that the voltage of both the reference and the perturbed trajectories of any neuron will cross the threshold only upon receiving spikes either from the external or excitatory recurrent input. Hence, the neurons fire precisely at the arrival time of the external or excitatory recurrent input spikes. There may be a group of neurons that fire at a particular time but with a distinct firing sequence amongst these neurons (see our firing dynamics described in the event driven algorithm in section Materials and Methods). Note that these simultaneously firing neurons only spike at the time when an external spike is received by some of the neurons within this group (possibly by all the neurons in the group) for the type of I&F neuronal networks with delta-pulse coupling. When we meet simultaneous firings during the sorting process, suppose these simultaneous firings are τp1q1, τp2q2, ···, τpaqa such that τp1q1 = τp2q2 = ··· = τpaqa, then we perform the following strategy to make the spike timing list unique.

1. p1 is chosen to be the smallest neuron label among the neurons that receive an external input spike in the simultaneously-firing group;

2. the remaining sequence of τp1q1, τp2q2, ···, τpaqa is sorted and reordered such that p2 < p3 < ··· < pa.

The same rule is applied in the perturbed network. Then we can obtain two unique increasing spike time sequences

τp1q1τp2q2···τpMqM,τ~p~1q~1τ~p~2q~2···τ~p~M~q~M~,    (10)

where τprqr (τ~p~r~q~r~) denotes the qr(q~r~)th spike of the pr(p~r~)th neuron in the reference (perturbed) trajectory, pr(p~r~)=1,2,,NE stands for an excitatory neuron and pr(q~r~)=NE+1, NE + 2, …, NE + NI stands for an inhibitory neuron. M and M~ are the total number of spikes in the reference and perturbed trajectories, respectively.

3.2.2. Proposition and Proof

We next turn to the following proposition,

PROPOSITION For each reference voltage trajectory v(t)=(v1(t),v2(t),,vNE(t),vNE+1(t),vNE+2(t) ,,vNE+NI(t))(v1E(t),v2E(t),,vNEE(t),v1I(t), v2I(t),,vNII(t)) described by Equation (1), and its perturbed voltage trajectory (t)=(1(t),2(t),,NE(t),NE+1(t), NE+2(t),,NE+NI(t))(1E(t), 2E(t),,NEE(t),1I(t), 2I(t), ,NII(t)) described by Equation (8), with their initial difference ϵ = |(t0) − v(t0)| at time t = t0, one can obtain the increasing spiking lists (10) for both trajectories over any fixed finite time T according to the sorting process described in the above. If the initial perturbation is sufficiently small, i.e., ϵ ≪ 1, then M~=M, τ~p~rq~r=τprqr, p~r=pr and q~r=qr for any r in the lists (10).

Note that, before the first spike of reference and perturbed trajectories, the time evolution of the perturbation δv(t)=v~(t)-v(t) can be obtained from the system of equations

ddtδvi=-gLδvi(t)  for i=1,2,,NE+NI.    (11)

That is, for a sufficiently small initial perturbation size ϵ, the distance between v(t) and (t) decays exponentially. Then we use mathematical induction to prove this proposition:

1. For the first firing event r = 1, either τp1q1 or τ~p~1q~1 is equal to some external input spike time (the first spike must be induced by the external input). Because the reference and perturbed trajectories receive the same external input, for a sufficiently small initial perturbation size ϵ, we have τ~p~1q~1=τp1q1, p~1=p1 and q~1=q1;

2. Suppose we have τ~p~rq~r=τprqr, p~r=pr and q~r=qr for rm, we then show that τ~p~m+1q~m+1=τpm+1qm+1, p~m+1=pm+1 and q~m+1=qm+1;

3. In fact, before the (m+1)th spike of both the reference and perturbed trajectories, we can find that δv(t) is still governed by Equation (11) because they receive the same external and recurrent input spikes. For a sufficiently small initial perturbation size ϵ, Equation (11) also ensures p~m+1=pm+1 and q~m+1=qm+1. Then, if the qm+1th spike of the pm+1th neuron is caused by the arrival of an excitatory recurrent input spike, this spike can always be traced back to the spike that is caused directly by an external input spike. By the ordering rule of neurons in the simultaneously firing group above, we have τpm+1qm+1 = τpsqs and τ~p~m+1q~m+1=τ~p~sq~s, where either the spike of τpsqs or the spike of τ~p~sq~s is caused by an external spike to the psth neuron (sm). Consequently, τ~p~m+1q~m+1=τpm+1qm+1. If the qm+1th spike of the pm+1th neuron is caused directly by the arrival of an external spike, clearly we have τ~p~m+1q~m+1=τpm+1qm+1 because the reference and perturbed trajectories receive the same external input spike;

4. Finally, we can readily obtain M~=M.

Therefore, τ~p~rq~r=τprqr holds for r = 1, 2, …, M with p~r=pr and q~r=qr. Next, according to the Equation (9), it is obvious that the largest Lyapunov exponent is always negative for the dynamics described by Equation (11). That is, the I&F network with delta-pulse coupling exhibits no chaotic dynamics while receiving pulse-like external input. Note that the voltage of both the reference and perturbed trajectories of any neuron, for example, the ith neuron in the kth population, will be reset to ϵRk after its first firing. Because both the external input spike time and the synaptic input spike time are the same in the reference and perturbed networks, the two trajectories will no longer separate from each other. This means that the difference between the reference and perturbed trajectories of the ith neuron will converge after its first firing event. Thus, the largest Lyapunov exponent will approach negative infinity. As show in Figure 3, the numerical results for the given fixed network size N and mean degree connectivity K indeed show that the total perturbation of all neurons' voltages always exponentially decay with time and eventually converge to zero in a finite time.

FIGURE 3
www.frontiersin.org

Figure 3. Time evolution of the total perturbation of all neurons' voltages. The total perturbation of all neurons' voltages |δv(t)| in either region described in Figure 2A always exponentially decay with time and eventually converge to zero in a finite time. (A) The parameters are marked by the black asterisk in the area of red dots in Figure 2A; (B) The parameters are marked by the black asterisk in the area of green dots in Figure 2A; (C) The parameters are marked by the black asterisk in the area of blue dots in Figure 2A. The parameters are chosen to make the mean firing rate of spiking neurons around 30 Hz. The total perturbation of all neurons' voltages at initial time in each case is chosen to be |δv(0)| = 5 × 10−4. The slope of the solid part of the curve in each case is −gL. The vertical dashed part of each curve indicates that the difference between the reference and perturbed trajectories vanishes, namely, the largest Lyapunov exponent approaches negative infinity.

As a matter of fact, the proof can be extended to a more generalized system. For example, the membrane potential vik of the ith neuron in the kth population obeys equation with refractory periods

dvikdt=-giLk(vik-ϵiRk)+fksδ(t-ζisk)+JkEj=1NECijkEs        δ(t-τjsE-τiE)-JkIj=1NICijkIsδ(t-τjsI-τiI),    (12)

but the neuron will come into a refractory period after vik crosses the threshold ϵiTk. That is when vik(τisk)ϵiTk, one has vik(t)=ϵiRk for τisk<tτisk+τi,refk, where τi,refk is the refractory period of the ith neuron in the kth population. τiE and τiI are the excitatory and the inhibitory synaptic delay of the ith neuron, respectively. In this system, the leakage conductance, the threshold, the resting potential, the refractory period and the synaptic delay of one neuron can be different from that of another neuron. Note that a neuron in this system can generate a spike also only upon receiving spikes either from the external or excitatory recurrent input, and the neuron fires precisely either at the arrival time of the external spike or at the arrival time of the excitatory recurrent spike with a fixed synaptic delay. Therefore, we can extend the proof above directly here to conclude that there is no chaotic dynamics in this generalized system. Despite the demonstration of irregular firing of neurons in the balanced state of this current-based I&F network as above, the irregular activity is not induced by chaos.

4. Discussion

Cortical neurons often exhibit spiking dynamics that are highly irregular. It is believed that the irregular neural spiking activity can be generated from a balance between excitatory and inhibitory inputs to a neuron. In the early work (Vreeswijk and Sompolinsky, 1998), the neuronal network consisting of binary neurons was applied to explain the balanced neuronal state. Meanwhile, the dynamics of the binary neuronal system was found to be chaotic in the balanced state. Therefore, chaos was then often thought to be closely related to the irregular firing activity in a balanced network. Note that the binary neuron is a highly idealized model for describing a neuron.

Here we address the issue of chaotic origin of spiking irregularity using a current-based leaky integrate-and-fire (I&F) model with delta-pulse coupling and pulse-like external input. By examining the defining characteristics of a balanced state (Vreeswijk and Sompolinsky, 1998), and by exploring a wide range of parameter values, we have found that the balanced state persists robustly in this I&F system.

We then apply the analysis of the largest Lyapunov exponent to characterize the dynamics of this I&F system. We mathematically demonstrate that the largest Lyapunov exponent is always negative in the current-based I&F system with delta-pulse coupling and pulse-like external input, in which each neuron can have its own distinct values of leakage conductance, resting potential, spiking threshold, refractory period, and synaptic decay. The reference and perturbed realizations of a neuronal network trajectory can converge in finite time. It is worthwhile to point out that the non-chaotic property in our proof holds for neuronal networks of any size.

To understand under what conditions cortical-like firing irregularity can be generated, several efforts have been made. For instance, an early study (Renart et al., 2007) investigated the possible existence of multiple balanced steady states with persistent activity. It shows that the Poisson-like irregular spiking activity can arise from balanced regimes with sustained persistent activity. The neuronal network system used in the study (Renart et al., 2007) is the same as that in our work, but with different scaling of synaptic connection strengths, e.g., homogeneous or heterogeneous multicolumnar architecture. Note that the proof of non-chaotic mechanism in our system is independent of scaling structure of connection strengths as long as the synaptic interaction is delta-pulse coupled. Therefore, it is expected that the non-chaotic mechanism also works in persistent activity states.

In addition, another recent study found that the I&F neuronal network endowed with probabilistic synaptic transmission can underlie the Poisson-like spiking variability over a wide range of firing rates (Moreno-Bote, 2014). Here, we show that the neuronal network consisting of I&F neurons with delta-pulse couplings is always non-chaotic over a wide range of firing rates. It could be interesting to investigate whether our non-chaotic neuronal network with probabilistic synaptic transmission can achieve Poisson-like variability of the spiking activity as observed in cortex. For example, we can consider that each pair of connected neurons in our system has multiple delta-pulse synaptic contacts, and each contact is activated with certain probability.

As is shown here in our works, the irregular spiking dynamics is always stable in the current-based I&F model with delta-pulse coupling and pulse-like external input. Therefore, the irregular activity cannot simply arise from the underlying chaotic dynamics. We point out that the proof of the non-chaotic dynamics in our system relies on the facts that both the external input and the synaptic interactions are in the form of delta-pulse coupling. Therefore, for any fixed network size N and the mean degree connectivity K, the neuron fires precisely at the arrival time of the external or excitatory recurrent spikes. As a result, the difference between the reference and perturbed trajectories of each neuron will converge after its first firing event. Thus, the largest Lyapunov exponent will approach negative infinity and the dynamics is non-chaotic. As for the case of smooth synaptic coupling, e.g., the current Iik(t) in Equation (2) is an α-like function (Dayan and Abbott, 2001) with the rise time constant τr and the decay time constant τd as Iik(t)=(e-t/τr-e-t/τd)/(τr-τd), it is expected that the non-chaotic mechanism still holds when the time constants τr and τd are both relatively small. However, if these time constants are not sufficiently small, e.g., greater than 2 ms, it is known that the network dynamics is not always stable and can be chaotic (Zhou et al., 2010; Harish and Hansel, 2015). In general, the phenomenon of chaos is model-dependent (Brette, 2004; Zhou et al., 2009; Sun et al., 2010) and could not be the ultimate source of irregularity in neuronal activity in the brain.

Author Contributions

QG, ZT, GK, DZ, and DC: Conceived and designed the research, Performed experiments and analyzed data, Wrote the paper.

Conflict of Interest Statement

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.

Acknowledgments

This work is supported by NYU Abu Dhabi Institute G1301 (QG, ZT, DZ, and DC); NSFC-11671259, NSFC-11722107, NSFC-91630208, and Shanghai Rising-Star Program-15QA1402600 (DZ); NSF DMS-1009575, NSFC-31571071 (DC); Shanghai 14JC1403800, 15JC1400104, and SJTU-UM Collaborative Research Program (DC and DZ). The remaining authors dedicate this paper to their late coauthor and mentor DC.

References

Braitenberg, V., and Schuz, A. (1998). “Comparison between synaptic and neuronal density,” in Cortex: Statistics and Geometry of Neuronal Connectivity (Berlin; Heidelberg: Springer), 37–38.

Google Scholar

Brette, R. (2004). Dynamics of one-dimensional spiking neuron models. J. Math. Biol. 48, 38–56. doi: 10.1007/s00285-003-0223-9

PubMed Abstract | CrossRef Full Text | Google Scholar

Brette, R., Rudolph, M., Carnevale, T., Hines, M., Beeman, D., Bower, J. M., et al. (2007). Simulation of networks of spiking neurons: a review of tools and strategies. J. Comput. Neurosci. 23, 349–398. doi: 10.1007/s10827-007-0038-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Britten, K. H., Shadlen, M. N., Newsome, W. T., and Movshon, J. A. (1993). Responses of neurons in macaque mt to stochastic motion signals. Vis. Neurosci. 10, 1157–1169. doi: 10.1017/S0952523800010269

PubMed Abstract | CrossRef Full Text | Google Scholar

Cinlar, E. (1972). “Superposition of point processes,” in Stochastic Point Processes: Statistical Analysis, Theory, and Applications, ed P. A. W. Lewis (New York, NY: John Wiley), 549–606.

Compte, A., Constantinidis, C., Tegnér, J., Raghavachari, S., Chafee, M. V., Goldman-Rakic, P. S., et al. (2003). Temporally irregular mnemonic persistent activity in prefrontal neurons of monkeys during a delayed response task. J. Neurophysiol. 90, 3441–3454. doi: 10.1152/jn.00949.2002

PubMed Abstract | CrossRef Full Text | Google Scholar

Dayan, P., and Abbott, L. F. (2001). Theoretical Neuroscience. Cambridge, MA: MIT Press.

Google Scholar

Gütig, R., and Sompolinsky, H. (2006). The tempotron: a neuron that learns spike timing–based decisions. Nat. Neurosci. 9, 420–428. doi: 10.1038/nn1643

PubMed Abstract | CrossRef Full Text | Google Scholar

Haider, B., Duque, A., Hasenstaub, A. R., and McCormick, D. A. (2006). Neocortical network activity in vivo is generated through a dynamic balance of excitation and inhibition. J. Neurosci. 26, 4535–4545. doi: 10.1523/JNEUROSCI.5297-05.2006

PubMed Abstract | CrossRef Full Text | Google Scholar

Harish, O., and Hansel, D. (2015). Asynchronous rate chaos in spiking neuronal circuits. PLoS Comput. Biol. 11:e1004266. doi: 10.1371/journal.pcbi.1004266

PubMed Abstract | CrossRef Full Text | Google Scholar

Hertz, J., and Prügel-Bennett, A. (1996). Learning short synfire chains by self-organization*. Network 7, 357–363. doi: 10.1088/0954-898X_7_2_017

PubMed Abstract | CrossRef Full Text | Google Scholar

Holmgren, C., Harkany, T., Svennenfors, B., and Zilberter, Y. (2003). Pyramidal cell communication within local networks in layer 2/3 of rat neocortex. J. Physiol. 551, 139–153. doi: 10.1113/jphysiol.2003.044784

PubMed Abstract | CrossRef Full Text | Google Scholar

Holt, G. R., Softky, W. R., Koch, C., and Douglas, R. J. (1996). Comparison of discharge variability in vitro and in vivo in cat visual cortex neurons. J. Neurophysiol. 75, 1806–1814. doi: 10.1152/jn.1996.75.5.1806

PubMed Abstract | CrossRef Full Text | Google Scholar

Jahnke, S., Memmesheimer, R.-M., and Timme, M. (2008). Stable irregular dynamics in complex neural networks. Phys. Rev. Lett. 100:048102. doi: 10.1103/PhysRevLett.100.048102

PubMed Abstract | CrossRef Full Text | Google Scholar

Jahnke, S., Memmesheimer, R.-M., and Timme, M. (2009). How chaotic is the balanced state? Front. Comput. Neurosci. 3:13. doi: 10.3389/neuro.10.013.2009

PubMed Abstract | CrossRef Full Text | Google Scholar

Jin, D. Z. (2002). Fast convergence of spike sequences to periodic patterns in recurrent networks. Phys. Rev. Lett. 89:208102. doi: 10.1103/PhysRevLett.89.208102

PubMed Abstract | CrossRef Full Text | Google Scholar

Litwin-Kumar, A., and Doiron, B. (2012). Slow dynamics and high variability in balanced cortical networks with clustered connections. Nat. Neurosci. 15, 1498–1505. doi: 10.1038/nn.3220

PubMed Abstract | CrossRef Full Text | Google Scholar

London, M., Roth, A., Beeren, L., Häusser, M., and Latham, P. E. (2010). Sensitivity to perturbations in vivo implies high noise and suggests rate coding in cortex. Nature 466, 123–127. doi: 10.1038/nature09086

PubMed Abstract | CrossRef Full Text | Google Scholar

Miura, K., Tsubo, Y., Okada, M., and Fukai, T. (2007). Balanced excitatory and inhibitory inputs to cortical neurons decouple firing irregularity from rate modulations. J. Neurosci. 27, 13802–13812. doi: 10.1523/JNEUROSCI.2452-07.2007

PubMed Abstract | CrossRef Full Text | Google Scholar

Monteforte, M., and Wolf, F. (2012). Dynamic flux tubes form reservoirs of stability in neuronal circuits. Phys. Rev. X 2:041007. doi: 10.1103/PhysRevX.2.041007

CrossRef Full Text | Google Scholar

Moreno-Bote, R. (2014). Poisson-like spiking in circuits with probabilistic synapses. PLoS Comput. Biol. 10:e1003522. doi: 10.1371/journal.pcbi.1003522

PubMed Abstract | CrossRef Full Text | Google Scholar

Newhall, K. A., Kovacic, G., Kramer, P. R., Zhou, D., Rangan, A. V., and Cai, D. (2010). Dynamics of current-based, poisson driven, integrate-and-fire neuronal networks. Commun. Math. Sci. 8, 541–600. doi: 10.4310/CMS.2010.v8.n2.a12

CrossRef Full Text | Google Scholar

Oseledec, V. I. (1968). A multiplicative ergodic theorem. lyapunov characteristic numbers for dynamical systems. Trans. Moscow Math. Soc. 19, 197–231.

Google Scholar

Ostojic, S. (2014). Two types of asynchronous activity in networks of excitatory and inhibitory spiking neurons. Nat. Neurosci. 17, 594–600. doi: 10.1038/nn.3658

PubMed Abstract | CrossRef Full Text | Google Scholar

Ott, E. (2002). Chaos in Dynamical Systems. Cambridge, UK: Cambridge University Press.

Parker, T. S., and Chua, L. (2012). Practical Numerical Algorithms for Chaotic Systems. Springer Science & Business Media.

Google Scholar

Peters, A. (1987). “Number of neurons and synapses in primary visual cortex,” in Cerebral Cortex (Springer), 267–294.

Google Scholar

Renart, A., De La Rocha, J., Bartho, P., Hollender, L., Parga, N., Reyes, A., et al. (2010). The asynchronous state in cortical circuits. Science 327, 587–590. doi: 10.1126/science.1179850

PubMed Abstract | CrossRef Full Text | Google Scholar

Renart, A., Moreno-Bote, R., Wang, X.-J., and Parga, N. (2007). Mean-driven and fluctuation-driven persistent activity in recurrent networks. Neural Comput. 19, 1–46. doi: 10.1162/neco.2007.19.1.1

PubMed Abstract | CrossRef Full Text | Google Scholar

Shadlen, M. N., and Newsome, W. T. (1998). The variable discharge of cortical neurons: implications for connectivity, computation, and information coding. J. Neurosci. 18, 3870–3896. doi: 10.1523/JNEUROSCI.18-10-03870.1998

PubMed Abstract | CrossRef Full Text | Google Scholar

Shu, Y., Hasenstaub, A., and McCormick, D. A. (2003). Turning on and off recurrent balanced cortical activity. Nature 423, 288–293. doi: 10.1038/nature01616

PubMed Abstract | CrossRef Full Text | Google Scholar

Sun, Y., Zhou, D., Rangan, A. V., and Cai, D. (2010). Pseudo-lyapunov exponents and predictability of hodgkin-huxley neuronal network dynamics. J. Comput. Neurosci. 28, 247–266. doi: 10.1007/s10827-009-0202-2

PubMed Abstract | CrossRef Full Text | Google Scholar

Sussillo, D., and Abbott, L. F. (2009). Generating coherent patterns of activity from chaotic neural networks. Neuron 63, 544–557. doi: 10.1016/j.neuron.2009.07.018

PubMed Abstract | CrossRef Full Text | Google Scholar

Troyer, T. W., and Miller, K. D. (1997). Physiological gain leads to high isi variability in a simple model of a cortical regular spiking cell. Neural Comput. 9, 971–983. doi: 10.1162/neco.1997.9.5.971

PubMed Abstract | CrossRef Full Text | Google Scholar

van Vreeswijk, C., and Sompolinsky, H. (1996). Chaos in neuronal networks with balanced excitatory and inhibitory activity. Science 274, 1724–1726. doi: 10.1126/science.274.5293.1724

PubMed Abstract | CrossRef Full Text | Google Scholar

Vogels, T. P., and Abbott, L. F. (2005). Signal propagation and logic gating in networks of integrate-and-fire neurons. J. Neurosci. 25, 10786–10795. doi: 10.1523/JNEUROSCI.3508-05.2005

PubMed Abstract | CrossRef Full Text

Vogels, T. P., Rajan, K., and Abbott, L. F. (2005). Neural network dynamics. Annu. Rev. Neurosci. 28, 357–376. doi: 10.1146/annurev.neuro.28.061604.135637

PubMed Abstract | CrossRef Full Text | Google Scholar

Vreeswijk, C. v., and Sompolinsky, H. (1998). Chaotic balanced state in a model of cortical circuits. Neural Comput. 10, 1321–1371. doi: 10.1162/089976698300017214

PubMed Abstract | CrossRef Full Text | Google Scholar

Wallace, E., Maei, H. R., and Latham, P. E. (2013). Randomly connected networks have short temporal memory. Neural Comput. 25, 1408–1439. doi: 10.1162/NECO_a_00449

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhou, D., Rangan, A. V., Sun, Y., and Cai, D. (2009). Network-induced chaos in integrate-and-fire neuronal ensembles. Phys. Rev. E 80:031918. doi: 10.1103/PhysRevE.80.031918

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhou, D., Sun, Y., Rangan, A. V., and Cai, D. (2010). Spectrum of lyapunov exponents of non-smooth dynamical systems of integrate-and-fire type. J. Comput. Neurosci. 28, 229–245. doi: 10.1007/s10827-009-0201-3

PubMed Abstract | CrossRef Full Text | Google Scholar

Zillmer, R., Livi, R., Politi, A., and Torcini, A. (2006). Desynchronization in diluted neural networks. Phys. Rev. E 74:036203. doi: 10.1103/PhysRevE.74.036203

PubMed Abstract | CrossRef Full Text | Google Scholar

Keywords: balanced state, irregular activity, chaotic dynamics, delta-pulse coupling, largest Lyapunov exponent

Citation: Gu QL, Tian ZK, Kovačič G, Zhou D and Cai D (2018) The Dynamics of Balanced Spiking Neuronal Networks Under Poisson Drive Is Not Chaotic. Front. Comput. Neurosci. 12:47. doi: 10.3389/fncom.2018.00047

Received: 17 August 2017; Accepted: 30 May 2018;
Published: 28 June 2018.

Edited by:

David Hansel, Université Paris Descartes, France

Reviewed by:

Ruben Moreno-Bote, Universidad Pompeu Fabra, Spain
Marc Timme, Technische Universität Dresden, Germany
Carl Van Vreeswijk, Centre National de la Recherche Scientifique (CNRS), France

Copyright © 2018 Gu, Tian, Kovačič, Zhou and Cai. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

*Correspondence: Douglas Zhou, zdz@sjtu.edu.cn

Disclaimer: All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article or claim that may be made by its manufacturer is not guaranteed or endorsed by the publisher.