Collective irregular dynamics in balanced networks of leaky integrate-and-fire neurons

We extensively explore networks of weakly unbalanced, leaky integrate-and-fire (LIF) neurons for different coupling strength, connectivity, and by varying the degree of refractoriness, as well as the delay in the spike transmission. We find that the neural network does not only exhibit a microscopic (single-neuron) stochastic-like evolution, but also a collective irregular dynamics (CID). Our analysis is based on the computation of a suitable order parameter, typically used to characterize synchronization phenomena and on a detailed scaling analysis (i.e. simulations of different network sizes). As a result, we can conclude that CID is a true thermodynamic phase, intrinsically different from the standard asynchronous regime.


I. INTRODUCTION
The use of simple models proves often very helpful to identify and characterize the mechanisms underlying general dynamical phenomena. Computational neuroscience is a field where this approach is potentially very powerful, given the myriad of interactions involved in the functioning of the mammalian brain [1,2]. However, setting the appropriate level of simplicity is not a priori obvious. A particularly enlightening example is the reproduction of the background neural activity. Most of the numerical and theoretical studies are based on the so-called rate models, where each neuron is characterized by a single coarse-grained variable representing the strength of the ongoing activity [3,4]. However, it is well known that neurons work by emitting single spikes, so that it is more natural to represent them as (nonlinear) oscillators. This is, indeed the philosophy adopted by many studies based on pulse coupled units, such as leaky integrate-and-fire (LIF) neurons [5]. Accordingly, a general question arises as to whether the two approaches are consistent with one another and, in particular, to what extent spiking neurons reproduce the scenario observed in rate models [6,7].
In this paper we look at the evolution of the so-called balanced networks, where excitatory and inhibitory interactions compensate each other [8], since the singleneuron dynamics is rather irregular and reminiscent of the background neural activity. A detailed theory of this regime has been developed for rate models in the limit of large connectivity [9][10][11][12], but some theoretical studies have been also made for spiking LIF neurons in highly diluited networks [13]. Altogether, the scenario which emerges from these studies is that of an asynchronous regime, i.e. a dynamical state characterized by "microscopic" fluctuations, but a steady and constant firing rate at the macroscopic level (in the thermodynamic limit).
Nevertheless, evidence of irregular collective dynamics has been recently found in both rate models [14] and spiking neurons [15], so that a question arises about the conditions for the emergence of partial synchronization.
In this paper we focus on networks of LIF spiking neurons, as we believe that this is a more realistic and meaningful setup. We extend the preliminary analysis presented in [15], by including a study of the avalanches and of the scaling behavior of the fluctuations of the spiking activity. Furthermore, we explore several variants of the model to test the robustness of the results of our claims. More precisely, we refer to a model that has been repeatedly investigated in the literature and in particular in [13,16] (several parameter values are herein selected as in Ref. [16]). However, there are important differences that should be stressed and, which allow us drawing convincing conclusions about the existence of a CID: (i) instead of sparse networks, we consider massive ones (i.e. the connectivity K is assumed to be proportional to the network size N ); (ii) the coupling strength is chosen to be of the order of 1/ √ K. Furthermore, (iii) the unbalance is chosen to be of the same order as the statistical fluctuations of the input current, so as to avoid a complete dominance of either excitation, or inhibition. Equipped with such assumptions, we have simulated different network sizes, finding that the seemingly CID observed for finite networks survives in the thermodynamic limit, thereby validating its identification as of a true thermodynamic phase.
In section II we introduce the model, briefly review the numerical scheme adopted for its simulations, and discuss a problem related to its ill-defined structure. In fact, in this model strictly simultaneous events can unavoidably occur, which require an additional protocol to specify the way they should be treated. This phenomenon is related to the occurrence of avalanches which we show to marginally influence the thermodynamic limit. In sec-tion III, we deal with the dynamical properties of the single-neuron dynamics, mostly focusing on the statistics of the inter-spike intervals (ISI) and on the spectral properties of the spike trains emitted by a single neuron as well as of the post-synaptic input currents received by each neuron. In section IV, we discuss the collective dynamics, introducing a suitable order parameter to quantify the degree of synchronization, and characterizing its stochastic-like behavior by means of the power spectrum of the global neural activity. The perturbative approach developed by Brunel in [13] is implemented to perform a comparison with the numerical observations. A qualitative agreement is found.
Finally a fractal-dimension analysis is performed, which confirms the stochastic-like character of the dynamics: i.e. its high-dimensional nature. In section V, we explore the collective behaviour of the LIF model when some of the parameters are varied, notably, delay, refractoriness, connectivity and the presence of external noise. All of the simulations confirm the robustness of CID. Section VI contains a summary of the main results and a brief list of the main open problems.

II. THE MODEL
In this Section we define the model following Ref. [13,16], equipped with a suitable scaling of some parameters, in order to preserve the CID observed at finite sizes. Furthermore, we discuss some intrinsic ambiguities present in the definition of the model that are related to the unavoidable occurrence of strictly synchronous events (the simultaneous arrival of spikes emitted by different presynaptic neurons crossing the threshold at the same moment). The treatment of these events is somehow arbitrary and requires the definition of a specific protocol. Furthermore, these synchronous events propagate in the network as a sequence of events separated by the propagation-time of the single spikes. The resulting avalanches are analysed in detail for networks of increasing size.
We consider networks of spiking neurons composed of N supra-threshold LIF neurons, split in bN excitatory cells and in (1 − b)N inhibitory ones [13,16]. The membrane potential V i of the generic i-th neuron evolves according to where τ = 20 ms is the membrane time constant, RI 0 = 24 mV is a constant external DC "current", and RI i corresponds to the synaptic current due to the recurrent connections within the network, namely where j(n) is the label of the node firing the n-th spike. The synaptic connections among the neurons are random without autapses, but with a fixed in-degree K for each neuron. In particular, we consider a massive network, where the in-degree grows proportionally to the size as K = cN ; the proportionality parameter c is termed connectivity -unless stated otherwise we set c = 0.1 throughout the entire paper. The random connections are embodied in the matrix G, whose elements take the values G ij = J e (−J i ), if the pre-synaptic neuron j is excitatory (inhibitory), otherwise G ij = 0. Whenever at time t j(n) the membrane potential V j of the j-th neuron reaches the threshold V th = 20 mV for the n-th time, two events are triggered: (i) the membrane potential is reset to V r = 10 mV and it is then held fixed for a refractory period τ r = 0.5 ms; (ii) a spike is emitted and received τ d = 0.55 ms later by the post-synaptic cells connected to neuron j.
In order to maintain a fixed balance between excitation and inhibition irrespective of the in-degree (and, thereby, of the system size) we assume that the coupling strength scales as the inverse of the square root of the in-degree, as done in most of the literature on the balanced state [9][10][11][12]. More precisely, we assume J e = J 1000/K and J i = (4 + g 1 c/K)J e [15]. With these choices, and for g 1 = 100, we recover the setup studied in [16] for N = 10, 000, K = 1000 and b = 0.8. The coupling strength J is our main control parameter.

A. Integration scheme
The model equations (1,2) can be solved by either implementing an event driven integration scheme, such as described in [17,18], or a more standard clock-driven strategy [19]. The former scheme is, in principle, exact; provided that the integration time step is small enough, also the latter scheme is sufficiently accurate. In fact, the results discussed in this paper have been obtained by implementing either scheme without any specific preference. One exception is the above mentioned presence of ambiguities in the very definition of the model, which can be singled out only with reference to the exact eventdriven scheme. In general, two types of event break the smooth evolution of the membrane potential: (i) a neuron reaches the threshold; (ii) a neuron receives one (or more) post-synaptic potentials (PSPs), elicited by one or more pre-synaptic neurons τ d = 0.55 ms earlier. In between these events all neurons evolve as being uncoupled according Eq. (1) with I i = 0. In order to evolve the system, we first need to identify the next type of event and thereby update all the membrane potentials {V i } until the occurrence of the event itself. This is done by solving analytically Eq. (1) with I i = 0. Thereafter, we process either the threshold passing (i) or the spike-receiving (ii) event. Each emitted spike is received τ d = 0.55 ms later by K = cN other neurons. If the sending neuron j is excitatory, it may happen that the post-synaptic potentials triggers several threshold passings events at exactly the same time. Since the delay is the same for all synaptic connections, in the next round, more than one of the simultaneous spikes can reach the same neuron. It is, therefore, necessary to complement the definition of the model with a rule to handle perfectly synchronous spikes, since the outcome depends on the way such events are treated.
We have decided that the most "neutral" rule consists in first estimating the net effect of all the PSPs on the basis of the current value V i of the membrane potential. Afterwards, all neurons which passed threshold are reset to the common value V r and the emitted spikes are stored to be later received by the postsynaptic neurons connected to the firing ones. We have verified that also in the case of event driven integration schemes spike avalanches separated by exactly one time delay can occur. In the next section we discuss how to quantify their relevance.

B. Characterization of the avalanches
We first monitored the number of simultaneously emitted spikes E and the corresponding density R E per unit of time; these are shown in Fig. 1 for different system sizes and two different synaptic coupling values. The density R E obviously increases with the system size. As shown in the insets of Fig. 1, all curves collapse on the same one, once the abscissa has been rescaled as E/ √ N . This means that the number of simultaneous events E is of order O( √ N ), while the total number of emitted spikes is of order O(N ). Therefore, the strictly simultaneous events become less and less relevant, while approaching the thermodynamic limit.
The presence of avalanches is a consequence of instantaneous synapses and identical time-delays. Avalanches arise when a spike triggers a cascade of sub-sequent spikes occurring at times exactly separated by the time delay. We have monitored the time duration (called length L) of the avalanches as well as their consistency (size S), corresponding to the total number of spikes emitted during an avalanche. The densities R L and R S per unit of time of the avalanche length L and size S are reported in Fig. 2 for various system sizes. There is a clear increase of the length and size of the avalanches with the system size N and hence (due to the massive coupling) with the in-degree K = cN . Upon rescaling the R S densities by √ N , they collapse onto a same curve (that we expect to depend on the coupling strength J) as shown in the inset of Fig. 2 (b). Finally, no simple scaling behavior has been found for the avalanche length L. We tested different scalings assumptions but none of them yielded a convincing data collapse. We can nevertheless safely conclude that the length grows even slower than 1 10 100 E 10 the size with N . Altogether, our results show that the avalanches unavoidably appear also in the exact event driven approach but do not contribute significantly to the network dynamics in the thermodynamic limit.

III. THE MICROSCOPIC DYNAMICS
The CID is a macroscopically observable phenomenon originated by an orchestrated interplay of the microscopic oscillators. Before introducing appropriate indicators to characterize the collective phenomena (see the following Section IV), we first shed light on the microscopic dynamics. Each oscillator i is characterised by a membrane potential V i which evolves continuously in time, but is affected by discrete events associated with the emission times {t i(n) } of the single spikes and the corresponding arrival times. A time trace of the individual membrane potential V 3 is shown in Fig. 3(b) as a red line together with the mean potential V (t) = 1/N i V i (t) (black line) (here and in the following the symbol · stands for an ensemble average). The membrane potential of a single neuron exhibits significantly larger fluctuations than those exhibited by the mean value V with only a limited correlation among the two observables (see Fig. 3(b)). The raster plot in the same time interval is depicted in Fig. 3(c); it clearly reveals irregular population bursts, whose degree of synchronization can be appreciated by looking at the global firing activity F g , i.e. the number of spikes emitted in a fixed time window per neuron (see Fig. 3 (d)). This last entity, whose time average corresponds to the average firing rate, reveals clear irregular oscillations.
As explained in Section II and shown in Fig. 1, simultaneous spike events are intrinsic properties of the model and can lead to the simultaneous arrival of a mixture of excitatory and inhibitory PSPs. Hence, the net result p(t) = RI i /τ can be either positive or negative. Additional information on the variability of the spiking activity is contained in the probability distribution function (PDF) of the ISIs of a single generic neuron. A few PDFs obtained for N = 10, 000 for different coupling strengths are reported in Fig. 5. For sufficiently large coupling, the PDF is composed of an exponential tail characteristic of a Poissonian dynamics plus a peak at very short ISI, which is the consequence of the occasionally periodic bursting activity of the neurons. These results clearly indicate that the single neuron dynamics is increasingly dominated by fluctuations for large J.
There are at least two ways to explore the spectral properties of the single neurons: by looking at the evolution of the membrane potential or by recording the spike events. We have focussed on the latter one, since it allows for a comparison between the input stimulus and the output response, as well as for the analysis of the collective spiking activity (discussed in the next Section).
In order to characterize the output signal we counted the number of emitted spikes within a fixed time window (we set it equal to 0.11 ms). The resulting signal is a series of 0s interspersed with a few 1s. The input is instead determined from the values of p(t) (reported in Fig. 4) coarse-grained over time bins of 0.11 ms. Furthermore, the power spectrum of the input signal has been rescaled according to the number of excitatory and inhibitory connections and their strength, i.e. by the factor cN (bJ 2 e + (b − 1)J 2 i ), to be comparable with the spectrum S S (f ) corresponding to the output time series.  At lower frequencies, and especially for 50 ≤ f ≤ 500 Hz, the S S (f ) spectra of input and output differ from one another (the differences persist for N = 160, 000, where they have reached an asymptotic shape -data not shown): in particular the input spectra exhibit a clear peak at 75 Hz, while the output ones reveal just a shoulder. In the presence of an asynchronous regime, input and output spectra should coincide (except for a scaling factor). In fact, in a series of recent papers, a recursive method was developed to generate asymptotic spectra, exactly by imposing a perfect correspondence between input and output [20,21]. The clear difference shown by our numerical results (Fig. 6) provides a first indication of a collective dynamics or, otherwise stated, of nontrivial correlations among the different neurons.
This will be extensively elucidated in the next Section, where we computed various indicators, including the power spectrum of the overall activity, whose spectrum is not too different from that of the input (see Fig. 9).

IV. COLLECTIVE DYNAMICS
In this section we discuss the collective dynamics which emerges from the correlations in the microscopic activity of the single neurons. A qualitative evidence is already noticeable in the structure of a typical raster plot, which consists in an irregular alternation of regions of different density (see Fig. 3(c)). The evolution of V provides a more accurate representation. The fluctuations of V are in fact smaller than those of the individual membrane potential V i (t) (see the red line in Fig. 3(b)), but nevertheless definitely appreciable.
On a quantitative level, it is convenient to introduce the synchronisation measure ρ where the overbar denotes a time average. Perfectly synchronised neurons behave in exactly the same way, so that the numerator and the denominator are equal to one another and ρ = 1. If instead, they are statistically independent, ρ ≈ 1/ √ N . The progression of the running average of ρ can be appreciated in Fig. 7 for J = 0.5 mV, see the middle bunch of trajectories, labeled by (0.5 , 0.55). The order parameter approaches ρ ≈ 0.35, irrespective of the networks size, indicating that CID is not a finite-size effect, but survives in the thermodynamic limit as defined in Ref. [15]. The jumps observed at early times of the running average are caused by sudden drops of the mean potential V . One of the strong drops is shown in Fig. 3(a) around t = 9, 500 ms. These rare events appear randomly at all times, but their effect on the cumulative average obviously decreases as time progresses.
The coefficient of variation C v is another measure of irregularity of the dynamics, based on the fluctuations of the ISI, rather than of the membrane potential. More The lower collection of lines corresponds to no delay τ d = 0.0 ms but with refractoriness τr = 0.5 ms. The system sizes are color coded in ascending order: black, red, green, blue and orange for N = 10, 000, 20, 000, 40, 000, 80, 000 and 160, 000, respectively.
precisely, we calculate where σ S is the standard deviation of the single-oscillator ISI, while τ S is the corresponding mean ISI. For J = 0 the single-neuron activity is strictly periodic and thus C v is equal to zero. We expect it to increase when the coupling strength J is switched on. For small J we can indeed appreciate a power-law growth, C v ≈ J α (see the black squares in Fig. 8) with a value of the rate α close to 1.6. The growth of C v continues for stronger coupling strenghts, becoming larger than 1, the value expected for a Poisson statistics. For sufficiently large coupling, we observe bursting dynamics of the neurons, corresponding to C v > 1. A similar behavior of C v with the coupling strength has been reported for inhibitory sparse networks of LIF in the absence of delay but for sufficiently slow synaptic decays [22]. Interestingly, the C v values are substantially independent of the system size (we have tested values of N up to 640, 000).
The coefficient of variation measures the amplitude of the fluctuations, but it is insensitive to temporal correlations: C v is strictly larger than zero already for a very regular sequence of ISIs such as a periodic alternation of two values t 1 and t 2 , with t 1 = t 2 . In order to have a more accurate indicator, we have computed the following diffusion coefficient. Let T n denote the time of n-th spike emitted by a given neuron, so that T n − T n−1 is the n-th ISI. Let, then be the diffusive coefficient of the process T n . We finally define Ξ is plotted in Fig. 8 (see the red crosses). Above J = 0.1, it basically coincides with C v , indicating that T n is essentially a renewal process. For J ≤ 0.1, Ξ decreases more slowly (with a rate close to 1.14) than C v . This is clearly due to the increasingly periodic character of the dynamics. The overall scenario is reminiscent of the phase transition discussed by Ostojic in [16]. The power spectrum of the global activity S g sheds light on collective phenomena from yet a different perspective. Analogously to the single oscillators, the spike times have been converted into a single time series counting the number of spikes emitted in each time bin. We chose the same time bin of 0.11 ms as in the previous cases. An example of the global field F g is included in Fig. 3 in the bottom panel (d); it clearly shows an irregular behavior. The power spectrum of the global activity has been divided by N 2 to allow for a meaningful comparison amongst different system sizes and with the single-neuron spectrum (Fig. 6). The spectra obtained for different system sizes are plotted in Fig. 9. For f > 40 Hz, they collapse onto one another, suggesting that the dynamics remains irregular in the thermodynamic limit, i.e. that the fluctuations exhibited by the collective variables are not finite-size effects. It should be stressed that in the absence of collective effects, i.e. for asynchronous states, the spectrum of the global activity would be proportional to N rather than to N 2 .
Below 40Hz, the spectral amplitude decreases with the system size, suggesting that the zero-frequency peak eventually disappears (at least for J = 0.5). Altogether, the spectral power is mostly concentrated in two frequency ranges: (i) a broad peak around f ≈ 75 Hz, which corresponds to the peak observed in the single neuron spectra S S and is presumably related to a time scale of the order of the membrane time constant τ = 20 msec and (ii) a peak around f ≈ 1818 Hz (and its multiples), which corresponds to the inverse of the delay. A comparison with the spectrum of the single neuron activity (see Fig. 6), reveals that the latter one is characterized by a much stronger high-frequency component (of white-noise type) and weaker peaks in correspondence of the inverse delay time.
We have finally implemented the perturbative approach developed by Brunel [13], based on the assumption of a sparse coupling. The idea basically consists in solving a self-consistent noisy Fokker-Planck equation for the probability density P (v, t) of the membrane potential v, The first term in the r.h.s. is nothing but the deterministic current defined in Eqs. (1,2), with µ e = RI 0 . The second, diffusive contribution, accounts for the unavoidable statistical fluctuations of the input signal arising from the coupling with the other neurons and its amplitude is estimated under the assumption of being a Poisson process. Finally, the last term is a common noise due to the fact that different neurons partially share the same input signals, whenever they share the same afferent neurons. More precisely, the drift µ is defined as where ν(t) is the instantaneous firing rate, while an expression for the diffusion coefficient σ 2 can be obtained by assuming that the spike train follows a Poisson statistics, Finally, σ 0 is the value of σ corresponding to the stationary value of the firing rate in the asynchronous regime ν 0 and J = J √ K (see the appendix for a more precise definition). The power spectrumn(ω) 2 of the neural activity can be determined by solving perturbatively the Fokker-Planck equation. The technical details are presented in Appendix A: we practically follow the method introduced in [13], the main difference being the numerical strategy adopted to determine the spectrum. The resulting curve is shown in Fig. 9 (black line) after converting the angular frequency ω into the frequency f . The perturbative approach qualitatively reproduces the shape of the spectrum, including the position of the peaks. On the other hand, the height of the peaks strongly deviates from the numerical simulations. For large coupling, the agreement worsens and the perturbative approach fails even in reproducing qualitatively the spectra at low frequencies, as shown in [15].
Finally, we have studied the neural activity by implementing nonlinear-dynamics tools, to determine the (effective) fractal dimension D e of the mean potential V (t). Given the sequence of values V (t n ), obtained by sampling the original signal every ∆t = 1 ms over 500 s (i.e. resulting in 500, 000 data points), this is embedded into a space of dimension m, by building vectors of the type [ V (t n ), V (t n+1 ), . . . , V (t n+m−1 )]. The fractal dimension has then been estimated by using a variant of the nearest-neighbour method recently proposed in [23]. In particular, N r reference points are randomly selected (N r = 10 5 in our case), then each reference point is compared with other n randomly selected points (up to the number of data points available) determining the distance ε m (k, n) of the k-th neighbour for different values of m and k. The distance is herein estimated using the maximum norm. An established theory [24], implies that for large n, where the angular brackets denote the average over the reference points, while D e is the information dimension. In practice, the logarithmic derivative of ε varies with n before eventually converging to its asymptotic value. Accordingly, it can be interpreted as an effective, resolutiondependent dimension, which is, in fact, independent of the order k of the neighbour considered in the simulations. In practice, given ε, the inverse of the logarithmic derivative is first determined and then plotted versus the resolution ε. The results are reported in Fig. 10. They show that, independently of the network size, the effective dimension increases upon decreasing the resolution. The stochastic-like nature of the dynamics is further confirmed.

V. ROBUSTNESS
In this section we investigate the robustness of CID, by testing its properties when some of the model parameters are modified, notably refractoriness, delay, connectivity and finally after introducing an external noise which acts independently on each neuron.
We start by showing the dependence of the average firing rate ν 0 and of the coefficient of variation C v on the system size in a setup where either delay or refractoriness is missing. From the data reported in the Table I, which refer to J = 0.5mV, we observe that, strange enough, the firing rate slows down, when the delay is removed. This is because the absence of delay induces a more homogeneous firing activity, which, in turn, is more dominated by the inhibitory neurons due to the weak unbalance. This interpretation is confirmed by the lower degree of synchronization that can be appreciated by looking at the (0.5,0.0) curves in Figure 7, which refer to different network sizes. Although ρ decreases, CID is still present in the absence of delay. In fact ρ is substantially indepen-TABLE I. Mean firing rates ν0 and mean coefficients of variation of the ISI Cv in absence of delay (τ d = 0) or no refractoriness (τr = 0) for different system sizes N , for J = 0.5 mV. The last two column reference to the standard setup with delay τ d = 0.55 ms and with refractoriness τr = 0.5 ms. dent of N (it actually even slowly increases). Additional studies of the spectral properties confirm that the collective dynamics is irregular (data not shown). This is at variance with the setup studied in [25], a heterogeneous ensemble of fully coupled, inhibitory, LIF neurons. In that context, CID disappears as soon as the delay vanishes. It is still to be understood whether the qualitative difference is due to the heterogeneity (dispersion in the bare firing rates of the single neurons).
Refractoriness is less relevant. From the data in Table I, we see that its absence does neither significantly modify the firing rate (which naturally increases by a small amount), nor the degree of irregularity of the single neuron. Even though C v slowly decreases upon increasing N , it remains substantially larger than 1, the expected value for a Poisson statistics. As for the collective dynamics, we notice in Fig. 7 (see the curves labeled (0.0,0.55)) that synchronization increases upon removing refractoriness. The convergence is slower than in the previous case: this is because of the presence of several sudden burst of synchronizations (see the upward jumps exhibited by ρ(t)), which require longer time scales for them to be suitably averaged out.
After having verified that neither delay nor refractoriness are necessary ingredient for CID to be observed, we now explore the role of the connectivity c. In Fig. 11 we plot three key parameters (the firing rate ν 0 , C v and ρ) as a function of the coupling strength J for different c values. In panel (a), we see that the firing rate is almost independent of c in the small coupling limit, while it progressively decreases upon increasing c in the strong coupling regime. This is due to the fact that a strong connectivity reduces the fluctuations which are known to be responsible for the larger ν 0 observed for strong coupling [26]. The progressive regularization of the neural activity is confirmed in panel (b), where we see that C v decreases upon increasing the connectivity.
Quite interesting is the dependence of the order parameter ρ on c. The clean data collapse for J < 0.4, indicates that up to a 30% connectivity, ρ scales as √ c, in agreement with the perturbative theory developed in [13]  and briefly recalled in Appx. A, which predicts a power spectrum proportional to c. The strong coupling regime (J > 0.4) seems to be characterized by different scaling properties, but additional simulations for different network sizes are required to put the statement on a more firm basis.
So far, our simulations have been performed for a slight prevalence of the inhibitory activity. In fact (for N = 10 4 ) the ratio between the two coupling strengths is g ≡ J i /J e = 5, to be compared with a 1 : 4 ratio of the two corresponding populations. In order to investigate the role of the degree of unbalance, we have studied two additional g-values, g = 4, and 5, which, respectively, correspond to a perfect balance and a stronger prevalence of inhibition. The results are presented in Fig. 12. The most important point is that CID is present for both parameter values, confirming the robustness of this phase. On a more quantitative level, unsurprisingly, the perfectly balanced state is characterized by a much stronger firing activity.
Finally, we analyse the role of noise, by adding iid white noise terms ξ(t) to the single neuron dynamics (Eq. (1)), such that ξ(t + τ )ξ(t) = 2Dδ(τ ). The dependence of ρ on D is reported in Fig. 13, for two different network sizes. The noise tends obviously to decrease the strength of the collective dynamics, without, however, killing it. In fact, CID survives even for moderately strong noise amplitudes, as it is appreciated by seeing that ρ does not vary significantly upon increasing N .
Altogether, CID is a very robust property, which survives even when noise is added, the connectivity is decreased, the balance is changed, the delay or refractoriness removed from the model equations. FIG. 12. The mean firing rate ν0, the mean coefficient of variation of the ISI Cv and the synchronisation measure ρ versus the coupling strength J for different balance factors g for a system size N = 10, 000. The dashed black line reference to the slight unbalance g = 5 used throughout the paper, whereas the solid green line shows the situation for stronger unbalance (g = 6) and the solid red line corresponds to the perfect balanced setup (g = 4). The blue crosses refer to a system size N = 40, 000 for a perfectly balanced situation, i.e. g = 4.

VI. CONCLUSIONS AND OPEN PROBLEMS
In this paper we have presented an extensive analysis of the collective dynamics emerging in a quasi-balanced network of LIF neurons. The irregularity of the collective dynamics is testified not only by the power spectra of the neural activity but also by a fractal-dimension analysis. The detailed simulations performed for different parameter values confirm that irregular dynamics is very ubiquitous. Several questions are, however, still open. Here we list the main ones.
(A) To what extent is this irregular dynamics related to the similar regime observed in globally coupled, heterogeneous neurons [23,25]? In those setups, which are reminiscent of the Kuramoto model, the heterogeneity seems to be a crucial ingredient, since CID disappears when the diversity among the neurons is removed. Here, it seems that the collective, stochastic-like dynamics is the result of a microscopic pseudo-chaotic evolution, which percolates up to macroscopic scales, as a consequence of the quasi-balanced regime. Whether this is really the correct explanation it is however still unclear.
(B) All the models so far explored assume δ-pulses, but this is obviously an approximation. The limit of infinitely narrow PSPs is singular, as shown, for instance, while investigating the stability of the splay state [27]. Furthermore, we have seen that the presence of strictly δlike pulses induces unavoidable synchronous events whose treatment requires additional ad-hoc hypotheses. It will therefore, be instructive to explore networks characterized by PSPs of finite duration, e.g by considering exponential or α-pulses.
(C) The numerical analysis has revealed that collective dynamics arises also for a very small coupling strength. The weak-coupling limit is typically amenable to a perturbative treatment. Accordingly, it is plausible that a model of Kuramoto-Daido phase-oscillators might be able to reproduce a similar regime and, at the same time, allow for an analytical treatment.
(D) The only limit where the irregular collective dynamics vanishes is that of a sparse network, where K/N → 0 for N → ∞. However, this statement refers to random Erdös-Rényi-type networks. It would be interesting to explore different more elaborated network structures as well as the role of heterogeneity.
(E) Qualitatively speaking, it looks like some differences exist between the weak and strong coupling regime. In the former case, the single neuron spiking activity is strongly correlated (being far from a renewal process) as shown in Fig. 8 and the strength of the collective dynamics is reproduced as expected by the perturbative theory (see the nice overlap among the curves reported in panel (c) of Fig. 11). In the latter case, the neuronal activity is very well approximated by a renewal process and, at the same time, the perturbative theory seems to fail already for a 1% connectivity. These differences suggest that at small coupling the neuronal dynamics is mean driven, i.e. is dominated by the mean value of the DC currents, while at large coupling it is fluctuation driven, i.e. the neurons are in proximity or below the threshold and the firings are triggered by fluctuations of the input currents. Similar transitions from mean to fluctuation driven dynamics has been recently reported in sparse inhibitory heterogeneous networks made of LIF neurons in [22] and composed of realistic models of striatal medium spiny neurons in [28]. An additional finite-size analysis is necessary to test whether this is a true transition that persists in the thermodynamic limit, as claimed by Ostojic [16] for strongly diluted networks.
(F) In all of our simulations, excitatory and inhibitory neurons have been assumed to be equal to one another. This implies that the same combination of excitatory and inhibitory fields is automatically consistent with the evolution of both types of neurons. This strong limitation should be lifted before drawing yet more general conclusions about the ubiquity of collective irregular dynamics.
So long as the fluctuating term in Eq. (6) can be neglected, the dynamics relaxes towards a stationary state which can be interpreted as an asynchronous regime characterized by a constant current ν 0 and constant fluctuations σ 0 , which can be determined self-consistently using the following expression for ν 0 , Vr −µ 0 −µe σ 0 du e u 2 (1 + erf(u)) (A3) It is now convenient to introduce the following changes of variables The threshold and reset potentials become, while the Fokker-Planck equation can be rewritten as accompanied by the boundary conditions ∂Q ∂y (y θ ) = − 1 + n(t) 1 + n(t − τ d ) and ∂Q ∂y (y + r ) − ∂Q ∂y (y − r ) = 1 + n(t − τ r ) 1 + n(t − τ d ) The stationary solution can be expressed as Q 0 (y) = e −y 2 F (y) y > y r Q 0 (y) = e −y 2 F (y r ) y < y r where F (y) = y θ y due u 2 .
Finally, from the definition of n it follows that n 0 = 0. As a next step, we linearize the Fokker-Planck equation around the stationary solution to treat the fluctuations in a perturbative way. Upon assuming Q = Q 0 + q and neglecting nonlinear terms in q and n,