Death and rebirth of neural activity in sparse inhibitory networks

In this paper, we clarify the mechanisms underlying a general phenomenon present in pulse-coupled heterogeneous inhibitory networks: inhibition can induce not only suppression of the neural activity, as expected, but it can also promote neural reactivation. In particular, for globally coupled systems, the number of firing neurons monotonically reduces upon increasing the strength of inhibition (neurons' death). However, the random pruning of the connections is able to reverse the action of inhibition, i.e. in a sparse network a sufficiently strong synaptic strength can surprisingly promote, rather than depress, the activity of the neurons (neurons' rebirth). Thus the number of firing neurons reveals a minimum at some intermediate synaptic strength. We show that this minimum signals a transition from a regime dominated by the neurons with higher firing activity to a phase where all neurons are effectively sub-threshold and their irregular firing is driven by current fluctuations. We explain the origin of the transition by deriving an analytic mean field formulation of the problem able to provide the fraction of active neurons as well as the first two moments of their firing statistics. The introduction of a synaptic time scale does not modify the main aspects of the reported phenomenon. However, for sufficiently slow synapses the transition becomes dramatic, the system passes from a perfectly regular evolution to an irregular bursting dynamics. In this latter regime the model provides predictions consistent with experimental findings for a specific class of neurons, namely the medium spiny neurons in the striatum.


Introduction
The presence of inhibition in excitable systems induces a rich dynamical repertoire, which is extremely relevant for biological [13], physical [32], and chemical systems [84]. In particular, inhibitory coupling has been invoked to explain cell navigation [87], morphogenesis in animal coat pattern formation [46], and the rhythmic activity of central pattern generators in many biological systems [27,45]. In brain circuits, the role of inhibition is fundamental to balance massive recurrent excitation [73] in order to generate physiologically relevant cortical rhythms [12,72].
Inhibitory networks are important not only for the emergence of rhythms in the brain, but also for the fundamental role they play in information encoding in the olfactory system [39] as well as in controlling and regulating motor and learning activity in the basal ganglia [5,15,47]. Furthermore, stimulus-dependent sequential activation of a group of neurons, reported for asymmetrically connected inhibitory cells [33,52], has been suggested as a possible mechanism to explain sequential memory storage and feature binding [67].
These explain the long-term interest for numerical and theoretical investigations of the dynamics of inhibitory networks. The study of globally coupled homogeneous systems have already revealed interesting dynamical features, ranging from full synchronization to clustering appearance [23,83,85], and from the emergence of splay states [90] to oscillator death [6]. The introduction of disorder, e.g. random dilution, noise, or other forms of heterogeneity in these systems leads to more complex dynamics, ranging from fast global oscillations [9] in neural networks and self-sustained activity in excitable systems [37], to irregular dynamics [3, 29-31, 42, 49, 56, 82, 90]. In particular, inhibitory spiking networks, due to stable chaos [63], can display extremely long erratic transients even in linearly stable regimes [3,29,30,42,49,82,89,90].
One of the most studied inhibitory neural populations is represented by medium spiny neurons (MSNs) in the striatum (which is the main input structure of the basal ganglia) [36,57]. In a series of papers, Ponzi and Wickens have shown that the main features of MSN dynamics can be reproduced by considering a randomly connected inhibitory network of conductance based neurons subject to external stochastic excitatory inputs [64][65][66]. Our study has been motivated by an interesting phenomenon reported for this model in [66]: namely, upon increasing the synaptic strength the system passes from a regularly firing regime, characterized by a large part of quiescent neurons, to a biologically relevant regime where almost all cells exhibit a bursting activity, characterized by an alternation of periods of silence and of high firing. The same phenomenology has been recently reproduced by employing a much simpler neural model [1], thus suggesting that this behavior is not related to the specific model employed, but is indeed a quite general property of inhibitory networks. However, the origin of the phenomenon and the minimal ingredients required to observe the emergence of this effect remain unclear.
In order to exemplify the problem addressed in this paper, we report in figure 1 the fraction of active neurons n A (i.e. the ones emitting at least one spike during the simulation time) as a function of the strength of the synaptic inhibition g in an heterogeneous network. For a fully coupled network, n A has a monotonic decrease with g ( figure 1(a)), while for a random sparse network, n A has a non-monotonic behavior, displaying a minimum at an intermediate strength g m ( figure 1(b)). In fully coupled networks the effect of inhibition is simply to reduce the number of active neurons (neuronal death). However, quite counter-intuitively, in the presence of dilution by increasing the synaptic strength the previously silenced neurons can return to firing (neuronal rebirth). Our aim is to clarify the physical mechanisms underlying neuronal death and rebirth, which are at the origin of the behavior reported in [1,66].
In particular, we consider a deterministic network of purely inhibitory pulse-coupled Leaky Integrate-and-Fire (LIF) neurons with an heterogeneous distribution of excitatory DC currents, accounting for the different level of excitability of the neurons. The evolution of this model is studied for fully coupled and for random sparse topology, as well as for synapses with different time courses. For the fully coupled case, it is possible to derive, within a self-consistent mean field approach, the analytic expressions for the fraction of active neurons and for the average firing frequency n as a function of coupling strength g. In this case, the monotonic decrease of n A Figure 1. Fraction of active neurons n A as a function of the inhibitory synaptic strength g for (a) a globally coupled system, where = -K N 1, and (b) a randomly connected (sparse) network with K=20. In (a),the asymptotic value n A calculated after a time =t 1 10 S 6 is reported. Conversely in (b), n A is reported at successive times: namely, t S = 985 (red squares), =t 1.1 10 S 4 (brown stars), =t 5 10 S 5 (blue diamonds), and =t 1 10 S 6 (green triangles). An estimation of the times needed to reach n A = 1 can be obtained by employing equation (13); these values range from =t 5 10 s 9 for g = 0.1, to5 10 5 for g=50. Insets in (b) depict the probability distributions n ( ) P of the single neuron firing rate ν for the sparse network for a given g at two different times: t S = 985 (red filled histograms) and =t 1 10 S 6 (thick empty green histograms). The histograms are calculated by considering only active neurons. The reported data refer to instantaneous synapses, to a system size N=400, and to an uniform distribution P(I) with with g can be interpreted as a Winner Takes All (WTA) mechanism [16,21,88], where only the most excitable neurons survive the inhibition increase. For random sparse networks, neuronal rebirth can be interpreted as a re-activation process induced by erratic fluctuations in synaptic currents. Within this framework it is possible to semi-analytically obtain, for instantaneous synapses, a closed set of equations for n A as well as for the average firing rate and coefficient of variation as a function of the coupling strength. In particular, the firing statistics of the network can be obtained via a mean field approach by extending the formulation derived in [70] to account for synaptic shot noise with constant amplitude. The introduction of a finite synaptic time scale does not modify the overall scenario provided that this is shorter than the membrane time constant. As soon as the synaptic dynamics become slower, the phenomenology of the transition is modified. At < g g m we have a frozen phase where n A does not evolve in time on the explored time scales, since the current fluctuations are negligible. Above g m we have a bursting regime, which can be related to the emergence of correlated fluctuations induced by slow synaptic times, as discussed in the framework of the adiabatic approach in [50,51].
The remainder of this paper is organized as follows: In section 2 we present the models that will be considered in the paper as well as the methods adopted to characterize their dynamics. In section 3 we consider the globally coupled network where we provide analytic self-consistent expressions accounting for the fraction of active neurons and the average firing rate. Section 4 is devoted to the study of sparsely connected networks with instantaneous synapses, and to the derivation of the set of semi-analytic self-consistent equations providing n A , the average firing rate, and the coefficient of variation. In section 5 we discuss the effect of synaptic filtering with a particular attention on slow synapses. Finally in section 6, we briefly discuss the obtained results with a focus on the biological relevance of our model.

Model and methods
We examine the dynamical properties of a heterogeneous inhibitory sparse network made of N LIF neurons. The time evolution of the membrane potential v i of the i-th neuron is ruled by the following first-order ordinary differential equation: where > g 0 is the inhibitory synaptic strength, I i is the neuronal excitability of the i-th neuron encompassing both intrinsic neuronal properties and the excitatory stimuli originating from areas outside the considered neural circuit, and E i (t) represents the synaptic current due to the recurrent interactions within the considered network. The membrane potential v i of neuron i evolves accordingly to equation (1) until it overcomes a constant threshold q = 1, which leads to the emission of a spike (action potential) transmitted to all connected post-synaptic neurons while v i is reset to its resting value v r = 0. The model in (1) is expressed in adimensional units, this amounts to assume a membrane time constant t = 1 m ; for the conversion to dimensional variables see appendix A. The heterogeneity is introduced in the model by assigning to each neuron a different value of excitability I i drawn from a flat distribution P(I), whose support is Î [ ] I l l , 1 2 with  q l 1 ; therefore, all neurons are supra-threshold.
The synaptic current E i (t) is given by the linear superposition of all the inhibitory post-synaptic potentials (IPSPs) h ( ) t emitted at previous times < t t n j by the pre-synaptic neurons connected to neuron i, namely where K is the number of pre-synaptic neurons. C ij represents the elements of the N×N connectivity matrix associated with an undirected random network, whose entries are 1 if there is a synaptic connection from neuron j to neuron i, and 0 otherwise. For the sparse network, we randomly select the matrix entries; however, to reduce the sources of variability in the network, we assume that the number of pre-synaptic neurons is fixed, namely å = for each neuron i, where autaptic connections are not allowed. We have verified that the results do not change if we randomly choose the links accordingly to an Erdös-Renyi distribution with a probability K/N. For a fully coupled network we have = -K N 1. The shape of the IPSP characterizes the type of filtering performed by the synapses on the received action potentials. We have considered two kinds of synapses: instantaneous ones, where h d = ( ) ( ) t t , and synapses where the PSP is an α-pulse, namely with H denoting the Heaviside step function. In this latter case the rise and decay time of the pulse are the same, namely t a = a 1 , and therefore the pulse duration t P can be assumed to be twice the characteristic time t a . Model equations (1) and (2) are integrated exactly in terms of the associated event driven maps for different synaptic filtering, which correspond to Poincaré maps performed at the firing times (for details see appendix A) [53,91].
For instantaneous synapses, we have usually considered system sizes N=400 and N = 1400, and in the sparse case in-degrees   K 20 80 and   K 20 600, respectively, with integration times up to =t 1 10 S 6 . For synapses with a finite decay time we have limited the analysis to N=400 and K=20, and to maximal integration times =t 1 10 S 5 . Finite size dependencies on N are negligible with these parameter choices, as we have verified.
In order to characterize the network dynamics, we measure the fraction of active neurons ( ) n t A S at time t S , i.e. the fraction of neurons emitting at least one spike in the time interval [ ] t 0, S . Therefore a neuron will be considered silent if it has a frequency smaller than t 1 S , and with our choices of = t 10 10 Furthermore, for each neuron we estimate the time averaged inter-spike interval (ISI) T ISI , the associated firing frequency n = T 1 ISI , as well as the coefficient of variation CV, which is the ratio of the standard deviation of the ISI distribution divided by T ISI . For a regular spike train CV=0, and for a Poissonian distributed one CV = 1, while > CV 1 is an indication of bursting activity. The indicators reported in the following to characterize network activity are ensemble averages over all active neurons, which we denote asā for a generic observable a.
To analyze the linear stability of the dynamical evolution we measure the maximal Lyapunov exponent λ, which is positive for chaotic evolution, and negative (zero) for stable (marginally stable) dynamics [4]. In particular, by following [2,55], λ is estimated by linearizing the corresponding event driven map.

Fully coupled networks: WTA
In the fully coupled case we observe that the fraction of active neurons n A saturates, after a short transient, to a value that remains constant in time. In this case, it is possible to derive a self-consistent mean field approach to obtain analytic expressions for the fraction of active neurons n A and for the average firing frequency n of neurons in the network. In a fully coupled network each neuron receives the spikes emitted by the other = -K N 1 neurons; therefore, each neuron is essentially subject to the same effective input> μ, apart from the finite size corrections ( ) N 1 . The effective input current, for a neuron with an excitability I, is given by is the number of active pre-synaptic neurons assumed to fire with the same average frequency n .
In a mean field approach, each neuron can be seen as isolated from the network and driven by the effective input current μ. Taking into account the distribution of excitabilities P(I), one obtains the following selfconsistent expression for the average firing frequency: is the fraction of active neurons. In (5) we have used the fact that for an isolated LIF neuron with constant excitability C, the ISI is simply given by [11]. An implicit expression for n A can be obtained by estimating the neurons with effective input m q > ; in particular, the fraction of silent neurons is given by where l 1 is the lower limit of the support of the distribution, while * n q = + l g n A . By solving self-consistently equations (5) and (6), one can obtain the analytic expression for n A and n for any distribution P(I).
In particular, for excitabilities distributed uniformly in the interval [ ] l l , 1 2 , the expression for the average frequency equation with the constraint that n A cannot be larger than one. The analytic results for these quantities compare quite well with the numerical findings estimated for different distribution intervals [ ] l l , 1 2 , coupling strengths, and system sizes, as shown in figure 2. For definitively large coupling g>10, some discrepancies between the mean field estimations and the simulation results are observable (see figure 2(b)). These differences are probably due to the discreteness of the pulses, which cannot be neglected for very large synaptic strengths.
As a general feature we observe that n A is steadily decreasing with g, thus indicating that a group of neurons with higher effective inputs (winners) silence the other neurons (losers) and that the number of winners eventually vanishes for sufficiently large coupling in the limit of large system sizes. Furthermore, the average excitability of active neurons (the winners)Ī A increases with g, as shown in the inset of figure 2(a), thus revealing that only neurons with higher excitabilities survive the silencing action exerted by the other neurons. At the same time, as an effect of the growing inhibition, the average firing rate of the winners dramatically slows down. Therefore, despite the increase ofĪ A , the average effective input m indeed decreases for increasing inhibition. This represents a clear example of the WTA mechanism obtained via (lateral) inhibition, which has been shown to have biological relevance for neural systems [20,60,88].
It is important to understand what is the minimal coupling value g c for which the firing neurons start to die. In order to estimate g c it is sufficient to set n A = 1 in equations (7) and (8). In particular, one gets q n = -( )¯( ) g l , 9 c 1 and thus, for even an infinitesimally small coupling is in principle sufficient to silence some neurons. Furthermore, from figure 3(a) it is evident that whenever the excitabilities become homogeneous, i.e. for  l l 1 2 , the critical synaptic coupling g c diverges toward infinity. Thus, heterogeneity in the excitability distribution is a necessary condition in order to observe a gradual neuronal death, as shown in figure 2(a). This is in agreement with the results reported in [7], where homogeneous fully coupled networks of inhibitory LIF neurons have been examined. In particular, for finite systems and slow synapses, the authors in [7] revealed the existence of a sub-critical Hopf bifurcation from a fully synchronized state to a regime characterized by oscillator death occurring at some critical g c . However, in the thermodynamic limit  ¥ g c for fast as well as slow synapses, which is in agreement with our mean field result for instantaneous synapses.
We also proceed to investigate the isolines corresponding to the same critical g c in the ( ) l l , 1 2 -plane, and the results are reported in figure 3(b) for three selected values of g c . It is evident that the l 1 and l 2 -values associated with the isolines display a direct proportionality among them. However, despite lying on the same g c -isoline, different parameter values induce a completely different behavior of n A as a function of the synaptic strength, as shown in the inset of figure 3 Direct simulations of the network at finite sizes, namely for N=400 and N=1400, show that sufficiently large coupling neurons with similar excitabilities tend to form clusters, similarly to what was reported in [42], but with a delayed pulse transmission. However, in contrast to [42], the overall macroscopic dynamics is asynchronous, and no collective oscillations can be detected for the whole range of considered synaptic strengths.

Sparse networks: neuronal rebirth
In this section we will consider a network with sparse connectivity, namely each neuron is supra-threshold and it receives instantaneous IPSPs from  K N randomly chosen neurons in the network. Due to the sparseness, the input spike trains can be considered as uncorrelated, and at a first approximation it can be assumed that each spike train is Poissonian with a frequency n corresponding to the average firing rate of the neurons in the network [8,9]. Usually, the mean activity of a LIF neural network is estimated in the context of the diffusion approximation [69,80]. This approximation is valid whenever the arrival frequency of the IPSPs is high with respect to the firing emission, while the amplitude of each IPSPs (namely, = G g K) is small with respect to the firing threshold θ. This latter hypothesis in our case is not valid for sufficiently large (small) synaptic strength g (in-degree K ), as can be appreciated by the comparison shown in figure 13 in appendix B. Therefore, the synaptic inputs should be treated as shot noise. In particular, here we apply an extended version of the analytic approach derived by Richardson and Swabrick in [70] to estimate the average firing rate and the average coefficient of variation for LIF neurons with instantaneous synapses subject to inhibitory shot noise of constant amplitude (for more details see appendices B and C).
Contrasting with the fully coupled case, the fraction of active neurons n A does not saturate to a constant value for sufficiently short times. Instead, n A increases in time due to the rebirth of losers previously silenced by the firing activity of winners, as shown in figure 1(b). This effect is clearly illustrated by considering the probability distributions n ( ) P of the firing rates of the neurons at successive integration times t S . These are reported in the insets of figure 1(b) for two coupling strengths and two times: namely, t S = 985 (red lines) and =t 1 10 S 6 (green lines). From these data is evident that the fraction of neurons with low firing rate (the losers) increases with time, while the fraction of high firing neurons remains almost unchanged. Moreover, the variation of n A slows down for increasing t S , and n A approaches some apparently asymptotic profile for sufficiently long integration times. Furthermore, n A has a non-monotonic behavior with g, which is opposite to the fully coupled case. In particular, n A reveals a minimum n A m at some intermediate synaptic strength g m followed by an increase toward n A = 1 at large g. As we have verified, as long as <  K N 1 , finite size effects are negligible and the actual value of n A depends only on the in-degree K and the considered simulation time t S . In the following we will try to explain the origin of such a behavior.
Despite the model being fully deterministic, due to the random connectivity the rebirth of silent neurons can be interpreted in the framework of activation processes induced by random fluctuations. In particular, we can assume that each neuron in the network will receive n K A independent Poissonian trains of inhibitory kicks of constant amplitude G characterized by an average frequency n ; thus, each synaptic input can be regarded as a single Poissonian train with total frequency n =R n K A . Therefore, each neuron, characterized by its own excitability I, will be subject to an average effective input m ( ) I (as reported in equation (4)) plus fluctuations in the synaptic current of intensity 1 0 A Indeed, we have verified that (10) gives a quantitatively correct estimation of the synaptic current fluctuations over the whole range of synaptic coupling considered (as shown in figure 4). A closer analysis of the probability distributions P(IAT) of the inter-arrival times (IATs) reveals that these are essentially exponentially distributed, as expected for Poissonian processes, with a decay rate given by R, as evident from figure 5 for two different synaptic strengths. However, all these indications are not sufficient to guarantee that the IAT statistics are indeed Poissonian. In particular, as pointed out in [40], a superposition of uncorrelated spike trains generated by the same non-Poissonian renewal process can result in a peculiar non-renewal process characterized by exponentially distributed and uncorrelated IATs with a non-flat power spectrum. In our case we have indeed verified that for small coupling the power spectrum associated with the IATs deviates from the flat one at low frequencies, which is similar to the results reported in [40]; meanwhile, at large g the spectrum recovers a Poissonian shape. Therefore the hypothesis that the neuronal input is Poissonian in our case should be considered only as a first-order approximation, in particular for small synaptic couplings. For instantaneous IPSP, the current fluctuations are due to stable chaos [63] since the maximal Lyapunov exponent is negative for the whole range of coupling, as we have verified. Therefore, as reported by many authors, erratic fluctuations in inhibitory neural networks with instantaneous synapses are due to finite amplitude instabilities, while at the infinitesimal level the system is stable [3,29,30,42,49,82,90].
In this circumstance, the silent neurons stay in a quiescent state corresponding to the minimum of the , and in order to fire they should overcome a barrier  q m The average time t A required to overcome such a barrier can be estimated accordingly to the Kramers' theory for activation processes [25,80], namely where t 0 is an effective time scale taking in account the intrinsic non-stationarity of the process, i.e. the fact that the number of active neurons increases during the time evolution.
It is important to stress that the expression (11) will remain valid also in the limit of large synaptic couplings, not only because of s 2 , but also because the barrier height will increase with g. Furthermore, both these quantities grow quadratically with g at sufficiently large synaptic strength, as it can be inferred from equations (4) and (10).
It is reasonable to assume that at a given time t S all neurons with < t t A S will have fired at least once, and that the more excitable will fire first. Therefore, by assuming that the fraction of active neurons at time t S is ( ) n t A S , the last neuron that has fired should be characterized by the following excitability: Here, excitabilities I are uniformly distributed in the interval [ ] l l , 1 2 . In order to obtain an explicit expression for the fraction of active neurons at time t S , one should solve the equation (11) for the neuron with excitabilityÎ by setting t S = t A , thus obtaining the following solution   (4) (equation (10)) and averaged only over the active neurons. The data refer to N=1400, K=140, , and to a simulation time =t 1 10 S 6 . In both panels, the red continuous line indicates the exponential distribution corresponding to a purely Poissonian process with an arrival rate given by n =R n K A , and the dashed blue vertical lines refer to the average IAT for the Poissonian distribution, namely R 1 . The distributions have been evaluated for the arrival of5 10 5 IPSPs. Other parameters used for the simulation are as in figure 1 Equation (13) gives the dependence of n A on the coupling strength g for a fixed integration time t S and time scale t 0 whenever we can provide the value of the average frequency n . A quick inspection of equations (11) and (13) shows that, setting n A = 1, we obtain two solutions for the critical couplings g c1 (g c2 ) below (above) that all neurons will fire at least once in the considered time interval. The solutions are reported in figure 6. In particular, we observe that whenever q  l 1 the critical coupling g c1 will vanish, which is analogous to the fully coupled situation. These results clearly indicate that n A should display a minimum for some finite coupling strength Î ( ) g g g , m c 1 c 2 . Furthermore, as shown in figure 6 the two critical couplings approach one another for increasing t S and finally merge, indicating that at sufficiently long times all neurons will be active at any synaptic coupling strength g.
The average frequency n can be obtained analytically by following the approach described in appendix B for LIF neurons with instantaneous synapses subject to inhibitory Poissonian spike trains. In particular, the selfconsistent expression for the average frequency reads as  (13) and (14) provides a theoretical estimation of n A and n for the whole considered range of synaptic strength, once the unknown time scale t 0 is fixed. This time scale has been determined via an optimal fitting procedure for sparse networks with N=400 and K=20, 40, and 80 at a fixed integration time =t 1 10 S 6 . The results for n A are reported in figure 7(a). The estimated curves reproduce reasonably well the numerical data for K=20 and 40, while for K=80 the agreement worsens at large coupling strengths. This could be due to the fact that by increasing g and K, the spike trains stimulating each neuron cannot be assumed to be completely independent, as done for the derivation of equations (13) and (14). Nevertheless, the average frequency n is quantitatively well reproduced for the considered K values over the entire range of synaptic strengths, as is evident from figures 7(b)-(d). A more detailed comparison between the theoretical estimations and the numerical data can be obtained by considering the distributions n ( ) P of the single neuron firing rate for different coupling strengths reported in figures 7(k)-(m) for K=40. The overall agreement can be considered as more than satisfactory, while observable discrepancies are probably due to the fact that our approach neglected a further source of disorder present in the system and related to heterogeneity in the number of active pre-synaptic neurons [9].
We have also analytically estimated the average coefficient of variation of the firing neurons CV by extending the method derived in [70] to obtain the response of a neuron receiving synaptic shot noise input. The analytic expressions of the coefficient of variation for LIF neurons subject to inhibitory shot noise with fixed postsynaptic amplitude are obtained by estimating the second moment of the associated first-passage-time distribution; the details are reported in appendix C. The coefficient of variation can be estimated once the self- consistent values for n A and n have been obtained via equations (13) and (14). The comparison with the numerical data, reported in figures 7(e)-(g), reveals a good agreement over the whole range of synaptic strengths for all considered in-degrees.
At sufficiently small synaptic coupling the neurons fire tonically and almost independently, as shown by the raster plot in figure 7(h) and by the fact that n approaches the average value for the uncoupled system (namely, 0.605) and  CV 0. Furthermore, the neuronal firing rates are distributed toward finite values, indicating that the inhibition has a minor role in this case, as shown in figure 7(k). By increasing the coupling, n A decreases, and as an effect of the inhibition more and more neurons are silenced (as evident from figure 7(l)) and the average firing rate decreases; at the same time, the dynamics become slightly more irregular, as shown in figure 7(i). At large coupling > g g m , a new regime appears, where almost all neurons become active but with extremely slow dynamics that are essentially stochastic with  CV 1, as also testified by the raster plot reported in figure 7(j) and by the firing rate distribution shown in figure 7(m).
Furthermore, from figure 7(a) it is clear that the minimum value of the fraction of active neurons n Am decreases by increasing the network in-degree K, while g m increases with K. This behavior is further investigated in a larger network, namely N=1400, and reported in the inset of figure 8. It is evident that n A stays close to the globally coupled solutions over larger and larger intervals for increasing K. This can be qualitatively understood by the fact that the current fluctuations in equation (10), responsible for the rebirth of silent neurons, are proportional to g and scales as K 1 . Therefore, at larger in-degrees the fluctuations have similar intensities only for larger synaptic coupling. of the single neuron firing rate ν, for the same values of g in the panels above. Empty-discontinuous bars correspond to the theoretical prediction while filled bars indicate the histogram calculated with the simulation. The remaining parameters are as in figure 1 The general mechanism behind neuronal rebirth can be understood by considering the values of the effective neuronal input and current fluctuations as a function of g. As shown in figure 4, the effective input current m A , averaged over the active neurons, essentially coincides with the average excitabilityĪ A for  g 0, where the neurons can be considered as independent from others. The inhibition leads to a decrease of m A , and to a crossing of the threshold θ at exactly = g g m . This indicates that at < g g m the active neurons, being on average supra-threshold, fire almost tonically inhibiting the losers via a WTA mechanism. In this case the firing neurons are essentially mean-driven and the current fluctuations play a role in the rebirth of silent neurons only on extremely long time scales; this is confirmed by the low values of s in such a range, as evident from figure 4. On the other hand, for > g g m , the active neurons are now on average below threshold while fluctuations dominate dynamics. In particular, the firing is now extremely irregular mainly due to re-activation processes. Therefore, the origin of the minimum in n A can be understood as a transition from a mean-driven to a fluctuation-driven regime [68].
A quantitative definition of g m can be given by requiring that the average input current of the active neurons m A crosses the threshold θ at = g g m , namely whereĪ A is the average excitability of the firing neurons, while n Am and n m are the fraction of active neurons and the average frequency at the minimum. For a uniform distribution P(I), a simpler expression for g m can be derived, namely We have compared the numerical measurements of g m with the estimations obtained by employing equation (15), where n A m and n are obtained from the simulations. As shown in figure 8, for a network of size N = 1400, the overall agreement is more than satisfactory for in-degrees ranging over almost two decades (namely, for   K 20 600). This confirms that our guess (that the minimum n A m occurs exactly at the transition from mean-driven to fluctuation-driven dynamics) is consistent with the numerical data for a wide range of in-degrees.
It should be stressed that, as we have verified for various system sizes (namely, N = 700, 1400, and 2800) and for a constant average in-degree K=140, for instantaneous synapses the network is in an heterogeneous asynchronous state for all considered values of synaptic coupling. This is demonstrated by the fact that the intensity of the fluctuations of average firing activity, measured by considering the low-pass filtered linear superposition of all the spikes emitted in the network, vanishes as N 1 [85]. Therefore, the observed transition at = g g m is not associated with the emergence of irregular collective behaviors as reported for globally coupled heterogeneous inhibitory networks of LIF neurons with delay [42] or pulse-coupled phase oscillators [82].

Effect of synaptic filtering
In this section we will investigate how synaptic filtering can influence the previously reported results. In particular, we will consider non-instantaneous IPSP with an α-function profile (3), whose evolution is ruled by a single time scale t a .

Fully coupled networks
Let us first examine the fully coupled topology. In this case we analogously observe the δ-pulse coupling by increasing the inhibition so that the number of active neurons steadily decreases toward a limit where only few neurons (or eventually only one) will survive. At the same time the average frequency also decreases monotonically, as shown in figure 9 for two different t a differing by almost two orders of magnitude. Furthermore, the mean field estimations (7) and (8) obtained for n A and n represent a very good approximation for α-pulses (as shown in figure 9). In particular, the mean field estimation essentially coincides with the numerical values for slow synapses, as evident from the data reported in figure 9 for t = a 10 (black filled circles). This can be explained by the fact that for sufficiently slow synapses, with t >T P I S I , the neurons feel the synaptic input current as continuous because each input pulse has essentially no time to decay between two firing events. Therefore, the mean field approximation for the input current(4) works well in this case. This is particularly true for t = a 10, where t = 20 P and - T 2 6 ISI in the range of the considered coupling. While for t = a 0.125, we observe some deviation from the mean field results (red squares in figure 9). The reason for these discrepancies resides in the fact that t <T P I S I for any coupling strength, and therefore the discreteness of the pulses cannot be completely neglected, particularly for large amplitudes (large synaptic couplings). This is analogous to what is observed for instantaneous synapses.

Sparse networks
For sparse networks, n A has the same qualitative behavior as a function of the synaptic inhibition observed for instantaneous IPSPs, as shown in figure 9(b) and figure 10(a). The value of n A decreases with g and reaches a minimal value at g m ; afterwards, it increases towards n A = 1 at larger coupling. The origin of the minimum of n A as a function of g is the same as for instantaneous synapses. For < g g m , active neurons are subject to, on average, a supra-threshold effective input m A while at larger coupling m q < A , as shown in the inset of figure 10(b). This is true for any value of t a ; however, this transition from mean-to fluctuation-driven dynamics becomes dramatic for slow synapses. As evidenced from the data for the average output firing rate n and the average coefficient of variation CV , these quantities have almost discontinuous jumps at = g g m , as shown in figure 11. Therefore, let us first concentrate on slow synapses with t a larger than the membrane time constant, which is one for adimensional units. For < g g m the fraction of active neurons is frozen in time, at least on the considered time scales, as revealed by the data in figure 9(b). Furthermore, for < g g m , the mean field approximation obtained for the fully coupled case works almost perfectly both for n A and n , as reported in figure 10(a). The frozen phase is characterized by extremely small values of current fluctuations s (as shown figure 10(b)) and a quite high firing rate n - 0.4 0.5 with an associated average coefficient of variation CV of almost zero (see black circles and red squares in figure 11). Instead, for > g g m the number of active neurons increases in time similarly to what is observed for instantaneous synapses, while the average frequency becomes extremely small n - 0.04 0.09 and the value of the coefficient of variation becomes definitely larger than one. These effects can be explained by the fact that below g m the active neurons (the winners) are subject to an effective input m q > A that induces a quite regular firing, as testified by the raster plot displayed in figure 10(c). The supra-threshold activity of the winners joined together with the filtering action of the synapses guarantee that on average each neuron in the network receives an almost continuous current with small fluctuations in time. These results explain why the mean field approximation still works in the frozen phase, where fluctuations in synaptic currents are essentially negligible and unable to induce any neuronal rebirth, at least on realistic time  scales. In this regime the only mechanism in action is the WTA, and fluctuations begin to have a role for slow synapses only for > g g m . Indeed, as shown in figure 10(b), the synaptic fluctuations s for t = a 10 (black circles) are almost negligible for < g g m and show an enormous increase of almost two orders of magnitude at = g g m . Similarly, at t = a 2 (red square) a noticeable increase of s is observable at the transition. In order to better understand the abrupt changes in n and CV observable for slow synapses at = g g m , let us consider the case t = a 10. As shown in figure 11(c), t > - T 2 3 P I S I for < g g m . Therefore, for these couplings the IPSPs have no time to decay between a firing emission and the next one, and thus the synaptic fluctuations s are definitely small in this case. At g m an abrupt jump is observable to large values where t > T ISI P , which is due to the fact that now the neurons display bursting activities, as evident from the raster plot shown in figure 10(d). The bursting is due to the fact that, for > g g m , the active neurons are subject to an effective input, which is on average sub-threshold; therefore, the neurons tend to be silent. However, due to current fluctuations, a neuron can pass the threshold and silent periods can be interrupted by bursting phases where the neuron fires almost regularly. In fact, the silent (inter-burst) periods are very long -700 900 compared to the duration of the bursting periods, namely -25 50, as shown in figure 11(c). This explains the abrupt decrease of the average firing rate reported in figure 11(a). Furthermore, the inter-burst periods are exponentially distributed with an associated coefficient of variation -0.8 1.0, which clearly indicates the stochastic nature of the switching from the silent phase to the bursting phase. The firing periods within the bursting phase are instead quite regular, with an associated coefficient of variation 0.2, and with a duration similar toT ISI measured in the frozen phase (shaded gray circles in figure 11(c)). Therefore, above g m the distribution of the ISI exhibits a long exponential tail associated with the bursting activity, and this explains the very large values of the measured coefficient of variation. By increasing coupling, the fluctuations in the input current become larger, and thus the fraction of neurons that fires at least once within a certain time interval increases. At the same time, n , the average inter-burst periods, and the firing periods within the bursting phase remain almost constant at > g 10, as shown in figure 11(a). This indicates that the decrease of m A and increase of s due to the increased inhibitory coupling essentially compensate for each other in this range. Indeed, we have verified that for t = a 10 and t = a 2, m A (s ) decreases (increases) linearly with g with similar slopes, namely m - For faster synapses, the frozen phase is no longer present. Furthermore, due to rebirths induced by current fluctuations, n A is always larger than the fully coupled mean field result (8), even at < g g m . It is interesting to notice that by decreasing t a , we are now approaching the instantaneous limit, as indicated by the results reported for n A in figure 10(a) and CV in figure 11(b). In particular, for t = a 0.125 (green triangles) the data almost collapses on the corresponding values measured for instantaneous synapses in a sparse network with the same characteristics and over a similar time interval (dashed line). Furthermore, for fast synapses with t < a 1 the bursting activity is no longer present, as can be appreciated by the fact that at most CV approaches one in the very large coupling limit.
For sufficiently slow synapses, the average firing rate n can be estimated by applying the so-called adiabatic approach developed by Moreno-Bote and Parga in [50,51]. This method applies to LIF neurons with a synaptic time scale longer than the membrane time constant. In these conditions, the output firing rate can be reproduced by assuming that the neuron is subject to an input current with time-correlated fluctuations, which can be represented as colored noise with a correlation time given by the pulse duration t t = a 2 P (for more details see appendix D). In this case we are unable to develop a self-consistent approach to obtain at the same time n A and the average frequency. However, once n A is provided by simulations, the estimated solution to (45) obtained with the adiabatic approach gives very good agreement with the numerical data for sufficiently slow synapses, namely for  t 1 P , as shown in figure 11(a) for t = a 10, 2, and 0.5. The theoretical expression (45) is even able to reproduce the jump in average frequencies observable at g m , and can therefore capture the bursting phenomenon. By considering t < 1 P , as expected, the theoretical expression fails to reproduce the numerical data, particularly at large coupling (see the dashed green line in figure 11(a) corresponding to t = a 0.125). By following the arguments reported in [50], the bursting phenomenon observed for t > a 1 and > g g m can be interpreted at a mean field level as the response of a sub-threshold LIF neuron subject to colored noise with correlation t P . In this case, the neuron is definitely sub-threshold, but in the presence of a large fluctuation it can lead to firing, and due to the finite correlation time, it can remain supra-threshold regularly firing for a period t  P . The validity of this interpretation is confirmed by the fact that the measured average bursting periods are of the order of the correlation time t t = a 2 P , namely, -27 50 ( -7 14) for t = a 10 (t = a 2). As a final point, to better understand the dynamical origin of the measured fluctuations in this deterministic model, we have estimated the maximal Lyapunov exponent λ. As expected from previous analysis, for noninstantaneous synapses we can observe the emergence of regular chaos in purely inhibitory networks [1,30,89]. In particular, for sufficiently fast synapses, we typically note a transition from a chaotic state at low coupling to a linearly stable regime (with l < 0) at large synaptic strengths, as shown in figure 12(a) for t = a 0.125. This is despite the fact that current fluctuations are monotonically increasing with synaptic strength. Therefore, fluctuations are due to chaos at small coupling, while at larger g they are due to finite amplitude instabilities, as expected for stable chaotic systems [3]. However, the passage from positive to negative values of the maximal Lyapunov exponent is not related to the transition occurring at g m from mean-driven to a fluctuation-driven dynamics in the network.
For slow synapses, λ is essentially zero at small coupling in the frozen phase, characterized by tonic spiking of the neurons, but becomes positive by approaching g m . For larger synaptic strengths λ, after reaching a maximal value, it decreases and eventually becomes negative at  g g m , as reported in figures 12(b)-(c). Only for extremely slow synapses, as shown in figure 12(c) for t = a 10, the chaos onset seems to coincide with the transition occurring at g m . These findings are consistent with recent results concerning the emergence of asynchronous rate chaos in homogeneous inhibitory LIF networks with deterministic [26] and stochastic [31] evolution. However, a detailed analysis of this aspect goes beyond the scope of the present paper.

Discussion
In this paper we have shown that the effect reported in [1,66] is observable whenever two sources of quenched disorder are present in the network: namely, a random distribution of neural properties and a random topology. In particular, we have shown that neuronal death due to synaptic inhibition is observable only for heterogeneous distributions of neural excitabilities. Furthermore, in a globally coupled network the less excitable neurons are silenced for increasing synaptic strength until only one or few neurons remain active. This scenario corresponds to the WTA mechanism via lateral inhibition, which has often been invoked in neuroscience to explain several brain functions [88]. WTA mechanisms have been proposed to model hippocampal CA1 activity [16], as well as to be at the basis of visual velocity estimates [24], and to be essential for controlling visual attention [28].
However, most brain circuits are characterized by sparse connectivity [10,38,60], and in these networks we have shown that an increase in inhibition can lead from a phase dominated by neuronal death to a regime where neuronal rebirths occur. Therefore, the growth of inhibition can have the counter-intuitive effect to activate silent neurons due to the enhancement of current fluctuations. The reported transition is characterized by a passage from a regime dominated by the almost tonic activity of a group of neurons, to a phase where subthreshold fluctuations are at the origin of the irregular firing of a high number of neurons in the network. For instantaneous synapses, the first and second moment of the firing distributions have been obtained together with the fraction of active neurons using a mean field approach, where the neuronal rebirth is interpreted as an activation process driven by synaptic shot noise [70].
For a finite synaptic time smaller than the characteristic membrane time constant, the scenario is similar to that observed for instantaneous synapses. However, the transition from mean-driven to fluctuation-driven dynamics becomes dramatic for sufficiently slow synapses. In this situation one observes for low synaptic strength a frozen phase, where synaptic filtering washes out the current fluctuations, thus leading to extremely regular dynamics controlled only by a WTA mechanism. As soon as the inhibition is sufficiently strong to lead the active neurons below threshold, neuronal activity becomes extremely irregular, exhibiting long silent phases interrupted by bursting events. The origin of these bursting periods can be understood in terms of the emergence of correlations in current fluctuations induced by the slow synaptic time scale, as explained in [50].
In our model, the random dilution of network connectivity is a fundamental ingredient to generate current fluctuations, whose intensity is controlled by the average network in-degree K. A natural question is whether the reported scenario will remain observable in the thermodynamic limit. On the basis of previous studies we can affirm that this depends on how K scales with the system size [22,41,75]. In particular, if K stays finite for  ¥ N the transition will remain observable. For K diverging with N, the fluctuations become negligible for All the reported data were calculated for a system size N=400 and K=20 and for simulation times  t 5 10 7 10 4 S 5 , thus ensuring a good convergence of λ to its asymptotic value. The other parameters are as in figure 9.
sufficiently large system sizes, impeding neuronal rebirths, and the dynamics will be controlled only by the WTA mechanism.
An additional source of randomness present in the network is related to the variability in the number of active pre-synaptic neurons. In our mean field approach we have assumed that each neuron is subject to n K A spike trains; however, this is true only on average. The number of active pre-synaptic neurons is a random variable binomially distributed with average n K A and variance - . Future developments of the theoretical approach reported here should include also such variability in modeling network dynamics [9].
Furthermore, we show that the considered model is not chaotic for instantaneous synapses; in such a case, we observe irregular asynchronous states due to stable chaos [63]. The system can become truly chaotic for only finite synaptic times [3,30]. However, we report clear indications that for synapses faster than the membrane time constant t m the passage from mean-driven to fluctuation-driven dynamics is not related to the onset of chaos. Only for extremely slow synapses do we have numerical evidence that the appearance of the bursting regime could be related to a passage from a zero Lyapunov exponent to a positive one. This is in agreement with the results reported in [26,31] for homogeneous inhibitory networks. These preliminary indications demand more detailed investigations of deterministic spiking networks in order to relate fluctuation-driven regimes and chaos onsets. Moreover, we expect that it will be hard to distinguish whether the erratic current fluctuations are due to regular chaos or stable chaos on the basis of network activity analysis, as also pointed out in [30].
Concerning the biological relevance of the presented model, we can attempt a comparison with experimental data obtained for MSNs in the striatum. This population of neurons is fully inhibitory with sparse lateral connections (connection probability ;10%-20% [77,81]) that are unidirectional and relatively weak [78]. Furthermore, for MSNs within the same collateral network the axonal propagation delays are quite small ;1-2 ms [76] and can be safely neglected. The dynamics of these neurons in behaving mice reveals a low average firing rate with irregular firing activity (bursting) with an associated large coefficient of variation [47]. As we have shown, these features can be reproduced by sparse networks of LIF neurons with sufficiently slow synapses at > g g m and t t > a m . For values of the membrane time constant that are comparable to those measured for MSNs [59,61] (namely, t  -10 20 m msec), the model is able to capture some of the main aspects of MSNs dynamics, as shown in table 1. We obtain a reasonable agreement with the experiments for sufficiently slow synapses, where the interaction among MSNs is mainly mediated by GABA A receptors, which are characterized by IPSP durations of the order of ;5-20 ms [34,81]. However, apart the burst duration, which is definitely shorter, all other aspects of the MSN dynamics can be already captured for t t = a 2 m (with t = 10 m ms), as shown in table 1. Therefore, we can safely affirm, as also suggested in [66], that the fluctuation-driven regime emerging at > g g m is the most appropriate in order to reproduce the dynamical evolution of this population of neurons.
Other inhibitory populations are present in the basal ganglia. In particular, two coexisting inhibitory populations, arkypallidal (Arkys) and prototypical (Protos) neurons, have been recently discovered in the external globus pallidus [43]. These populations have distinct physiological and dynamical characteristics, and have been shown to be fundamental for action suppression during the performance of behavioral tasks in rodents [44]. Protos are characterized by a high firing rate 47 Hz and a not too large coefficient of variation (namely,  CV 0.58) both in awake and slow wave sleep (SWS) states; meanwhile, Arkys have clear bursting dynamics with  CV 1.9 [18,44]. Furthermore, the firing rate of Arkys is definitely larger in the awake state (namely, 9 Hz) with respect to the SWS state, where firing rates are  -3 5 Hz [44]. On the basis of our results, on the one hand, Protos can be modeled as LIF neurons with reasonably fast synapses in a mean-driven regime, namely with synaptic coupling < g g m . On the other hand, Arkys should be characterized by IPSP with definitely longer durations, and should be in a fluctuation-driven phase as suggested from the results reported in figure 11. Since, as shown in figure 11(a), the firing rate of inhibitory neurons decreases by increasing synaptic strength g, we expect that the passage from awake to SWS should be characterized by a reinforcement of Arkys synapses. Our conjectures about Arkys and Protos synaptic properties based on their dynamical behaviors ask for experimental verification, which we hope will happen shortly. Besides the straightforward applicability of our findings to networks of pulse-coupled oscillators [48], it has been recently shown that LIF networks with instantaneous and non-instantaneous synapses can be transformed into the Kuramoto-Daido model [17,35,62]. Therefore, we expect that our findings should extend to phase oscillator arrays with repulsive coupling [79]. This will allow for a wider applicability of our results, due to the relevance of limit-cycle oscillators not only for modeling biological systems [86], but also for many scientific and technological applications [19,58,71,74].

Acknowledgments
Some preliminary analysis on the instantaneous synapses has been performed in collaboration with A Imparato; the complete results will be reported elsewhere [54]. We thank J Berke, B Lindner, G Mato, G Giacomelli, S Gupta, A Politi, MJE Richardson, R Schmidt, and M Timme for useful discussions. This work has been partially supported by the European Commission under the program 'Marie Curie Network for Initial Training', through project N.

Appendix A. Event driven maps
By following [53,91] the ordinary differential equations (1) and (2) describing the evolution of the membrane potential of neurons can be rewritten exactly as discrete time maps connecting successive firing events occurring in the network. In the following we will explicitly report such event driven maps for the case of instantaneous and α-synapses.
For instantaneous PSPs, the event driven map for neuron i takes the following expression: is the ISI associated with two successive neuronal firings. This latter quantity is calculated from the following expression: For α-pulses, the evolution of the synaptic current E i (t), stimulating the i-th neuron can be expressed in terms of a second-order differential equation, namely å å The model so far introduced contains only adimensional units; however, the evolution equation for the membrane potential (1) can be easily re-expressed in terms of dimensional variables as follows: For the other parameters/variables the transformation to physical units is simply given by = + - th r where = -V 60 r mV and = -V 50 th mV are realistic values of the membrane reset and threshold potential. The isolated i-th LIF neuron is supra-threshold whenever > I V r where r ( ) t is the instantaneous firing rate of the neuron. The flux can be decomposed in an average drift term plus the inhibitory part, namely = - The set of equations (28) to (30) is complemented by the boundary conditions q r , 0 , 0 , inh and by the requirement that the membrane potential distribution should be normalized, i.e ò = q -¥ ( ) P v t dv , 1 .

By introducing bilateral Laplace transforms
and by performing some algebra along the lines described in [70] it is possible to derive the analytic expression for the average firing rate ò n = -