A robust balancing mechanism for spiking neural networks

Dynamical balance of excitation and inhibition is usually invoked to explain the irregular low firing activity observed in the cortex. We propose a robust nonlinear balancing mechanism for a random network of spiking neurons, which works also in absence of strong external currents. Biologically, the mechanism exploits the plasticity of excitatory-excitatory synapses induced by short-term depression. Mathematically, the nonlinear response of the synaptic activity is the key ingredient responsible for the emergence of a stable balanced regime. Our claim is supported by a simple self-consistent analysis accompanied by extensive simulations performed for increasing network sizes. The observed regime is essentially fluctuation driven and characterized by highly irregular spiking dynamics of all neurons.

Dynamical balance of excitation and inhibition is usually invoked to explain the irregular low firing activity observed in the cortex.We propose a robust nonlinear balancing mechanism for a random network of spiking neurons, which works also in absence of strong external currents.Biologically, the mechanism exploits the plasticity of excitatory-excitatory synapses induced by short-term depression.Mathematically, the nonlinear response of the synaptic activity is the key ingredient responsible for the emergence of a stable balanced regime.Our claim is supported by a simple self-consistent analysis accompanied by extensive simulations performed for increasing network sizes.The observed regime is essentially fluctuation driven and characterized by highly irregular spiking dynamics of all neurons.
Neurons in the cortex fire irregularly and with a low firing rate despite being subject to a continuous stimulation (bombardment) from thousands of pre-synaptic neurons.This seemingly counterintuitive evolution has been explained in terms of the theory of dynamical balance of excitation and inhibition [1], considered as one of the major contributions of theoretical physics to neuroscience.However, this theory has been recently criticized, because it requires unphysically large external currents, experimentally unjustified.Here, we propose a nonlinear balancing mechanism based on a biologically plausible form of synaptic plasticity (short-term depression) which works also with weak external currents.
Inferring the collective behavior of large ensembles of oscillators is a highly challenging task, since it requires combining concepts and tools of nonlinear dynamics with those of statistical mechanics [2].Already the classification of proper "thermodynamic" phases and of the conditions for their emergence is a non trivial task: what are the qualitative differences among the various regimes?A preliminary difficulty is posed by the identification of appropriate model classes.Typically (but not exclusively), the coupling is assumed to result from the linear composition of two-body interactions.However, already within this simplified setup, a question arises in the case of massive coupling, when the number K of connections is proportional to the number N of oscillators (mean-field models being the ultimate example).In fact, a meaningful thermodynamic limit requires the coupling term to be finite for K → ∞.This is typically ensured by assuming a coupling constant of order 1/K: we call this type-I coupling and the Kuramoto model is perhaps the most famous example [3].An alternative approach can be adopted, when the single two-body coupling terms are, on average, equal to zero.In such cases, it makes sense to assume that the coupling constant is of order 1/ √ K.We call it type-II coupling and the XY spin-glass model [4,5] is one of the most famous representatives of this class.
The Hamiltonian-mean-field [6] is an enlightening example, which encompasses both options.This extension of the Kuramoto model, if equipped with type-I coupling, proves useful to describe chaotic properties of the synchronized (magnetized) phase [7]; if equipped with type-II coupling, it has helped to discover nontrivial properties of the asynchronous (unmagnetized) regime [8].
In this context, for mean-field models of globally coupled oscillators with type-I coupling, the stationary regime is often found to be asynchronous, i.e. characterized by constant collective features (such as the local field potential) possibly accompanied by tiny fluctuations resulting from the finiteness of the neuronal population.Partial synchrony may manifest itself as either periodic macroscopic oscillations [13,18], or irregular fluctuations [19,20].Anyway, the corresponding single-neuron firing activity is typically regular (the coefficient of variation (CV) of the interspike intervals is small), contrasting the experimental evidence that cortical neurons in vivo operate erratically and with a relatively low firing rate [21] in spite of receiving stimulations from thousands of pre-synaptic neurons [22][23][24].
An irregular firing activity is generated if the neurons operate in the so-called fluctuation-driven regime, when they stay in proximity of the firing threshold, crossed at random times thanks to self-generated fluctuations [25].This can happen when inhibition is strong and accompanied by a random connectivity which suppresses coherence across the neuronal population [26].
Altogether, it is widely accepted that all these features can dynamically emerge whenever the underlying dynamics is in the so-called balanced regime [1], observed in excitatory-inhibitory networks characterized by type-II coupling, an assumption consistent with optogenetic experiments in vitro [27].A balanced state can be, e.g., found, by assuming: (i) a sufficiently large in-degree K; (ii) coupling strengths of order 1/ √ K; (iii) external currents of order √ K [1, [28][29][30][31][32][33][34][35].If the external currents are of order O(1), excitation and inhibition still balance each other, but the firing activity decreases as 1/ √ K, suggesting that strong external currents are a necessary ingredient.However, this latter hypothesis has recently received several criticisms [36,37] based on experimental evidences that the external inputs are O(1) [38][39][40].
In this paper, we show that the introduction of nonlinearities can robustly sustain and stabilize a balanced regime, where the irregular firing of excitatory and inhibitory neurons compensate each other, without neither the inclusion of strong external currents, nor the ad hoc adjustment of parameter values.The nonlinear mechanism, herein invoked is the well known short-term synaptic depression (STD) [41], arising from the finitude of available resources [42].It has been shown that depression has a prevalent effect on excitatory synapses in the visual cortex, inducing dynamical variations of the balance between excitation and inhibition [43].
More precisely, we consider a plastic network of pulsecoupled phase-oscillators, where STD modifies nonlinearly the synaptic inputs.For the sake of simplicity and consistently with experimental indications [43,44], STD is assumed to act only on the synapses connecting excitatory neurons.We show that this little adjustment suffices to ensure the self-sustainement of an irregular activity.
The model We consider two coupled populations each composed of N neurons.The evolution of the membrane potential v e/i j of the j-th neuron within the excitatory/inhibitory population follows from the equation, Whenever v e/i j reaches the threshold 1, it is reset to 0, and, simultaneously, a smooth post-synaptic α-pulse p α (t) = α 2 te −αt is delivered to all the connected neurons, mimicking a non-istantaneous synaptic transmission [12,45,46].For large α-values, p α is well approximated by a δ-pulse [13].F (v) > 0 describes the bare neuron velocity field under the action of a weak constant external current, such that it operates slightly suprathreshold; C e/i j denotes the incoming synaptic recurrent current (see below for its definition); H(v) gauges the impact of the current, which may depend on the value of the membrane potential; finally, G ′ denotes the overall coupling strength.For the sake of simplicity, F (v) and H(v) are taken to be the same for all neurons; the difference between excitatory and inhibitory neurons is encoded in their mutual couplings.
Without loss of generality, F (v) can be assumed to be constant.In fact, if we introduce the new variable ϕ obtained by integrating the ode dϕ/dv = ω/F (v) (under the condition ϕ(0) = 0), Eq. (1) rewrites as φe/i j = ω + GZ(ϕ e/i j )C e/i j . ( where Z(ϕ) = ηωH(v(ϕ))/F (v(ϕ)) and η is a normalization constant suitably chosen to set the maximum of Z(ϕ) equal to one (hence, G = G ′ /η).The value of ω is determined by imposing ϕ(1) = 1 (this condition is equivalent to ω = 1/T , where T is the period of the bare neuron activity).Hence, ϕ is a phase-like variable, while Z(ϕ) can be viewed as an effective phase response curve (PRC) [47,48].
Leaky Integrate-and-Fire (LIF) neurons are among the most popular models used in computational neuroscience (see, e.g.[49]).For current based coupling, they are characterized by the LIF model can be recast in the equivalent representation (15), where Z LIF (ϕ) = e (ϕ−1) (see the blue curve in Fig. 1(a)).
In this paper, we have considered also 6 ] (see the red curve in Fig. 1(a)) for its continuity at threshold, as usually assumed in realistic PRCs [50], and its resemblence to PRC for type I membrane excitability [51].
Finally, the incoming synaptic currents are defined as where the coefficients (g e e , g e i , g i e , g i i ) quantify the specific intra and inter synaptic strengths of excitatory and inhibitory populations, while K e/i is the average in-degree, and I j and E e/i j represent the incoming inhibitory and (effective) excitatory fields.The inhibitory field obeys the differential equation where α is the inverse pulse-width, while t i (n, m) denotes the delivery time of the m-th spike from the n-th inhibitory neuron to the j-th neuron.The sum ′ n is restricted to the K i neighbours of neuron j.This representation amounts to assuming that I j (t) is the linear superposition of the α-pulses received by the neuron j from inhbitory neurons until time t.The excitatory field is treated in a slightly different way, where n identifies the excitatory neuron sending the mth spike to the j-th neuron; x e/i ∈ [0, 1] represents the synaptic efficacy.If the receiving neuron is inhibitory x i ≡ 1, while x e is affected by the STD acting on excitatory-to-excitatory connections.By following [52], its evolution can be written as where t e (n, m) identifies the time of the m-th spike emitted by the n-th neuron itself.Whenever the neuron spikes, the synaptic efficacy x e n is reduced by a factor u, representing the fraction of resources consumed to produce a post-synaptic spike.So long as the n-th excitatory neuron does not spike, the variable x e b increases towards 1 over a time scale τ d .
The in-degrees of the two populations are distributed as in a massively coupled Erdös-Renyi random graph [53], i.e. setting K e/i = p e/i N , where p e/i ∈ [0, 1] is the probability to have a pre-synaptic connection.
Self-consistent approach Before discussing the direct numerical simulations, we present a simple self-consistent approach to explain how STD can actually stabilize a balanced state even in absence of strong external currents.In the above defined setup, E e/i j and I j , being proportional to the in-degree, are also proportional to N , so that the two terms in Eq. ( 3) are both proportional to √ N .It is useful to make this dependence explicit, by writing where represent the average firing rates of the inhibitory (excitatory) pre-synaptic spike trains stimulating the j-th neuron; finally, E e j = E e j /(p e N ) is the effective firing rate of an excitatory neuron scaled to account for the reduced efficacy due to STD.The approximation consists in neglecting neuron-to-neuron fluctuations, as well as temporal variations so that we can drop the j dependence of both the fields and the input currents and assume that they are constant.Within this approximation, a balanced regime can exist if C e/i remains finite for N → ∞, or, equivalently, if the terms in square brackets in Eq. ( 16) converge to 0 (as 1/ √ N ).Accordingly, where the subscript "•" here and elsewhere means that the variable refers to the N → ∞ limit.Hence, In the absence of STD, since E e • = E i • , Eq. ( 9) is satified only when the rightmost h.s. is set equal to 1 from the outset: hence, the balanced regime is highly non-generic.In the literature, a way out is typically found by assuming that the external current ω is of order √ N , so that it must be included in the balance conditions (8) as additional given terms.As a result, the two equations compose a generically solvable, linear, inhomogeneous system [1], the only condition for its validity being that the fields must turn out positive.
In the present context, instead, the novelty is that the ratio between E i and E e is not a priori equal to 1, but depends on the activity of the excitatory neurons.In fact, E e = θE i , where θ is the value of the synaptic efficacy when the excitatory neurons reach threshold during their periodic firing activity.Hence, in the thermodynamic limit, the balance condition requires where the inequality follows from the fact that θ • must, by definition, be smaller than 1.In other words, a balanced regime is generic as it can arise for a broad range of coupling constants.In this paper, since we set g e e = 1, g e i = 1/2, g i e = 1, and g i i = 2, the inequality is satisfied (1/4 < 1).
The equality between the first two terms in (10) allows determining the value of synaptic efficacy and, in turn, of the interspike interval T e • of the excitatory neurons.In fact, from the integration of Eq. ( 6) in between two consecutive spikes, The still unknown initial condition x e (0) can be determined by imposing x e (0) = ux e (T e ) (u is the depletion factor), obtaining By then solving for T e , we find, in the thermodynamic limit, which means that the ISI and thereby the amplitude of the excitatory field E i • = 1/T e • , are independent of the PRC (within this approximation).Finally, the amplitude I • of the inhibitory field is obtained from the first of Eq. (8).Mean Firing Rates The self-consistent analysis is useful to identify the necessary conditions for the onset of a balanced regime, but it unavoidably predicts a currentdriven regime.In order to analyse the actual network behavior it is necessary to perform numerical simulations.
Here below, we report the results for a homogenous network, where the bare firing rate is set to ω = 50 Hz and the PRC is Z I (ϕ) (see [54] for the other parameter values).
In Fig. 1(b) we plot the population firing rates of excitatory ν e = ⟨E i j ⟩ and inhibitory ν i = ⟨I j ⟩ neurons (⟨•⟩ denotes an average over all neurons of a given population) versus N .The data are well fitted by the law ν e/i 0 + µ e/i / √ N , with ν e 0 ≃ ν i 0 ≃ 5.78 Hz, µ e ≃ 399 Hz, and µ i ≃ 762 Hz (curves not shown).This indicates that the single-neuron activity remains finite for N → ∞, a clear signature of a balanced regime.This conclusion is confirmed by the N -dependence of the average unbalance ∆ e/i = ⟨∆ e/i j ⟩, reported in Fig. 1(c), where a clear 1/ √ N decrease is visible, implying that the average value of C e/i stays constant for N → ∞.
It is instructive to compare the numerical results with the semi-analytical perturbative implementation of the self-consistent approximation.In the previous section, we have shown how to determine the values of the excitatory and inhibitory fields for N → ∞.For our parameter values, I • = E i • = 1/ log(7/6)Hz ≃ 6.487 Hz [55].In the Supplementary Material, we show how to go beyond the asymptotic values, determining finite-size corrections.Here, we sketch the procedure.From the knowledge of the asymptotic fields, one can determine the input currents C e/i • responsible for those fields, by integrating Eq. ( 15) under the assumption of a constant current.Then, the definition (16) of C e/i can be used as a consistency relationship to determine the finite-size corrections for both fields, which turn out to be proportional to 1/ √ N .The resulting analytic expressions are reported in the Supplementary Material and plotted in Fig. 1(b) (see the solid curves).They overestimate the numerical values, but are not too far from them.
Fluctuations We start investigating whether the collective dynamics of the network remains asynchronous even for large system sizes, as expected in brain circuits [28].The raster plot for N = 16000, reported in the inset of Fig. 2(a), does not reveal any population oscillation.A more quantitative analysis has been made by computing the time-average of the standard deviation of the incoming fields (i.e. of the instantaneous population firing rates), here denoted with σ e/i .The values computed for different network sizes, reported in Fig. 2(a), decrease consistently with the 1/ √ N scaling expected from the central limit theorem for an asynchronous dynamics.
Next, we focus on temporal fluctuations at the singleneuron level.They are disregarded a priori by the self-consistent approach, but the probability distribution density (PDF) of the coefficient of variations CV j reported in Fig. 2(b) for the excitatory neurons gives a clear evidence of irregularity.Some neurons are charac-terized by a CV even larger than 1, the value expected for a Poisson distribution, and the irregularity tends to increase with N .A similar scenario is exhibited by inhibitory neurons (data not shown).
Finally, we turn our attention to ensemble fluctuations.The firing rates themselves are broadly distributed from nearly vanishing values (almost silent neurons) up to 50-60 Hz, with a pronounced peak around 5-10 Hz.When N is increased, the PDF widths remain finite and their shapes appear to converge to a given asymptotic form, as clearly visible in Fig. 2(c) where the data refer to excitatory neurons.This manifestation of heterogeneity is not surprising in a massively coupled Erdös-Renyi network.In fact, the single-neuron connectivity is expected to exhibit fluctuations of order √ N , which transform themselves into fluctuations of O(1) for the C e/i j , and therefore for the firing rates.
The distribution of firing rates ν e j induces a distribution of synaptic efficacies θ j (taken in correspondence of the spiking times).Under the approximation of negligible temporal fluctuations, one can reformulate Eq. ( 12) as The inset in Fig. 2(c) reveals a good agreement with the numerical simulations.
The PDF shapes reported in Fig. 2(c), are similar to those measured experimentally in the cortex and hippocampus [56][57][58][59][60], with many neurons exhibiting a low firing rate and a high-frequency long tail, akin to a log-normal distribution.These shapes are typically interpreted as an indication of fluctuation-driven dynamics [61].It is therefore convenient to test whether the neurons, in our model, operate either above or below threshold.This can be done as follows.From Eq. ( 15), since the maximum of Z(v e/i ) is 1, ve/i may have a stable zero, only if ω + GC e/i < 0. Hence, a neuron characterized by an average current C e/i is typically below threshold if C e/i < −ω/G = −50 Hz.The data reported in Fig. 2(d) reveal a mixed behavior: depending on their effective firing rate, neurons may be either fluctuationor current-driven.By further averaging over the entire population, we find that while the inhibitory neurons are significantly fluctuation-driven with ⟨C i ⟩ = −101.5Hz, on average the excitatory neurons operate slightly below threshold being ⟨C e ⟩ = −55.5 Hz.Altogether, the network self-stabilizes in a regime, where fluctuations play a major role, consistently with the observation of a pseudo log-normal distribution of the firing rates [61].
Robustness of the mechanism Additional simulations performed for different parameter values and even introducing heterogeneity in the external currents have shown the generality of the mechanism.Details can be found in the Supplementary Material.Here we focus  b-c) different system sizes coded as in Fig. 2 (b-c).The results here reported refer to the LIF model.
on the most important test, made by using the PRC of LIF neurons, Z LIF .All parameters have been left unchanged, except for a faster synaptic transmission α [54] to check the specificity of the pulsewidth.The simulations are performed by integrating Eq. ( 15) with Z LIF .As shown analytically (and verified numerically) this model is fully equivalent to a standard current-driven spiking LIF network (1) with a = e/(e − 1) ≃ 1.582 and G ′ = 1/(e − 1) ≃ 0.582.In Fig. 3 (a) we report the firing rates of the two populations (black and red dots), together with the outcome of the self-consistent approach (solid curves).The theoretical curves converge, as they should, to the same value 1/T • , which coincides with the previous asymptotic value since, from its definition, it does not depend on the PRC shape.Also the numerically determined firing rates converge to the same value, again following the law ν e/i ≃ ν 0 + µ e/i / √ N , where ν 0 ≃ 5.72 Hz, µ e ≃ 480 Hz and, and µ i = 803 Hz (curves not shown).The smallness of the discrepancy with the previous asymptotic value (≃ 1%) indicates that the irrelevance of the PRC extends to the complete model, where all kinds of fluctuations are automatically included.
In Fig. 3 we also report the PDF of the firing rates of inhibitory neurons ν i j (panel b) and of the corresponding coefficients of variations CV j (panel c) for different network sizes.The distributions of the firing rates display a long tail towards vanishing rates, a typical feature of neurons which operate below threshold and that are therefore fluctuation driven, as confirmed by the values of the corresponding CV j .Analogous results have been found for the excitatory neurons (data not shown).We can safely state that the data obtained for the LIF network confirm the scenario of a balanced regime as for the Z I (ϕ) PRC.

Conclusions
The typical setup studied in the literature to discuss dynamically balanced regimes requires the presence of strong external inputs [1].An alternative layout, which does not suffer this limitation, was proposed in Ref. [37], together with the concept of sparse balance.It, however, requires an anomalously broad distribution of synaptic strengths and leads to a vanishing fraction of active neurons (in the thermodynamic limit).
The mechanism proposed here is more robust and general: it exploits the dynamical adjustment of the synaptic currents, resulting from short-term synaptic depression (STD).STD is a much studied mechanism, already invoked to explain fundamental cognitive functions, such as working memory [62][63][64] and the internal representation of spatio-temporal information [65][66][67][68].By lowering the strength of highly active excitatory connections, STD eventually binds the activity of excitatory neurons.
Mathematically, the balanced regime is made possible by the nonlinear dependence of one excitatory current on the amplitude of the corresponding field.In practice, the homogeneous set of linear conditions (8), which emerges for weak external currents, is transformed into a nonlinear one.While the former one admits a meaningful solution only for a special combination of the coupling constant, the latter is generically solvable for a broad range of parameter values.Once this has been understood, it is straightforward to infer that other nonlinear mechanisms may play the same role as STD in the absence of strong inputs.In fact, other sources of naturally expected nonlinearities have been recently investigated in computational neuroscience although, always in the presence of strong external currents.Spike-frequency adaptation is one such mechanism, studied in networks with highly heterogenous in-degrees [69].Similarly, facilitation has been found to promote the emergence of bistable balanced regimes [70,71].It is time to move on and to investigate all such mechanisms in the absence of strong external currents.

Self-Consistent Analysis
Here we develop a perturbative analysis to determine the mean field properties of the asynchronous regime, for a homogeneous setup, via a self-consistent approach which neglects neuron-to-neuron fluctuations.Within this approximation, the neuron label j can be dropped into the corresponding evolution equations φe/i = ω + GZ(ϕ e/i )C e/i ; (15) where and the field E e/i (I) represents the mean firing rate of the pre-synaptic excitatory (inhibitory) neurons, which, for an asynchronous regime, is constant.In this case the firing period is In the paper, we have shown how to determine the asymptotic values (i.e. for N → ∞) of all the fields.By integrating Eq. ( 15 For the LIF model we simply have where v = exp ϕ and Y e/i = (ωe)/(C For the values selected in our paper [54], T e/i • = ln 7/6 = 0.15415 s, so that .487 Hz and these values do not depend on the chosen PRC.By inserting our parameter values, we find that Y e/i = −2.7204,which corresponds to C e/i • = −49.96Hz.Hence, in the LIF case, the neurons families operate (in this mean field approximation) even closer to threshold.
Now, we can focus on Eq. ( 16) to obtain the first-order corrections to the fields since C e • and C i • are known, β e e θE i − β e i I = where we have exploited the equality E e = θE i , θ being the value of the synaptic efficacy when the excitatory neurons reach threshold during their periodic firing activity.The value of θ can be determined from Eq. ( 12) in the main text, here rewritten as with T = T e /τ d , to simplify the notations.By solving the equation ( 20) for and replacing it in (21), By employing the explicit expression (22) for θ, the above equation can be rewritten as and hence, where is the asymptotic (N → ∞) θ value as defined in Eq. ( 10) in the main text.Now, we expand At first order, it is sufficent to expand the numerator of the l.h.s.
and finally The perturbative expression for I is and where we have set C where the first order corrections are slightly larger than in the case of Z I (ϕ).

Heterogeneous Network and Further Dynamical Regimes
In the paper, we have investigated a homogeneous setup, where all neurons are characterized by the same bare frequency, but this is definitely an idealization of reality.In order to further validate the robustness of the proposed mechanism, we have also studied a model where the bare frequencies are uniformly distributed over two different ranges (for the excitatory and inhibitory neurons).A balanced regime emerges again, where the network-induced heterogeneity adds up to the spontaneous single-neuron heterogeneity.In Fig. 4(a), we plot the difference between the effective ν i j and the bare frequency ω i j , versus the bare frequency itself (for each inhibitory neuron) to highlight the impact of the network coupling.An overall inhibitory effect is noticeable, which increases upon increasing the bare frequency.
We have also performed various simulations for different sets of parameter values, confirming that an asynchronous balanced state is a generic regime in networks accompanied by STD.The only qualitative change we have found is the onset of collective oscillations, typically when some time scales are varied.As an example, in Fig. 4(b) one can see relatively wide oscillations of both the excitatory and inhibitory fields (see the red dots), to be confronted with the small blue spot which corresponds to the asynchronous regime studied in detail in this Letter.The structure of the attractor is suggestive of an underlying possibly low-dimensional dynamics, but a simulation performed for a single (actually a few) networks sizes does not suffice to conclude whether the fatness of the attractor is either a finite-size effect, or the manifestation of an intrinsic high-dimensionality.

FIG. 1 .
FIG. 1.(a) The PRC Z(ϕ) vs the phase-like variable ϕ for ZI (red curve) and ZLIF (blue curve); (b) the population firing rates ν e = ⟨E i j ⟩ and ν i = ⟨Ij⟩ versus N ; (c) the average unbalance |∆ e/i | versus N , the blue dashed line denotes a power law decay as 1/ √ N .In (b-c) the black (red) color refers to excitatory (inhibitory) neurons and the solid lines in (b) to the self-consistent approximations.The results in (b-c) refer to ZI.

FIG. 2 .
FIG. 2. (a) Standard deviation σ e/i of the population firing rates versus N , the blue dashed line indicates a power law decay as 1/ √ N .The inset displays a raster plot for N = 16000 over a time window of 20 ms.PDF of the coefficients of variation CVi (b) and of the firing rates ν e j (c) for excitatory neurons.In the inset of (c), the synaptic efficacy at threshold θj is displayed versus the corresponding firing rate ν e j : the dashed line is the self-consistent estimation (14).(d) Time averaged values C e j (C i j ) versus the corresponding firing rates ν e j (ν i j ) for N = 16000.The dashed line denotes the value C e/i = −50 Hz discriminating fluctuation from current driven dynamics.The black (red) color in (a-d) refers to excitatory (inhibitory) neurons, while the colors in (b-c) to different system sizes : namely, N = 8000 (black); N = 16000 (red); N = 32000 (green); N = 64000 (blue) and N = 128000 (magenta).The reported data refer to ZI.

FIG. 3 .
FIG. 3. (a)The average population firing rates ν e (black circles) and ν i (red circles) versus N ; PDF of the inhibitory firing rates ν i j (b) and of the corresponding CVj (c).In (a) the black (red) solid line indicates the self-consistent approximation for excitatory (inhibitory) neurons and the colors in (b-c) different system sizes coded as in Fig.2 (b-c).The results here reported refer to the LIF model.
), we can now determine C e/i • as as those values which yield the expected firing period T e/i • (i.e. the expected field).For the PRC Z I (ϕ) we find C e • = C i • = −49.108Hz, i.e. both families of neurons operate slightly above threshold since we set ω = 50 Hz.

=
T • , since they coincide for our choice of the parameters.For Z I (ϕ) we have found numerically T • and C • and therefore also the values for the following perturbative expansion to the first order:E i = 6.487+643.61/√ N , I = 6.487+817.23/√ N ,where E i and I are expressed in Hz.For the LIF model we have instead found analytically that :

FIG. 4 .
FIG. 4. (a) Difference between the single neuron firing rates and their bare fequencies ν i j − ω i j versus ω i j for an heterogeneous network with N = 16000 neurons with ω e j (ω i j ) uniformly distributed in the interval [15, 65] Hz ([35, 85] Hz).(b) Inhibitory fields Ij versus the corresponding excitatory field E i j for a homogeneous network with N = 16000 and τ d = 1 s (blue dots) and τ d = 0.1 s (red dots).The values of the fields for τ d = 1 s ar arbitrarily shifted for a better visualization.The reported results refer to the PRC ZI(ϕ).