Collective dynamics in the presence of finite-width pulses
TThe effect of pulse width on the dynamics of pulse-coupled oscillators
Afifurrahman, a) Ekkehard Ullner, b) and Antonio Politi c) Institute for Pure and Applied Mathematics and Department of Physics (SUPA), Old Aberdeen, Aberdeen, AB24 3UE,United Kingdom (Dated: 9 February 2021)
The idealisation of neuronal pulses as δ -spikes is a convenient approach in neuroscience but can sometimes lead toerroneous 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 phaseoscillators coupled through exponential pulses with different widths. We consider three coupling functions, inspiredby 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 fullysynchronous regime if the inhibitory pulses are sufficiently wider than the excitatory ones. The transition to synchronyis 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. Aconjectured first-order phase transition emerging for δ -spikes is smoothed out for finite-width pulses. Neuronal networks with a nearly balanced excita-tory/inhibitory activity evoke significant interest in neu-roscience due to the resulting emergence of strong fluc-tuations akin to those observed in the resting state of themammalian brain. While most studies are limited to a δ -like pulse setup, much less is known about the collectivebehavior in the presence of finite pulse-widths. In this pa-per, we investigate exponential pulses, with the goal of test-ing the robustness of previously identified regimes such asthe spontaneous emergence of collective irregular dynam-ics (CID). Moreover, the finite-width assumption paves theway to the investigation of a new ingredient, present inreal neuronal networks: the asymmetry between excita-tory and inhibitory pulses. Our numerical studies confirmthe emergence of CID also in the presence of finite pulse-width, although with a couple of warnings: (i) the ampli-tude of the collective fluctuations decreases significantlywhen the pulse-width is comparable to the interspike in-terval; (ii) CID collapses onto a fully synchronous regimewhen the inhibitory pulses are sufficient longer than theexcitatory ones. Both restrictions are compatible with therecorded behavior of real neurons. Additionally, we findthat a seemingly first-order phase transition to a (quasi)-synchronous regime disappears in the presence of a finitewidth, confirming the peculiarity of the δ -spikes. A tran-sition to synchrony is instead observed upon increasingthe ratio between the width of inhibitory and excitatorypulses: this transition is accompanied by a hysteretic re-gion, which shrinks upon increasing the network size. In-terestingly, for a connectivity comparable to that of themammalian brain, such a finite-size effect is still sizable.Our numerical studies might help to understand abnor-mal synchronisation in neurological disorders. a) Electronic mail: [email protected] b) Electronic mail: [email protected] c) Electronic mail: [email protected]
I. INTRODUCTION
Irregular firing activity is a robust phenomenon observedin certain areas of mammalian brain such as hippocampus orcortical neurons . It plays a key role for the brain functioningin the visual and prefrontal cortex. This behavior does notonly emerge from intrinsic properties of the single neuronalcells, but it can also arise from the combined action of manyinteracting units .Several numerical and theoretical studies of minimal mod-els have been performed to uncover the underlying mecha-nisms. In binary networks, for instance, the irregular activ-ity was hypothesized to be the result of a (statistical) bal-ance between excitation and inhibition . The resulting bal-anced state is characterized by a constant (on average) firingrate accompanied by irregular fluctuations, which can treatedas stochastic-like. In fact, this idea allowed Brunel to de-velop a powerful analytical method to describe a more real-istic model of leaky integrate-and-fire (LIF) neurons. Alto-gether this regime can be seen as a fluctuating asynchronousstate, with negligible correlations among the various neurons.A numerical study of an ensemble of heterogeneous LIFneurons with delay has independently revealed that irregularactivity can also emerge in a strictly-inhibitory environment .This self-sustained regime, that we call “collective irregulardynamics" (CID), differs from the asynchronous regime inthat the temporal fluctuations persist in the so-called thermo-dynamic limit, when the number N of neurons is let diverge toinfinity . In other words, we are in the presence of a collec-tive phenomenon, akin to partially synchronized states . Therobustness of the mechanism has been confirmed by furtherstudies of pulse-coupled phase oscillators in the absence ofdelay.The stochastic-like dynamics in CID can be of differenttypes: (i) the manifestation of low-dimensional chaos in thepresence of two distinct populations ; (ii) high-dimensionalchaos for a single population of oscillators ; (iii) stable chaosfound in a diluted network of LIF and in more disorderedsetups . Stable chaos is a form of irregular dynamics,which emerges in the presence of negative Lyapunov expo- a r X i v : . [ q - b i o . N C ] F e b nents akin to the behavior of chaotic cellular automata.Recently, it has been shown that CID arises also in the typ-ical balanced setup, under the condition of small unbalanceand in the limit of massive coupling (finite connectivity den-sity) . It is therefore of utmost importance to test the ro-bustness of such a phenomenon, by exploring more realisticsetups. In this paper we explore the role played by the spike-width. Most of the numerical studies found in the literaturedeal with δ -spikes. This is an idealization, typically promptedby the easiness of the related simulations, but it has its ownlimits. First, real pulses have a small but finite width . Sec-ond, it has been shown that the stability of some synchronousregimes of LIF neurons may qualitatively change when ar-bitrarily short pulses are considered (in the thermodynamiclimit) . For this reason, it is wise to extend the analysis ofCID to networks of neurons interacting via finite-width pulses.A preliminary study has been already published inRef. [20], where the authors have not performed any finite-size scaling analysis and, more important, no any test of thepresence of CID has been carried out. Here we study a sys-tem composed of two populations of (identical) excitatory andinhibitory neurons, which interact via exponential pulses ofdifferent width, as it happens in real neurons .Handling pulses with a finite width requires two additionalvariables per single neuron, in order to describe the evolu-tion of the incoming excitatory and inhibitory fields. The cor-responding mathematical setup has been recently studied inRef. [22] with the goal of determining the stability of the fullysynchronous state in a sparse network. The presence of twodifferent pulse-widths leads to non-intuitive stability proper-ties, because the different time dependence of the two pulsesmay change the excitatory/inhibitory character of the overallfield perceived by each single neuron. Here, we basically fol-low the same setup introduced in Ref. [22] with the main dif-ference of a massively coupled network, to be able to performa comparative analysis of CID.The randomness of the network accompanied by the pres-ence of three variables per neuron, makes an analytical treat-ment quite challenging. For this reason we limit ourselves toa numerical analysis. However, we accompany our studieswith a careful finite-size scaling to extrapolate the behavior ofmore realistic (larger) networks. Our first result is that CID isobserved also in the presence of finite pulse-width, althoughwe also find a transition to full synchrony when the inhibitorypulses are sufficiently longer than excitatory ones. The transi-tion is first-order (discontinuous) and is accompanied by hys-teresis: there exists a finite range of pulse-widths where CIDand synchrony coexist.The finite-size analysis suggests that in the thermodynamiclimit CID is not stable when the pulse-width of inhibitory neu-rons is strictly longer than that of the excitatory ones. How-ever, the convergence is rather slow and we cannot excludethat the asymmetry plays an important role in real neuronalnetworks of finite size.The remainder of this paper is organized as follows. Insection II, we define the model, including the phase responsecurves (PRCs) used in our numerical simulation. In the samesection we also introduce the tools and indicators later used to characterize the dynamical regimes, notably an order parame-ter to quantify the degree of synchronization . In section III,we present some results obtained for strictly δ pulses to testrobustness of CID in our context of coupled phase oscillators.In Sec. IV we discuss the symmetric cases of identical finitepulse-widths. Sec. V is devoted to a thorough analysis of CIDby varying the finite pulse-widths. Sec. VI contains a discus-sion of the transition region, intermediate between standardCID and full synchrony. In the same section, the robustnessof the transition region is analysed, by considering differentPRCs. Finally, section VII is devoted to the conclusions and abrief survey of the open problems. II. MODEL
Our object of study is a network of N phase oscillators (alsoreferred to as neurons), the first N e = bN being excitatory,the last N i = ( − b ) N inhibitory (obviously, N e + N i = N ).Each neuron is characterized by the phase-like variable Φ j ≤ G with entries G j , k = (cid:40) , if k → j active0 , otherwisewhere ∑ N e k = G j , k = K e and ∑ Nk = N e + G j , k = K i , meaning thateach neuron j is characterized by the same number of incom-ing excitatory and inhibitory connections, as customary as-sumed in the literature ( K = K e + K i represents the connec-tivity altogether). Here, we assume that K is proportional to N , that is K = cN , i.e. we refer to massive connectivity. Fur-ther, the network structure is without autapse, i.e. G j , j = Φ j = + µ √ K Γ (cid:0) Φ j (cid:1) (cid:0) E j − I j (cid:1) , (1)where E j ( I j ) the excitatory (inhibitory) field perceived by the j th neuron, while Γ ( Φ ) represents the phase-response curve(PRC) assumed equal for all neurons; finally, µ is the cou-pling strength. Whenever Φ k reaches the threshold Φ th = Φ r = t r during which it stands still and is insensitive to the actionof both fields. The fields E j and I j are the linear superposi-tion of exponential spikes emitted by the upstream neurons.Mathematically,˙ E j = − α (cid:18) E j − ∑ n G j , k P k δ ( t − t kn ) (cid:19) (2)˙ I j = − β (cid:18) I j − g ∑ n G j , k ( − P k ) δ ( t − t kn ) (cid:19) , where α ( β ) denotes the inverse pulse-width of the excitatory(inhibitory) spikes and t kn is the emission time of the n th spikeemitted by the k th neuron. The coefficient g accounts for the Φ Φ Γ(Φ) Φ FIG. 1. Example of the phase response curves (PRCs): PRC with Φ = − . Φ = . (red dashed line), PRC (blue dashed and dot line). The vertical dot line refers to the resetmembrane potential ( Φ r = relative strength of inhibition compared to excitation. P k = k th neuron is excitatory, otherwise P k = α ( β ) → ∞ ( δ -spikes) both fields can be ex-pressed as simple sums E j = ∑ n G j , k P , k δ ( t − t kn ) (3) I j = g ∑ n G j , k ( − P k ) δ ( t − t kn ) . Let us now introduce the PRCs used later in our numericalsimulations. We consider three different shapes:• PRC Γ ( Φ j ) = (cid:40)(cid:0) Φ j − Φ (cid:1) if Φ < Φ j < Φ Γ ( Φ j ) = Φ j − Φ . − Φ if Φ < Φ j < . − (cid:16) Φ j − . Φ − . (cid:17) if 0 . < Φ j < Φ Γ ( Φ j ) = sin (cid:0) π Φ j (cid:1) (6)The various curves are plotted in Fig. 1. PRC (see theblack curve, which corresponds to Φ = − . Φ = . (see Eq. (6)), as aprototype of type I PRC .The network dynamics is simulated by implementing theEuler algorithm with a time step δ t = − . However, in somecases, smaller steps have been considered to avoid problemsof spurious synchronization. We typically initialize the phasesuniformly in the unit interval, while the fields are initially setequal to zero.In all cases, we have set b = . c = . g = + (cid:112) / K (following the existing literature). The last optionis justified by the motivation to maintain a balanced regime for K , N → ∞ . We have, instead, systematically explored the roleof α and β , as the pulse-width is the focal point of this paper.Additionally, the coupling strength µ has been varied, as wellas the network-size N , to test for the amplitude of finite-sizeeffects.The following statistical quantities are used to characterizethe emerging dynamical states.1. The mean firing rate is a widely used indicator to quan-tify the neural activity. It is defined as ν = lim t → ∞ tN N ∑ j = N j ( t ) (7)where N j ( t ) denotes the number of spikes emitted bythe neuron j over a time t .2. The coefficient of variations C v is a microscopic mea-sure of irregularity of the dynamics based on the fluctu-ations of the interspike intervals (ISIs). The average C v is defined as C v = N N ∑ j = σ j τ j , (8)where σ j is the standard deviation of the single-oscillator’s ISI, and τ j is the corresponding mean ISI. If C v >
1, then the neurons show a bursting activity, while C v < The order parameter , χ , is typically used to quantify thedegree of synchronisation of a population of neurons .It is defined as χ ≡ (cid:104) Φ (cid:105) − (cid:104) Φ (cid:105) (cid:104) Φ − Φ (cid:105) , (9)where (cid:104)·(cid:105) represents an ensemble average, while over-bar is a time average. The numerator is the variance ofthe ensemble average (cid:104) Φ (cid:105) , while the denominator is theensemble mean of the single-neuron’s variances. Whenall neurons behave in exactly the same way (perfectsynchronization), then χ =
1. If instead, they are in-dependent, then χ ≈ / √ N . Regimes characterized byintermediate finite values 0 < χ < χ -value strictly larger than zero;here we are interested in the so-called collective irreg-ular dynamics (CID), characterized by stochastic-likefluctuations of the average (cid:104) Φ (cid:105) . µ C v µ χ µ ν a) b)c) FIG. 2. Characterization of the global network dynamics with inter-actions through δ pulse, PRC and respectively PRC . The meanfiring rate ν , mean coefficient of variations C v , and order parameter χ vs. coupling strength µ are shown in the panels (a), (b) and (c),respectively. The network size N is encoded in the symbols and colorfor PRC : 10000 (black triangles), 20000 (red circles), 40000 (greencrosses), and 80000 (blue diamonds). Similarly for PRC : 10000(orange stars) and 40000 (green squares). The vertical dashed linerepresents a critical coupling µ c = . III. DELTA PULSE
Most spiking network models deal with δ -spikes, includingthose giving rise to CID . This paper is mostly focused onthe more realistic exponential spikes, but before proceedingin that direction we wish to briefly discuss the case of zeropulse-width. This is useful to confirm the ubiquity of CIDwithin the class of pulse-coupled oscillators and as referenceto identify possibly new phenomena induced by the presenceof finite widths. δ pulses correspond to the limiting case α , β → ∞ ; they canbe treated by invoking Eq. (3). Figure 2 shows the various in-dicators introduced in Section II, to characterize the collectivedynamics. As in previous papers, we explore the parame-ter space, by varying the coupling strength µ and the systemsize N .In panel (c) we can appreciate that CID emerges already forvery small coupling strength; it is accompanied by an increas-ing coefficient of variations C v , due to the coupling which in-duces increasing deviations from purely periodic behavior. Inparallel, the mean firing rate ν decreases as a result of theprevalent inhibitory character of the network. This weak-coupling emergence of CID is comparable to what observedin balanced LIF models with δ spikes .Above µ ≈ . = µ c (see the vertical dashed lines), atransition occurs towards a highly synchronous regime ( χ isslightly smaller than 1), accompanied by a larger firing rate.The corresponding firing activity is mildly irregular: C v issmaller than in Poisson processes (when C v = Φ th =
1. In fact, similar studies performedwith PRC do not reveal any evidence of a phase transition(see orange stars and green squares in Fig. 2) indicating thatsuch behavior is nothing else but a peculiarity of PRC with δ -pulses. We have not further explored this regime. It is never-theless worth noting that the sudden increase of the firing rateobserved when passing to the strong coupling regime is rem-iniscent of the growth observed in LIF neurons , althoughin such a case, the increase is accompanied by a significantlybursty behavior .More important is the outcome of the finite-size scalinganalysis, performed to investigate the robustness of the ob-served scenario. In Fig. 2 one can see that the various indica-tors under stimulation of PRC are size-independent deeplywithin the two dynamical phases, while appreciable devia-tions are observed in the transition region. This is customarywhen dealing with phase-transitions. It is not easy to con-clude whether the transition is either first or second order: the C v is reminiscent of the divergence of susceptibility seen incontinuous transitions, but this is an issue that would requireadditional work to be assessed. IV. IDENTICAL FINITE-WIDTH PULSES
In this section, we start our analysis of finite pulses, by as-suming the same width for inhibitory and excitatory neurons,i.e. α − = β − . The asymmetric case is discussed in the nextsection. All other system parameters are kept the same as inthe previous section (including the PRC shape).The resulting phase diagram is reported in Fig. 3. In panel(a) we plot the order parameter χ as a function of µ for differ-ent widths: from broad pulses (red stars correspond to α =
1, awidth comparable to the ISI), down to very short ones (greentriangles correspond to α = δ -pulses (incorrespondence of the vertical dashed line at µ c = . δ -spikes is confirmed by the relativelylarge deviations appearing already for α = α is presented inFig. 3(b), where we plot χ versus α for different couplingstrengths: µ = . N = N = µ χ α χ a)b) FIG. 3. Characterization of the global network dynamics in the pres-ence of identical pulse-width and PRC . Panel a): order parameter χ vs. µ for N = α = α = α =
10 (orange squares), and α = δ pulses (see Fig. 2 (c), N = µ c derived therein. Panel b): order parameter χ vs. α for N = µ = . µ = .
47 (red pluses), and µ = . α -axis is in logarithmic scale. Inboth panels the magenta circles show the results for N = V. FULL SETUP
In the previous section we have seen that the finite widthof the spikes does not kill the spontaneous emergence of CID.Here, we analyse the role of an additional element: the asym-metry between inhibitory and excitatory pulses. We proceedby exploring the two-dimensional parameter space spannedby the coupling strength µ and the asymmetry between pulsewidths. The latter parameter dependence is explored by set-ting α =
100 and letting β (the inverse width of inhibitorypulses) vary. All other network parameters, including the PRCshape, are assumed to be the same as in the previous section.We start by illustrating the collective dynamics in a systemof N = µ = .
95. Fig. 4refers to a weakly asymmetric case, β = t n of a subset of 100 neu-rons: there, one can easily spot the time intervals character-ized by a more coordinated action (see, for instance, aroundtime 80 ∼
90 in Fig. 4(a)). A more quantitative representationis presented in Fig. 4(c), where the instantaneous phase dis-tribution P ( Φ ) is plotted at three different times in correspon- dence of three qualitatively different regimes of the phase dy-namics (see the vertical lines in panel (a)). The peak at Φ = P ( Φ ) presented in Fig. 4(c), show that, atvariance with the classical asynchronous regime, the shape ofthe probability density changes with time. The narrowest dis-tribution (blue curve) corresponds to the region where strongregular oscillations of (cid:104) Φ (cid:105) are visible in panel (a): within thistime interval a “cloud" of neurons homogeneously oscillatesfrom reset to threshold and back.The results of a systematic numerical analysis are plottedin Fig. 5, where we report three indicators: the firing rate ν ,the mean coefficient of variation C v , and the order parame-ter χ , versus β for three different coupling strengths (see thedifferent columns), and four network sizes.All indicators reveal the existence of two distinct phases: asynchronous regime arising for small β values, and CID ob-served beyond a critical point which depends on the networksize: the transition is discontinuous. All panels reveal a sub-stantial independence of the network size, with the exceptionof the transition between them (we further comment on thisissue later in this section).The first regime is synchronous and periodic, as signalledby χ =
1, and C v =
0. The corresponding firing rate ν is a bitsmaller than 0.97, the value expected in the limit of uncoupledneurons (taking into account refractoriness). This is consistentwith the expected predominance of inhibition over excitationin this type of setup. A closer look shows that in the syn-chronous regime ν increases with β . This makes sense sincethe smaller β , the longer the time when inhibition prevailsthereby decreasing the network spiking activity. The weak de-pendence of ν on the coupling strength µ is a consequence ofsmall effective fields felt by neurons when the PRC is small.Finally, for intermediate β values (around 80) and large cou-pling strengths, χ is large but clearly smaller than 1. This thirdtype of regime will be discussed in the next section.CID is characterized by a significantly smaller order pa-rameter which, generally tends to increase with the couplingstrength. CID is also characterized by a significantly smallerfiring rate. This is due the prevalence of inhibition which isnot diminished by the refractoriness as in the synchronousregime. Finally, the coefficient of variation is strictly largerthan 0, but significantly smaller than 1 (the value of Poissonprocesses) revealing a limited irregularity of the microscopicdynamics. In agreement with our previous observations for δ -spikes, C v increases with the coupling strength.Our observation of an irregular regime for β < α are in linewith the numerical observations discussed in Ref. [20], withreference to LIF neurons and with the hypothesis that a slowerinhibitory activity might lead to rhythmic oscillations . Onthe other hand, our finite-size scaling analysis also shows thatthe degree of asymmetry (between pulse widths) compatiblewith CID progressively reduces upon increasing the numberof neurons. Although the N -dependence varies significantlywith the coupling strength, it is natural to conjecture that, t < Φ > Φ t n j a) b) c) P( Φ ) FIG. 4. Microscopic phase dynamics for nonidentical pulse-width. All presented data refers to PRC , µ = . α / β = /
95 and N = t n for 100 oscillators out of N = ( Φ ) at three different time points t =
10 (green), t =
40 (red), and t =
85 (blue) time units. The curves are normalized such that the area underneath is 1.
20 40 60 80 100 12000.20.40.60.81 χ
20 40 60 80 100 120 20 40 60 80 100 120 14000.20.40.60.8 C v ν β A B C
FIG. 5. Characterization of the global network dynamics for nonidentical pulse-width and PRC . Each column refers to different couplingstrengths: µ = . µ = .
47 (B), and µ = .
95 (C). Rows: mean firing rate ν , mean coefficient of variations C v , and order parameter χ versus β . Colours and symbols in all panels define various network sizes N : 10000 (black triangles), 20000 (red crosses), 40000 (orangecircles), and 80000 (blue stars). Each data point is based on a time series generated over 10000 time units and sampled every 1000 steps afterthe transient has sorted out. S(f) β =60 β =90 β =120 S(f) β =75 β =85 β =110 f S(f) β =92 β =95 β =120 a)b)c) [a.u.][a.u.][a.u.] FIG. 6. Power spectra of average membrane potential as functionof frequency. All presented data refers to PRC , α =
100 and N = µ : 0.3 (a), 0.47 (b), and0.95 (c). asymptotically, CID survives only for β ≥ α . This is not toosurprising from the point of view of self-sustained balancedstates. They are expected to survive only when inhibition andexcitation compensate each other: the presence of differenttime scales makes it difficult, if not impossible to ensure asteady balance.We conclude this section with a more quantitative charac-terization of the irregularity of the collective dynamics. InFig. 6, we plot the Fourier power spectrum S ( f ) obtainedfrom (cid:104) Φ (cid:105) ( t ) . The panels correspond to three different cou-pling strengths ( µ = .
3, 0.47 and 0.95, from top to bottom).For each value of µ , we have sampled three different pulse-widths.Altogether, one can notice a general increase of the powerwith µ . This is quite intuitive, as CID is the result of mutualinteractions. A less obvious phenomenon is the increase ofthe power observed when the inhibitory pulse-width β − isincreased. This is an early signature of a transition towardsfull synchronization, observed when β is decreased below acritical value. Anyway, the most important message conveyedby Fig. 6 is that all spectra exhibit a broadband structure, al-though most of the power is concentrated around a specificfrequency: f ≈ . f ≈ . f ≈ . S(f)
S(f) f S(f) a) b) c) [a.u.][a.u.][a.u.] FIG. 7. Power spectra of average membrane potential as function offrequency. Panel: a) β = µ = .
3, b) β = µ = . β = µ = .
95. Each colour in all panels refers to different networksize N : 20000 (red), 40000 (green), and 80000 (blue). All presenteddata refers to PRC and α = ν ≈ .
523 for µ = . ν ≈ .
44 for µ = . rections, by computing the power spectrum S ( f ) for differ-ent network sizes for three different parameter sets: µ = . β =
90 (panel a), µ = . β =
120 (panel b), and µ = . β =
95 (panel c). In all cases, the spectra are substantiallyindependent of the number of neurons, although in panel (b)we observe a weak decrease in the band f ∈ [ , . ] , while anew set of peaks is born in panel (c). Since the connectiv-ity K of the largest networks herein considered ( N = , K = ), we can at least conjecture that this phenomenon mayhave 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 theconnection with the microscopic firing rate is lost in panel (c). VI. TRANSITION REGION
In the previous section we have seen a clear evidence ofa first-order phase transition, when either the pulse-width orthe coupling strength is varied. So far, each simulation hasbeen performed by selecting afresh a network structure. Thestability of our results indicates that the transition does notsuffer appreciable sample-to-sample fluctuations.The main outcome of our numerical simulations is sum-marized in Fig. 8; the various lines identify the transition be-tween the two regimes, for different network sizes. The crit-ical points have been determined by progressively decreasing
40 50 60 70 80 90 100 β µ SYNCHRONY CID
FIG. 8. Phase diagram for different network sizes N : 10000 (blacktriangles), 20000 (red squares), 40000 (orange circles), and 80000(blue stars). β (see Fig. 5) and thereby determining the minimum β -valuewhere CID is stable. Upon increasing N , the synchronizationregion grows and the transition moves towards β = α .So far, the initial condition has been fixed by selecting inde-pendent, identically uniformly distributed random phases andzero fields. Since it is known that discontinuous transitionsare often accompanied by hysteretic phenomena, we now ex-plore this question. Fig. 9 combines results presented in theprevious section for a network with N = µ = . δ p (while the fields are setequal to zero and δ t = − ). The resulting values of thevarious observables for δ p = − are reported in Fig. 9: seethe black curves and triangles. The most relevant result is thediscovery of a wide range of β -values, where the order pa-rameter is either exactly equal to 1 or close to it. This impliesthe presence of a bistability between CID and a stable quasi-synchronous solution with a small basin of attraction.Linear-stability analysis of the periodic synchronous stateprovides some relevant informations. The conditional Lya-punov exponent λ c can be determined in a semi-analyticform: its expression for PRC is derived in Appendix (A):see Eq. (A8). The results are presented in Fig. 9(c) (see thered curve): there we see that when the synchronous solutionis transversally unstable, the quasi-synchronous regime sys-tematically arises. However, the reverse is not true: there are β values where the synchronous state is stable and yet CID isobserved. The inconsistency is partially solved by looking ata more accurate quantifier of the stability of the periodic so-lution: the exact maximal Lyapunov exponent λ . By follow-ing the approach developed in Ref. [22], λ can be determinedby computing the eigenvalues of a suitable random matrix.The implementation of this method leads to the green curve,slightly above the conditional Lyapunov exponent. The new curve partially closes the gap between the stability limit of theperiodic state and the onset of CID, but the general scenario isstill unclear.A careful look allows identifying the dynamical regimes ac-cording to different β intervals. The first one comes from theobservation of random ICs: [ , ] is characterized by stableCID (see blue circles in Fig. 9(c)). Below β ≈
42 the irreg-ular dynamics collapses into periodic states. As confirmedby the finite-width perturbation analysis (black triangles inFig. 9(c)), the periodic dynamics is stable over a range [ , ] followed by a hysteretic behaviour in [ , ] . In the lat-ter interval, depending on the initial condition, one can eitherobserve CID or periodic states. Understanding the nature ofthe transition from CID to synchrony would require additionalwork.More specifically, [ , ) is characterized by the fully syn-chronous dynamics. This is consistent with the linear stabil-ity analysis, which also shows that β ≈
46 is the value be-yond which the synchronous state becomes unstable (see the λ curve).Symmetrically, for β ∈ [ , ] , we observe the simulta-neous presence of CID and full synchrony in a region wherethe synchronous-state stability is restored (see again λ ). Thebasin of attraction for the fully and quasi-synchronous solu-tion is quite small in the range β ∈ [ , ] and relies ona sufficient small spread in the initial conditions δ p . For δ p ≥ . β ∈ [ , ] .In that range, the fully synchronous solution is linearly unsta-ble, but yet the phases remain confined to a small region andquasi-synchrony is stable. The phenomenon is a combined ef-fect of refractoriness and the shape of the PRC. As we suspectit to be rather peculiar, we do not investigate in detail but limitourselves to illustrate it in Fig. 10, where we report four snap-shots of the phase distributions close to the β value, wherethis regime ceases to exist. The width of the various distri-butions progressively shrinks while approaching the criticalpoint, where the synchronous state becomes stable: the greencurve, which corresponds to β = .
3, is almost δ -like. Themechanism responsible for the loss of stability of this quasi-synchronous regime at the other edge of the interval [ , ] appears to be different and further studies would be required.Finally, we have performed numerical studies with broaderpulses, to test the robustness of our findings. More precisely,now we assume the pulse-widths α − , β − to be longer thanthe refractory time t r as observed in real neurons . Theresults are displayed in Fig. 11 for α =
12 and µ = . β < α and that the transition point moves progressively to-wards β = α upon increasing the network size (see the dif-ferent curves). On the other hand the strength of CID is sig-nificantly low ( χ = . δ p = − ), reveals again bista-bility in a relatively wide interval of β -values, β (cid:39) . − . β = α : a result, compatible withthe transversal stability (see the red curve for λ c in Fig. 11).
20 40 60 80 120 β C v
20 40 60 80 100 120 β χ
20 40 60 80 120 β ν λ λ c a) b) c) FIG. 9. The emergence of a bistable regime for nonidentical pulse widths and PRC . The parameter set is the same as in Fig. 5A with N = ν , b) mean coefficient of variations C v , and c) order parameter χ versus β . The blue circles and blacktriangles in all panels correspond to the full network with a full range (uniform) and narrow ICs, respectively. The narrow ICs are chosen to bein the order of δ p = − . The green diamonds corresponds to the maximal Lyapunov exponent, and the red one is the conditional Lyapunovexponent as function of β . The magenta line (a) represents the semi analytic firing rate given in Eq.A2. The horizontal dashed line (c) is areference point ( λ =
0) in which the synchronous state changes its stability. Φ P( Φ ) FIG. 10. The instantaneous probability distribution P( Φ ) at particu-lar time point. Colour and shape define different β values orderedfrom the widest to the narrowest distribution: 70 (red), 50 (black),48 (blue), and 46.3 (green). All snapshots are corresponding to theblack triangles in Fig. 9. The vertical axis is in the logarithmic scale.The curves are normalized such that the area underneath is 1. A. Robustness
In the previous sections we have investigated the depen-dence of CID on the spike-width as well as on the coupling β χ λ c FIG. 11. Characterization of the global network dynamics for longpulse width and PRC . The parameters µ = . α =
12 arefixed. The blue circles ( N = N = N = N = N = δ p = − . The red curve is the con-ditional Lyapunov exponent. strength. Now, we examine the role of the PRC shape. Fol-lowing Fig. 1, we consider a couple of smoothened versions ofPRC , defined in section II. The results obtained for a network0of N = α = β has been again varied in the range [ , ] . In each panel,blue circles, orange stars and green diamonds have been ob-tained by setting µ = .
3; they correspond to PRC respec-tively. As a first general remark, the overall scenario is notstrongly affected by the specific shape of the PRC. The meanfiring rate is approximately the same in all cases, while the co-efficient of variation is substantially higher for the sinusoidal(and more realistic) PRC . Moreover, the order parameter forPRC is remarkably close to that for PRC (see panel c).The most substantial difference concerns the transitionfrom synchrony to CID, which occurs much earlier in PRC .On the other hand, the χ -behavior of PRC can be brought toa much closer agreement by increasing the coupling strength(the green asterisks in Fig. 12 refer to µ = . µ as shown in Fig. 12.Anyhow, these qualitative arguments need a more solid jus-tification. In fact, one can also notice that in this last case(PRC and µ = . C v is significantly larger (above 0.6 in-stead of below 0.2), consistently with the analysis carried outin Ref. [28], where it is shown that a large coupling strentghinduces a bursting phenomena in LIF neurons.Finally, we investigate the presence of hysteresis in the caseof PRC . The results, obtained by setting all parameters as inthe previous cases, are reported in Fig. 12 (see black trian-gles): they have been obtained by setting the initial spreadof phases δ p = − . Once again, there exists a wide pa-rameter range where CID coexists with a stable synchronousregime. At variance with the previous case (see Fig. 9), thesynchronous state is always stable over the range β ≤ , does not ex-hibit an “instability island". As from Eq. (A8), λ c is the sumof two terms. In the case of PRC , the second one is absentbecause the PRC amplitude is zero at the reset value Φ r = VII. CONCLUSION AND OPEN PROBLEMS
In this paper we have discussed the impact of finite pulse-widths on the dynamics of a weakly inhibitory neuronal net-work, especially with reference to the sustainment and stabil-ity of the balanced regime. We have considered the simplestcase, namely exponential pulses, whose simulation requiresjust one variable per field (excitatory/inhibitory) per neuron.We have found that the collective irregular dynamics persists,at least in so far as the width of the spikes is small comparedto the average interspike intervals (as it happens in reality).The stability analysis of the synchronous regime, togetherwith direct numerical simulations, reveals a crucial roleplayed by the asymmetry between excitatory and inhibitoryspikes. Indeed, the case of identical pulse-widths is very pe- culiar (degenerate), as the balance between inhibition and ex-citation is automatically maintained over different time scales.On the other hand, different pulse-widths imply that either ex-citation or inhibition may prevail over different time ranges,potentially breaking the balance. This is particularly evidentin the synchronous regime where the overall field is the su-perposition of two suitably weighted exponential shapes withopposite sign: depending on the time, the effective field maychange sign signalling a prevalence of either inhibition or ex-citation.Altogether CID is robust when inhibitory pulses are shorterthan excitatory ones (this is confirmed by the correspondinginstability of the synchronous regime). More intriguing is thescenario observed in the opposite case, when CID and syn-chrony maybe simultaneously stable. A finite-size analysisperformed by simulating increasingly large networks showsthat the hysteretic region progressively shrinks, although itis still prominent - especially for weak coupling - for N = , K = , δ -spike case, here additional variableswould be required to account for the inhibitory/excitatoryfields. ACKNOWLEDGMENTS
Afifurrahman was supported by the Ministry of Financeof the Republic of Indonesia through the Indonesia En-dowment Fund for Education (LPDP) (grant number: PRJ-2823/LPDP/2015).
DATA AVAILABILITY
The data that support the findings of this study are availablefrom the corresponding author upon reasonable request.
Appendix A: Mean field model for finite-width pulse
We investigate the stability of the period-1 synchronousstate through the conditional Lyapunov exponent. This regimeis characterised by a synchronous threshold-passing of all os-cillators leading to exactly the same exponentially decayingexcitatory and inhibitory field for all oscillators. The syn-chronous solution Φ ( t ) with a period T of Eqs. (1,2) is ob-1
20 40 60 80 120 β C v
20 40 60 80 100 120 β χ
20 40 60 80 120 β ν
60 80 100 120-0.500.51 -0.0200.020.04 a) b) c) χ λ c β FIG. 12. The emerging scenario for different PRCs with µ = . µ = . ), α = , N = ν , b) mean coefficient of variations C v , and c) order parameter χ vs. β . Blue circles correspond to PRC , orange starsPRC , green stars PRC with µ = .
7, green diamond PRC . Inset: a bistable regime emerges for PRC in a network of N = δ p = − (black triangels) accompanied by the conditional Lyapunov exponent (red curve). The magenta curve (panel a) represents thesemi-analytic firing rate for PRC according to Eq.A2. tained by integrating the equation ˙ Φ = + µ √ K Γ ( Φ )( E ( t ) − I ( t )) E ( t ) = E ◦ e − α t I ( t ) = I ◦ e − β t (A1)where E ◦ = K e α − e − α T , I ◦ = gK i β − e − β T . The fields follow an exponential decay with the initial am-plitudes E ◦ , I ◦ for the excitatory and inhibitory field, respec-tively. In order to determine the stability of the synchronousstate, we first need to find the period T via a self-consistentiterative approach. Setting the origin t = T as a timepoint in which the phase variable reaches its maximal valuei.e. Φ ( T ) =
1. We integrate the phase starting from Φ = Φ ( T ) by initially imposing arbitrary non-zero values for E ◦ and I ◦ . The period T is used to restart the procedure withrefined estimate of the initial field amplitudes E ◦ , I ◦ , until con-vergence to a fixed point is ensured. The firing rate is givenby, ˜ ν ≡ T (A2)The conditional (also known as transversal) Lyapunov ex-ponent is a simple tool to assess the stability of the syn-chronous regime. It quantifies the stability of a single neuron subject to the external periodic modulation resulting from thenetwork activity. The transversal, or conditional Lyapunovexponent is the growth rate λ c of an infinitesimal perturba-tion. Let us now denote with δ t r the time shift at the end of arefractory period. The corresponding phase shift is δ φ r = ˙ Φ ( t r ) δ t r = (cid:26) + µ √ K Γ ( )[ E ( t r ) − I ( t r )] (cid:27) δ t r . (A3)From time t r up to t the phase shift evolves according to,˙ δ φ = µ √ K Γ (cid:48) ( Φ )( E ( t ) − I ( t )) δ φ , (A4)where t is the minimum between the time when PRC dropsto zero and the time when the threshold is reached (in eithercase, we neglect the variation of field dynamics, since the fieldis treated as an external forcing). As a result, δ φ = e D δ φ r , (A5)where, D = (cid:90) tt r µ √ K Γ (cid:48) ( Φ )( E ( t ) − I ( t )) dt . (A6)The corresponding time shift is δ t = δ φ ˙ Φ Φ is the velocity at t . The shift δ t carries overunchanged until first the threshold Φ th = R of the time shift over one pe-riod (a sort of Floquet multiplier) can be written as R = δ t δ t r = + µ √ K Γ ( )[ E ( t r ) − I ( t r )] ˙ Φ e D (A7)This formula is substantially equivalent to Eq. (54) ofRef. [31] ( Λ ii corresponds to R ), obtained while studying asingle population under the action of α -pulses. An additionalmarginal difference is that while in Ref. [31] the single neurondynamics is described by a non uniform velocity field F ( x ) and homogeneous coupling strength, here we refer to a con-stant velocity and a phase-dependent PRC, Γ ( Φ ) .The corresponding conditional Lyapunov exponent is λ c = ln | R | T = D + ln (cid:12)(cid:12)(cid:12) [ + µ √ K Γ ( )( E ( t r ) − I ( t r ))] / ˙ Φ (cid:12)(cid:12)(cid:12) T . (A8)The formula (A8) is valid for all PRCs as long as t is replacedby T . The formula (A8) is the sum of two contributions: theformer one accounts for the linear stability of the phase evolu-tion from reset to threshold ( D / T ); the latter term arises fromthe different velocity (frequency) exhibited at threshold andat the end of the refractory period. Notice that in the limit ofshort pulses, the field amplitude at time t is effectively negli-gible, and one can thereby neglect the effect of the fields andassume ˙ Φ = B. Jarosiewicz, B. L. McNaughton, and W. E. Skaggs, “Hippocam-pal population activity during the small-amplitude irregular activitystate in the rat,” Journal of Neuroscience S. Shinomoto, H. Kim, T. Shimokawa, N. Matsuno, S. Funahashi,K. Shima, I. Fujita, H. Tamura, T. Doi, K. Kawano, and et al., “Relat-ing neuronal firing patterns to functional differentiation of cerebral cortex,”PLoS Computational Biology , e1000433 (2009). W. Truccolo, L. R. Hochberg, and J. P. Donoghue, “Collective dynamics inhuman and monkey sensorimotor cortex: predicting single neuron spikes,”Nature Neuroscience , 105–111 (2009). W. Gerstner, W. M. Kistler, R. Naud, and L. Paninski,
Neuronal Dynamics:From Single Neurons to Networks and Models of Cognition (CambridgeUniversity Press, USA, 2014). C. van Vreeswijk and H. Sompolinsky, “Chaos in neuronal networkswith balanced excitatory and inhibitory activity,” Science , 1724–1726(1996). C. v. Vreeswijk and H. Sompolinsky, “Chaotic balanced state in amodel of cortical circuits,” Neural Computation , 1321–1371 (1998),https://doi.org/10.1162/089976698300017214. N. Brunel, “Dynamics of sparsely connected networks of excitatory and in-hibitory spiking neurons,” Journal of Computational Neuroscience , 183–208 (2000). S. Luccioli and A. Politi, “Irregular collective behavior of heterogeneousneural networks,” Physical Review Letters (2010), 10.1103/Phys-RevLett.105.158104. E. Ullner, A. Politi, and A. Torcini, “Ubiquity of collective ir-regular dynamics in balanced networks of spiking neurons,” Chaos:An Interdisciplinary Journal of Nonlinear Science , 081106 (2018),https://doi.org/10.1063/1.5049902. A. Politi, E. Ullner, and A. Torcini, “Collective irregular dynamics in bal-anced networks of leaky integrate-and-fire neurons,” The European Physi-cal Journal Special Topics , 1185–1204 (2018). A. Pikovsky and M. Rosenblum, “Dynamics of globally coupled oscilla-tors: Progress and perspectives,” Chaos: An Interdisciplinary Journal ofNonlinear Science , 097616 (2015), https://doi.org/10.1063/1.4922971. E. Ullner and A. Politi, “Self-sustained irregular activity in an ensemble ofneural oscillators,” Phys. Rev. X , 011015 (2016). E. Montbrió, D. Pazó, and A. Roxin, “Macroscopic description for net-works of spiking neurons,” Phys. Rev. X , 021028 (2015). R. Zillmer, R. Livi, A. Politi, and A. Torcini, “Desynchronization in dilutedneural networks,” Phys. Rev. E , 036203 (2006). S. Jahnke, R.-M. Memmesheimer, and M. Timme, “Stable irregular dy-namics in complex neural networks,” Phys. Rev. Lett. , 048102 (2008). R. Zillmer, N. Brunel, and D. Hansel, “Very long transients, irregular fir-ing, and chaotic dynamics in networks of randomly connected inhibitoryintegrate-and-fire neurons,” Phys. Rev. E , 031909 (2009). A. Politi and A. Torcini, “Stable chaos,” in
Nonlinear dynamics and chaos:advances and perspectives , Underst. Complex Syst. (Springer, Berlin,2010) pp. 103–129. C. C. Canavier,
Encyclopedia of Computational Neuroscience , edited byD. Jaeger and R. Jung (Springer New York, New York, NY, 2013) pp. 1–11. R. Zillmer, R. Livi, A. Politi, and A. Torcini, “Stability of the splay state inpulse-coupled networks,” Phys. Rev. E , 046102 (2007). D.-P. Yang, H.-J. Zhou, and C. Zhou, “Co-emergence of multi-scale corti-cal activities of irregular firing, oscillations and avalanches achieves cost-efficient information capacity,” PLOS Computational Biology , 1–28(2017). W. Gerstner and W. Kistler,
Spiking Neuron Models: An Introduction (Cam-bridge University Press, New York, NY, USA, 2002). Afifurrahman, E. Ullner, and A. Politi, “Stability of synchronous states insparse neuronal networks,” Nonlinear Dynamics , 733–743 (2020). D. Golomb, D. Hansel, and G. Mato, “Chapter 21 mechanisms of syn-chrony of neural activity in large networks,” in
Neuro-Informatics and Neu-ral Modelling , Handbook of Biological Physics, Vol. 4, edited by F. Mossand S. Gielen (North-Holland, 2001) pp. 887 – 968. S. Ostojic, “Two types of asynchronous activity in networks of excitatoryand inhibitory spiking neurons,” Nature Neuroscience , 594–600 (2014). C. C. Canavier, “Phase response curve,” Scholarpedia , 1332 (2006). E. M. Izhikevich and B. Ermentrout, “Phase model,” Scholarpedia , 1487(2008). D. Golomb, “Neuronal synchrony measures,” Scholarpedia , 1347 (2007). E. Ullner, A. Politi, and A. Torcini, “Quantitative and qualitative analysisof asynchronous neural activity,” Phys. Rev. Research , 023103 (2020). X.-J. Wang, “Neurophysiological and computational principles of corti-cal rhythms in cognition,” Physiological Reviews , 1195–1268 (2010),pMID: 20664082, https://doi.org/10.1152/physrev.00035.2008. In the simulations, it is crucial to fix the time-step δ t at least ten timessmaller than δ p , in order to ensure that the spike times are properly handledduring the integration process. S. Olmi, A. Politi, and A. Torcini, “Linear stability in networks of pulse-coupled neurons,” Frontiers in Computational Neuroscience8