Collective dynamics in the presence of finite-width pulses

The idealisation of neuronal pulses as δ-spikes is a convenient approach in neuroscience but can sometimes lead to erroneous conclusions. We investigate the effect of a finite pulse-width on the dynamics of balanced neuronal networks. In particular, we study two populations of identical excitatory and inhibitory neurons in a random network of phase oscillators coupled through exponential pulses with different widths. We consider three coupling functions, inspired by leaky integrate-and-fire neurons with delay and type-I phase-response curves. By exploring the role of the pulse-widths for different coupling strengths we find a robust collective irregular dynamics, which collapses onto a fully synchronous regime if the inhibitory pulses are sufficiently wider than the excitatory ones. The transition to synchrony is accompanied by hysteretic phenomena (i.e. the co-existence of collective irregular and synchronous dynamics). Our numerical results are supported by a detailed scaling and stability analysis of the fully synchronous solution. A conjectured first-order phase transition emerging for δ-spikes is smoothed out for finite-width pulses.

The idealisation of neuronal pulses as δ -spikes is a convenient approach in neuroscience but can sometimes lead to erroneous conclusions.We investigate the effect of a finite pulse-width on the dynamics of balanced neuronal networks.
In particular, we study two populations of identical excitatory and inhibitory neurons in a random network of phase oscillators coupled through exponential pulses with different widths.We consider three coupling functions, inspired by leaky integrate-and-fire neurons with delay and type-I phase-response curves.By exploring the role of the pulsewidths for different coupling strengths we find a robust collective irregular dynamics, which collapses onto a fully synchronous regime if the inhibitory pulses are sufficiently wider than the excitatory ones.The transition to synchrony is accompanied by hysteretic phenomena (i.e. the co-existence of collective irregular and synchronous dynamics).Our numerical results are supported by a detailed scaling and stability analysis of the fully synchronous solution.A conjectured first-order phase transition emerging for δ -spikes is smoothed out for finite-width pulses.
Neuronal networks with a nearly balanced excitatory/inhibitory activity evoke significant interest in neuroscience due to the resulting emergence of strong fluctuations akin to those observed in the resting state of the mammalian brain.While most studies are limited to a δ -like pulse setup, much less is known about the collective behavior in the presence of finite pulse-widths.In this paper, we investigate exponential pulses, with the goal of testing the robustness of previously identified regimes such as the spontaneous emergence of collective irregular dynamics (CID), an instance of partial synchrony with a non-periodic macroscopic dynamics.Moreover, the finitewidth assumption paves the way to the investigation of a new ingredient, present in real neuronal networks: the asymmetry between excitatory and inhibitory pulses.Our numerical studies confirm the emergence of CID also in the presence of finite pulse-width, although with a couple of warnings: (i) the amplitude of the collective fluctuations decreases significantly when the pulse-width is comparable to the interspike interval; (ii) CID collapses onto a fully synchronous regime when the inhibitory pulses are sufficient longer than the excitatory ones.Both restrictions are compatible with the recorded behavior of real neurons.Additionally, we find that a seemingly first-order phase transition to a (quasi)-synchronous regime disappears in the presence of a finite width, confirming the peculiarity of the δ -spikes.A transition to synchrony is instead observed upon increasing the ratio between the width of inhibitory and excitatory pulses: this transition is accompanied by a hysteretic region, which shrinks upon increasing the network size.Interestingly, for a connectivity comparable to that of the mammalian brain, such a finite-size effect is still sizable.Our numerical studies might help to understand abnormal synchronisation in neurological disorders.

I. INTRODUCTION
Irregular firing activity is a robust phenomenon observed in certain areas of mammalian brain, such as hippocampus or cortical neurons 1,2 .It plays a key role for the brain functioning in the visual and prefrontal cortex.This behavior emerges from the combined action of many interacting units 3,4 .
This paper focuses on a regime called collective irregular dynamics (CID), which arises in networks of oscillators (neurons).Mathematically, CID is a non-trivial form of partial synchrony.Like partial synchrony, it means that the order parameter χ used to identify synchronization (see Sec. II for a precise definition) is strictly larger than 0 and smaller than 1.Moreover, it implies a stochastic like behavior of macroscopic observables such as the average membrane potential.
There are (at least) two mechanisms leading to CID: (i) the intrinsic infinite dimensionality of the nonlinear equations describing whole populations of oscillators; (ii) an imperfect balance between excitatory and inhibitory activity.
Within the former framework, no truly complex collective dynamics can arise in mean-field models of identical oscillators of Kuramoto type.In fact, the Ott-Antonsen Ansatz 5 implies a strong dimension reduction of the original equations.Nevertheless, in this and similar contexts, CID can arise either in the presence of a delayed feedback 6 , or when two interacting populations are considered 7 .Alternatively, it is sufficient to consider either ensembles of heterogeneous oscillators: e.g., leaky integrate-and-fire (LIF) neurons 8 , and pulse-coupled phase oscillators 9 (notice that in these cases, Ott-Antonsen Ansatz does not apply).
Within the latter framework, an irregular activity was first observed and described in networks of binary units, as a consequence of a (statistical) balance between excitation and inhibition 10 .This balanced regime 11 can be seen as an asynchronous state accompanied by statistical fluctuations.In fact, this interpretation led Brunel 12 to develop a powerful analytical method based on a self-consistent Fokker-Planck equation to describe an ensemble of LIF neurons.In the typical (sparse) setups considered in the literature, the fluctuations of the single neuron activity vanish when averaged over the whole pop-ulation, testifying to their statistical independence; in terms of order parameter, χ = 0.
However, it has been recently shown that a truly CID can be observed in the presence of massive coupling (finite connectivity-density) under the condition of small unbalance 13,14 .In this paper we test the robustness of these results, obtained while dealing with δ -pulses, by studying more realistic finite-width pulses.In fact, real pulses have a small but finite width 15 .Moreover, it has been shown that the stability of some synchronous regimes of LIF neurons may qualitatively change, when arbitrarily short pulses are considered (in the thermodynamic limit) 16 .
A preliminary study has been already published in Ref. [17], where the authors have not performed any finitesize scaling analysis and, more important, no any test of the presence of CID has been carried out.Here we study a system composed of two populations of (identical) excitatory and inhibitory neurons, which interact via exponential pulses of different width, as it happens in real neurons 18 .
Handling pulses with a finite width requires two additional variables per single neuron, in order to describe the evolution of the incoming excitatory and inhibitory fields.The corresponding mathematical setup has been recently studied in Ref. [19] with the goal of determining the stability of the fully synchronous state in a sparse network.The presence of two different pulse-widths leads to non-intuitive stability properties, because the different time dependence of the two pulses may change the excitatory/inhibitory character of the overall field perceived by each single neuron.Here, we basically follow the same setup introduced in Ref. [19] with the main difference of a massively coupled network, to be able to perform a comparative analysis of CID.
The randomness of the network accompanied by the presence of three variables per neuron, makes an analytical treatment quite challenging.For this reason we limit ourselves to a numerical analysis.However, we accompany our studies with a careful finite-size scaling to extrapolate the behavior of more realistic (larger) networks.Our first result is that CID is observed also in the presence of finite pulse-width, although we also find a transition to full synchrony when the inhibitory pulses are sufficiently longer than excitatory ones.The transition is first-order (discontinuous) and is accompanied by hysteresis: there exists a finite range of pulse-widths where CID and synchrony coexist.
The finite-size analysis suggests that in the thermodynamic limit CID is not stable when the pulses emitted by inhibitory neurons are strictly longer than those emitted by the excitatory ones.However, the convergence is rather slow and we cannot exclude that the asymmetry plays an important role in real neuronal networks of finite size.
More precisely in section II, we define the model, including the phase response curves (PRCs) used in our numerical simulation.In the same section we also introduce the tools and indicators later used to characterize the dynamical regimes, notably an order parameter to quantify the degree of synchronization 20 .In section III, we present some results obtained for strictly δ pulses to test robustness of CID in our context of coupled phase oscillators.In Sec.IV we discuss the symmet-ric cases of identical finite pulse-widths.Sec.V is devoted to a thorough analysis of CID by varying the pulse-widths.Sec.VI contains a discussion of the transition region, intermediate between standard CID and full synchrony.In the same section, the robustness of the transition region is analysed, by considering different PRCs.Finally, section VII is devoted to the conclusions and a brief survey of the open problems.

II. MODEL
Our object of study is a network of N phase oscillators (also referred to as neurons), the first N e = bN being excitatory, the last Each neuron is characterized by the phase-like variable Φ j ≤ 1 (formally equivalent to a membrane potential), while the (directed) synaptic connections are represented by the connectivity matrix G with entries where ∑ N e k=1 G j,k = K e and ∑ N k=N e +1 G j,k = K i , meaning that each neuron j is characterized by the same number of incoming excitatory and inhibitory connections, as customary assumed in the literature 21 (K = K e + K i represents the connectivity altogether).Here, we assume that K is proportional to N, that is K = cN, i.e. we refer to massive connectivity.Further, the network structure is without autapse, i.e.G j, j = 0.
The evolution of the phase of both excitatory and inhibitory neurons is ruled by the same equation, where E j (I j ) the excitatory (inhibitory) field perceived by the jth neuron, while Γ(Φ) represents the phase-response curve (PRC) assumed equal for all neurons; finally, µ is the coupling strength.Whenever Φ k reaches the threshold Φ th = 1, it is reset to the value Φ r = 0 and enters a refractory period t r during which it stands still and is insensitive to the action of both fields.The fields E j and I j are the linear superposition of exponential spikes emitted by the upstream neurons.Mathematically, where α (β ) denotes the inverse pulse-width of the excitatory (inhibitory) spikes and t k n is the emission time of the nth spike emitted by the kth neuron.The coefficient g accounts for the relative strength of inhibition compared to excitation.If the kth neuron is excitatory, P k = 1, otherwise P k = 0.
pressed as simple sums Let us now introduce the PRCs used later in our numerical simulations.We consider three different shapes: • PRC 3 The various curves are plotted in Fig. 1.PRC 1 (see the black curve, which corresponds to Φ L = −0.1 and Φ U = 0.9) has been introduced in Ref. [19] to study the stability of the synchronous regime; its shape has been proposed to mimic a network of leaky integrate-and-fire neurons in the presence of delay (see also Ref. [9]).
The two other PRCs have been selected so as to explore the effect of a progressive regularization of the neuronal response.In particular, we consider the smooth PRC 3 (see Eq. ( 6)), as a prototype of type I PRC 22,23 .
The network dynamics is simulated by implementing the Euler algorithm with a time step δ t = 10 −3 .However, in some cases, smaller steps have been considered to avoid spurious synchronization.We typically initialize the phases uniformly in the unit interval, while the fields are initially set equal to zero.
In all cases, we have set b = 0.8, c = 0.1 and g = 4 + 1000/K (following the existing literature 13 ).The last condition ensures that the balanced regime is maintained for K, N → ∞.Moreover, we have systematically explored the role of α and β , as the pulse-width is the focal point of this paper.Additionally, the coupling strength µ has been varied, as well as the network-size N, to test for the amplitude of finite-size effects.
The following statistical quantities are used to characterize the emerging dynamical states.
1.The mean firing rate is a widely used indicator to quantify the neural activity.It is defined as where N j (t) denotes the number of spikes emitted by the neuron j over a time t.
2. The coefficient of variations C v is a microscopic measure of irregularity of the dynamics based on the fluctuations of the interspike intervals (ISIs).The average C v is defined as where σ j is the standard deviation of the singleoscillator's ISI, and τ j is the corresponding mean ISI.
If C v > 1, then the neurons show a bursting activity, while C v < 1 means that the spike train is relatively regular.
3. The order parameter, χ, is typically used to quantify the degree of synchronisation of a population of neurons 24 .
It is defined as where • represents an ensemble average, while the over-bar is a time average.The numerator is the variance of the ensemble average Φ , while the denominator is the ensemble mean of the single-neuron's variances.When all neurons behave in exactly the same way (perfect synchronization), then χ = 1.If instead, they are independent, then χ ≈ 1/ √ N. Regimes characterized by intermediate finite values 0 < χ < 1 are referred to as instances of partial synchronization.However, χ > 0 does not necessarily imply that the collective dynamics is irregular: it is, e.g., compatible with a periodic evolution.In fact, here we report several power spectra to testify the stochastic-like dynamics of macroscopic (average) observables.

III. DELTA PULSE
Most spiking network models deal with δ -spikes, including those giving rise to CID 13,14 .This paper is focused on the more realistic exponential spikes, but before proceeding in that direction we wish to briefly discuss the case of zero pulse-width.This is useful to gauge the different PRCs used in this paper.Since δ pulses correspond to the limiting case α, β → ∞, they can be treated by invoking Eq. (3). Figure 2 shows the various indicators introduced in Section II, to characterize the collective dynamics.As in previous papers, 13,14 we explore the parameter space, by varying the coupling strength µ and the system size N.
In panel (c) we can appreciate that CID emerges already for very small coupling strength; it is accompanied by an increasing average coefficient of variations C v , due to the coupling which induces increasing deviations from purely periodic behavior.In parallel, the mean firing rate ν decreases as a result of the prevalent inhibitory character of the network.This weak-coupling emergence of CID is comparable to what observed in balanced LIF models with δ spikes 13 .
Above µ ≈ 0.537 ≡ µ c (see the vertical dashed lines), a transition occurs towards a highly synchronous regime (χ is slightly smaller than 1), accompanied by a larger firing rate.The corresponding firing activity is mildly irregular: C v is smaller than in Poisson processes (when C v = 1).A quick analysis suggests that this self-sustained regime emerges from the vanishing width of the pulses combined with the PRC shape, which is strictly equal to zero in a finite phase range below the threshold Φ th = 1.In fact, similar studies performed with PRC 3 do not reveal any evidence of a phase transition (see orange stars and green squares in Fig. 2) indicating that such behavior is nothing else but a peculiarity of PRC 1 with δpulses.We have not further explored this regime.It is nevertheless worth noting that the sudden increase of the firing rate observed when passing to the strong coupling regime is reminiscent of the growth observed in LIF neurons 21 , although in such a case, the increase is accompanied by a significantly bursty behavior 25 .
More important is the outcome of the finite-size scaling analysis, performed to investigate the robustness of the observed scenario.In Fig. 2 one can see that the various indicators under stimulation of PRC 1 are size-independent deeply within the two dynamical phases, while appreciable deviations are observed in the transition region.This is customary when dealing with phase-transitions.It is not easy to conclude whether the transition is either first or second order: the C v is reminiscent of the divergence of susceptibility seen in continuous transitions, but this is an issue that would require additional work to be assessed.

IV. IDENTICAL FINITE-WIDTH PULSES
In this section, we start our analysis of finite pulses, by assuming the same width for inhibitory and excitatory neurons, i.e. α −1 = β −1 .The asymmetric case is discussed in the next section.All other system parameters are kept the same as in the previous section (including the PRC shape).
Before discussing the macroscopic measures, we turn our attention to typical CID features.The average phase Φ (t) = 1 N ∑ j Φ j (t) (see Fig. 3(a)) exhibits stochastic-like oscillations, which represent a first evidence of a non-trivial collective dynamics.The raster plot presented in panel Fig. 3(b) contains the firing times t n of a subset of 100 neurons: there, one can easily spot the time intervals characterized by a more coordinated action (see, for instance, around the vertical green line at time 8374 in Fig. 3(a)).A more quantitative representation is presented in Fig. 3(c), where the instantaneous phase distribution P(Φ) is plotted at two different times in correspondence of qualitatively different regimes of the phase dynamics (see the vertical lines in panel (a)).The peak at Φ = 0 is due to the finite fraction of neurons standing still in the refractory period.A small amount of negative phases are also seen: they are due to prevalence of inhibition over excitation at the end of refractoriness.Moreover, the instantaneous phase distribution P(Φ) presented in Fig. 3(c), show that, at variance with the classical asynchronous regime, the shape of the probability density changes with time.The narrowest distribution (green curve) corresponds to the region where strong regular oscillations of Φ are visible in panel (a): within this time interval a "cloud" of neurons homogeneously oscillates from reset to threshold and back.
The resulting order parameter is reported in Fig. 4. In panel (a) we plot χ as a function of µ for different widths: from broad pulses (red stars correspond to α = 1, a width comparable to the ISI), down to very short ones (green triangles correspond to α = 1000).The general message is that partial synchrony is preserved.Nevertheless, it is also evident that increasing the width progressively decreases the ampli- tude of the order parameter.The main qualitative difference is the smoothening of the transition observed for δ -pulses (in correspondence of the vertical dashed line at µ c ).The singular behavior of δ -spikes is confirmed by the relatively large deviations appearing already for α = 1000.
A more direct illustration of the role of α is presented in Fig. 4(b), where we plot χ versus α for different coupling strengths: µ = 0.2 (black triangles), 0.47 (red crosses), and 0.95 (blue diamonds).An overall increasing trend upon shortening the pulse-width is visible for all coupling strengths, although the rate is relatively modest for weak coupling, becoming more substantial in the strong-coupling limit.
Finally, we have briefly investigated the presence of finitesize effects, by performing some simulations for N = 40000 (to be compared with N = 10000 used in the previous simulations): see magenta circles in both panels.We can safely conclude that the overall scenario is insensitive to the network size.

V. FULL SETUP
In the previous section we have seen that the finite width of the spikes does not kill the spontaneous emergence of CID.Here, we analyse the role of an additional element: the asymmetry between inhibitory and excitatory pulses.We proceed by exploring the two-dimensional parameter space spanned by the coupling strength µ and the asymmetry between pulse widths.The latter parameter dependence is explored by setting α = 100 and letting β (the inverse width of inhibitory pulses) vary.All other network parameters, including the PRC shape, are assumed to be the same as in the previous section.
The microscopic manifestation of CID in the setup with non-identical pulses is qualitatively the same as for identical pulses shown in Fig. 3.The results of a systematic numerical analysis are plotted in Fig. 5, where we report three indica-tors: the firing rate ν, the mean coefficient of variation C v , and the order parameter χ, versus β for three different coupling strengths (see the different columns), and four network sizes.
All indicators reveal the existence of two distinct phases: a synchronous regime arising for small β values, and CID observed beyond a critical point which depends on the network size: the transition is discontinuous.All panels reveal a substantial independence of the network size, with the exception of the transition between them (we further comment on this issue later in this section).
The first regime is synchronous and periodic, as signalled by χ = 1, and C v = 0.The corresponding firing rate ν is a bit smaller than 0.97, the rate of uncoupled neurons (taking into account refractoriness).This is consistent with the expected predominance of inhibition over excitation in this type of setup.A closer look shows that in the synchronous regime ν increases with β .This makes sense since the smaller β , the longer the time when inhibition prevails thereby decreasing the network spiking activity.The weak dependence of ν on the coupling strength µ is a consequence of small effective fields felt by neurons when the PRC is small.Finally, for intermediate β values (around 80) and large coupling strengths, χ is large but clearly smaller than 1.This third type of regime will be discussed in the next section.
CID is characterized by a significantly smaller order parameter which, generally tends to increase with the coupling strength.CID is also characterized by a significantly smaller firing rate.This is due the prevalence of inhibition which is not diminished by the refractoriness as in the synchronous regime.Finally, the coefficient of variation is strictly larger than 0, but significantly smaller than 1 (the value of Poisson processes) revealing a limited irregularity of the microscopic dynamics.In agreement with our previous observations for δ -spikes, C v increases with the coupling strength.
Our finite-size scaling analysis also shows that the degree of asymmetry (between pulse widths) compatible with CID progressively reduces upon increasing the number of neurons.
Although the N-dependence varies significantly with the coupling strength, it is natural to conjecture that, asymptotically, CID survives only for β ≥ α.This is not too surprising from the point of view of self-sustained balanced states.They are expected to survive only when inhibition and excitation compensate each other: the presence of different time scales makes it difficult, if not impossible to ensure a steady balance.Transition to synchrony upon lowering β was already observed in Ref. [17] in a numerical study of LIF neurons, where, however, no finite scaling analysis was performed.Interestingly, the onset of a synchronous activity when inhibition is slower than excitation is also consistent with experimental observations 26 .
We conclude this section with a more quantitative characterization of the irregularity of the collective dynamics.In Fig. 6, we plot the Fourier power spectrum S( f ) obtained from Φ (t).The panels correspond to three different coupling strengths (µ = 0.3, 0.47 and 0.95, from top to bottom).For each value of µ, we have sampled three different pulsewidths.
Altogether, one can notice a general increase of the power with µ.This is quite intuitive, as CID is the result of mutual interactions.A less obvious phenomenon is the increase of the power observed when the inhibitory pulse-width β −1 is increased.This is an early signature of a transition towards full synchronization, observed when β is decreased below a critical value.Anyway, the most important message conveyed by Fig. 6 is that all spectra exhibit a broadband structure, although most of the power is concentrated around a specific frequency: f ≈ 1.5 (panel a), f ≈ 1.4 (panel b), and f ≈ 0.93 (panel c).As a result, one can preliminarily conclude that the underlying macroscopic evolution is stochastic-like.A more detailed analysis could be performed by computing macroscopic Lyapunov exponents, but this is an utterly difficult task, as it is not even clear what a kind of equation one should refer to.
Additional evidence of the robustness of CID is given in Fig. 7, where we investigate the amplitude of finite-size corrections, by computing the power spectrum S( f ) for different network sizes for three different parameter sets: µ = 0.3, β = 90 (panel a), µ = 0.3, β = 120 (panel b), and µ = 0.95, β = 95 (panel c).In all cases, the spectra are substantially independent of the number of neurons, although in panel (b) we observe a weak decrease in the band f ∈ [1, 2.5], while a new set of peaks is born in panel (c).Since the connectivity K of the largest networks herein considered (N = 80000) is comparable to that of the mammalian brain (K = 8000 vs 10000) 4 , we can at least conjecture that this phenomenon may have some relevance in realistic conditions.
Finally, the low frequency peak clearly visible for small µ coincides with the mean firing rate (see Fig. 5(a)), while the connection with the microscopic firing rate is lost in panel (c).

VI. TRANSITION REGION
In Fig. 5 we have seen a clear evidence of a first-order phase transition, when either the pulse-width or the coupling strength is varied.So far, each simulation has been performed by selecting afresh a network structure.The stability of our results indicates that the transition does not suffer appreciable sample-to-sample fluctuations.
The main outcome of our numerical simulations is summarized in Fig. 8; the various lines identify the transition between the two regimes, for different network sizes.The critical points have been determined by progressively decreasing β (see Fig. 5) and thereby determining the minimum β -value where CID is stable.Upon increasing N, the synchronization region grows and the transition moves towards β = α.
So far, the initial condition has been chosen by selecting independent, identically uniformly distributed random phases and zero fields.Since it is known that discontinuous transitions are often accompanied by hysteretic phenomena, we now explore this possibility.We start fixing a different type of initial conditions: the phases are selected within a small interval of width δ p (while the fields are set equal to zero and δ t = 10 −4 ). 27Fig. 9 combines the scenario presented in the previous section for a network with N = 10000 neurons and µ = 0.3 (the blue dots correspond to the content of Fig. 5A), with the results of the new simulations obtained for δ p = 10 −3 (see black curves and triangles). is a clear bistability: the new simulations reveal that χ ≈ 1, much above the typical CID value.More precisely, χ < 1 for β ∈ I 2 ≈ [46, 91], while χ = 1 for β ∈ (I 1 − I 2 ).Since χ = 1 is accompanied by a vanishing C v , it is straightforward to conclude that this regime is the periodic synchronous state, whose linear stability can be assessed quantitatively.
The conditional Lyapunov exponent λ c provides a semianalytical approximate formula.In Appendix A we have derived Eq. (A8), whose implementation leads to the red curve presented in Fig. 9(c).It provides a qualitative justification of the phase diagram: for instance, we see that the synchronous solution is unstable in the interval I 2 , where χ < 1.By following the approach developed in Ref. [19], we can compute also the maximal Lyapunov exponent λ : it is given by the maximal eigenvalue of a suitable random matrix.The resulting values correspond to the green curve.The changes of sign of λ coincide almost exactly with the border of the intervals where the synchronous state ceases to be observed.
What is left to be understood is the regime observed within the interval I 2 : it differs from the perfectly synchronous state, but it is nevertheless nearly-synchronous.While approaching the left border of I 2 , where the synchronous state becomes stable, the width of the phase distribution progressively shrinks.This is clearly seen in Fig. 10, where four instantaneous phase distributions are plotted for decreasing β values (from red to green curve).The transition scenario occurring at the other edge of the interval I 2 appears to be different and further studies would be required.However, a comparative analysis of different models suggest that this regime follows from a suitable combination of refractoriness and the shape of the PRC.As we suspect not to be very general, we do not investigate it in further detail.
Finally, we have considered broader pulses, to test the robustness of our findings.More precisely, now we assume the pulse-widths α −1 , β −1 to be longer than the refractory time t r as observed in real neurons 4,21 .The results are displayed in Fig. 11 for α = 12 and µ = 0.3.Once again, we see that CID extends to the region where β < α and that the transition point moves progressively towards β = α upon increasing the network size (see the different curves).On the other hand the strength of CID is significantly low (χ = 0.11), pos- sibly due to the relative smallness of the coupling strength.Furthermore, the evolution of quasi-synchronous solutions (δ p = 10 −3 ), reveals again bistability in a relatively wide interval of β -values, β 8.5 − 14.3, which now extends beyond β = α: a result, compatible with the transversal stability (see the red curve for λ c in Fig. 11).

A. Robustness
In the previous sections we have investigated the dependence of CID on the spike-width as well as on the coupling strength.Now, we examine the role of the PRC shape.Following Fig. 1, we consider a couple of smoothened versions of PRC 1 , defined in section II.The results obtained for a network of N = 10000 neurons are reported in Fig. 12.
All simulations have been performed for α = 100, while β has been again varied in the range [20, 120].In each panel, blue circles, orange stars and green diamonds have been obtained by setting µ = 0.3; they correspond to PRC 1,2,3 respectively.As a first general remark, the overall scenario is not strongly affected by the specific shape of the PRC.The mean firing rate is approximately the same in all cases, while the coefficient of variation is substantially higher for the sinusoidal (and more realistic) PRC 3 .Moreover, the order parameter for PRC 3 is remarkably close to that for PRC 1 (see panel c).
The most substantial difference concerns the transition from synchrony to CID, which occurs much earlier in PRC 2 .On the other hand, the χ-behavior of PRC 2 can be brought to a much closer agreement by increasing the coupling strength (the green asterisks in Fig. 12 refer to µ = 0.7).This observation raises the issue of quantifying the effective amplitude of the coupling: PRCs are introduced in Sec.II are all functions whose maximum value is equal to 1.This does not exclude that the effective amplitude may be significantly different, deviation that can be partially removed by adjusting the value of the coupling constant µ as shown in Fig. 12. Anyhow, these qualitative arguments need a more solid justification.In fact, in this last case (PRC 2 and µ = 0.7) C v is significantly larger (above 0.6 instead of below 0.2), consistently with the analysis carried out in Ref. [25], where it is All snapshots correspond to the black triangles in Fig. 9.
shown that a large coupling strength induces a bursting phenomena in LIF neurons.Finally, we investigate the presence of hysteresis in the case of PRC 3 .The results, obtained by setting all parameters as in the previous cases, are reported in Fig. 12 (see black triangles): they have been obtained by setting the initial spread of phases δ p = 10 −3 .Once again, there exists a wide parameter range where CID coexists with a stable synchronous regime.At variance with the previous case (see Fig. 9), the syn- chronous state is always stable over the range β ≤ 110.This is consistent with the variation of the conditional Lyapunov exponent, which does not exhibit an "instability island".As from Eq. (A8), λ c is the sum of two terms.In the case of PRC 3 , the second one is absent because the PRC amplitude is zero at the reset value Φ r = 0.

VII. CONCLUSION AND OPEN PROBLEMS
In this paper we have discussed the impact of finite pulsewidths on the dynamics of a weakly inhibitory neuronal network, mostly with reference to the sustainment and stability of the balanced regime.
In computational neuroscience, both exponential 28 and αpulses 29,30 are typically studied.The former ones are simpler to handle, as they require one variable per neuron per field type (inhibitory/excitatory); the latter ones, being continuous, are more realistic, but require twice as many variables.In this paper we have selected exponential pulses to minimize the additional computational complexity.We have prioritized the analysis of short pulses (about hundredth of the interspike interval) in order to single out deviations from δ -spikes.However tests performed for relatively longer spikes suggest that the general scenario is substantially confirmed for tentimes longer pulses (a value compatible with the time scales of AMPA receptors 26,31 ).The main changes observed when decreasing α down to 12 (starting from our reference 100) is the disappearance of the quasi-synchronous regime for a small degree of asymmetry: this happens around α ≈ 60 ∼ 70.
Besides pulsewidth, the asymmetry between excitatory and inhibitory spikes (a parameter which does not make sense in the case of δ -pulses) plays a crucial role in the preservation of the balance between excitation and inhibition.In fact, upon changing the ratio between excitatory and inhibitory pulsewidth different regimes may arise.The role of time scales is particularly evident in the synchronous regime, where the overall field is the superposition of two suitably weighted exponential shapes with opposite sign: depending on the time of observation, the effective field may change sign signalling a prevalence of either inhibition or excitation.
Altogether CID is robust when inhibitory pulses are shorter than excitatory ones (this is confirmed by the corresponding instability of the synchronous regime).More intriguing is the scenario observed in the opposite case, when CID and synchrony maybe simultaneously stable.A finite-size analysis performed by simulating increasingly large networks shows that the hysteretic region progressively shrinks, although it is still prominent -especially for weak coupling -for N = 80000, when the connectivity of our networks (K = 8000) is comparable to that of the mammalian brain.Anyhow, on a purely mathematical level, one can argue that the transition from CID to synchrony eventually occurs for identical widths.
Further studies are definitely required to reconstruct the general scenario, since the dynamics depends on several parameters.Here, we have explored in a preliminary way the role of the PRC shape: so long as it is almost of Type I, the overall scenario is robust.
Finally, the transition from CID to synchrony requires more indepth studies.A possible strategy consists in mimicking the background activity as a pseudo-stochastic process, thereby writing a suitable Fokker-Planck equation.However, at variance with the δ -spike case, here additional variables would be required to account for the dynamics of the inhibitory/excitatory fields.
the different velocity (frequency) exhibited at threshold and at the end of the refractory period.Notice that in the limit of short pulses, the field amplitude at time t m is effectively negligible, and one can thereby neglect the effect of the fields and assume Φ(t m ) = 1.
In mean-field models, the conditional Lyapunov exponent coincides with the exponent obtained by implementing a rigorous theory which takes into account mutual coupling.

FIG. 2 .
FIG.2.Characterization of the global network dynamics with interactions through δ -pulses.Mean firing rate ν, mean coefficient of variations C v , and order parameter χ are plotted vs. the coupling strength µ in panels (a), (b) and (c), respectively.Black triangles, red circles, green crosses, and blue diamonds correspond to N = 10000, 20000, 40000, and 80000, respectively, all obtained with PRC 1 .Orange stars and green squares correspond to N = 10000 and 40000 obtained with PRC 3 .The vertical dashed line represents the critical coupling µ c = 0.537.

FIG. 3 .
FIG. 3. CID properties for PRC 1 , µ = 0.95, α = β = 100, and N = 10000.Panel a): time series of the mean field Φ .Panel b): raster plot of spiking times t n for 100 oscillators out of N = 10000.Panel c): instantaneous probability distribution of the phases P(Φ) at two different time points t = 8363 (red), and t = 8374 (green).The probability distributions are normalized such that the area underneath is 1.

FIG. 5 .
FIG. 5. Characterization of the global network dynamics for nonidentical finite pulse-width, obtained with α = 100 and PRC 1 .Each column refers to different coupling strengths: = 0.3 (A), µ = 0.47 (B), and µ = 0.95 (C).Rows: mean firing rate ν, mean coefficient of variations C v , and order parameter χ versus β .Colours and symbols define network sizes N: 10000 (black triangles), 20000 (red crosses), 40000 (orange circles), and 80000 (blue stars).Each data point is based on a time series generated over 10000 time units and sampled every 1000 steps after the transient has sorted out.

FIG. 9 .FIG. 10 .
FIG.9.The emergence of a bistable regime for nonidentical finite pulse-widths and PRC 1 .The parameter set is the same as in Fig.5Awith N = 10000.Panels: a) mean firing rate ν, b) mean coefficient of variations C v , and c) order parameter χ versus β .The blue circles and black triangles in all panels correspond to different initial conditions: fully random (circles), restricted to a tiny interval (triangles).The narrow ICs are chosen to be in the order of δ p = 10 −3 .The green diamonds corresponds to the maximal Lyapunov exponent, and the red one is the conditional Lyapunov exponent as function of β .The magenta line (a) represents the semi analytic firing rate given in Eq.A2.The horizontal dashed line (c) is a reference point (λ = 0) in which the synchronous state changes its stability.

FIG. 12 .
FIG.12.The robustness for other PRCs.The mean firing rate ν, the mean coefficient of variations C v , and the order parameter χ vs. inhibitory pulse widths β are shown in panel a), b) and c), respectively for N = 10000 and α = 100.PRC 1 with random ICs is shown for µ = 0.3 (blue circles) as reference to Figs.5(A) and 9.The PRC 2 with random ICs is depicted for µ = 0.3 (orange stars) and µ = 0.7 (green stars).Green diamond and black triangles result from PRC 3 and µ = 0.3.The former has been created with random ICs and the latter with strongly restricted ICs with δ p = 10 −3 within the narrow basin of attraction for the synchronous attractor.The magenta curve (panel a) represents the semi-analytic firing rate for PRC 3 according to Eq. A2.Panel c) shows on the alternative y-axis also the conditional Lyapunov exponent λ c (red curve) for synchronous solutions and PRC 3 .The horizontal red dashed line is the null line to the axis on the right.