Exact mean-field theory explains the dual role of electrical synapses in collective synchronization
EExact mean-field theory explains the dual role of electrical synapses in collective synchronization
Ernest Montbri´o and Diego Paz´o Department of Information and Communication Technologies, Universitat Pompeu Fabra, 08003 Barcelona, Spain. Instituto de F´ısica de Cantabria (IFCA), CSIC-Universidad de Cantabria, 39005 Santander, Spain. (Dated: September 4, 2020)Electrical synapses play a major role in setting up neuronal synchronization, but the precise mechanismswhereby these synapses contribute to synchrony are subtle and remain elusive. To investigate these mechanismsmean-field theories for quadratic integrate-and-fire neurons with electrical synapses have been recently put for-ward. Still, the validity of these theories is controversial since they assume that the neurons produce unrealistic,symmetric spikes, ignoring the well-known impact of spike shape on synchronization. Here we show that theassumption of symmetric spikes can be relaxed in such theories. The resulting mean-field equations reveal adual role of electrical synapses: First, they equalize membrane potentials favoring the emergence of synchrony.Second, electrical synapses act as “virtual chemical synapses”, which can be either excitatory or inhibitorydepending upon the spike shape. Our results offer a precise mathematical explanation of the intricate effectof electrical synapses in collective synchronization. This reconciles previous theoretical and numerical works,and confirms the suitability of recent low-dimensional mean-field theories to investigate electrically coupledneuronal networks.
Electrical coupling via gap junctions is broadly presentacross brain areas [1–5]. There is ample experimental ev-idence that gap junctions are involved in synchronizing in-hibitory networks [6–15]. Yet, and despite the fact that elec-trical synapses are recognized to constitute a critical compo-nent of the brain, the function of these junctions and the mech-anisms whereby they contribute to synchrony are subtle andnot well understood [4, 15].Unlike chemical synapses, the transmission of action po-tentials via gap junctions greatly depends upon the overallshape of the action potential. This makes the effects of elec-trical coupling far from trivial, and the theoretical and com-putational analysis of electrically coupled networks of spik-ing neurons a difficult task. Previous theoretical studies iden-tified both the shape of the spikes and the firing frequencyof the neurons as key parameters influencing synchrony [16–22]. Yet, a precise mechanistic explanation of the role ofthese parameters in synchrony is lagging, and different worksoften provide results that are difficult to reconcile. For ex-ample, studies on homogeneous, two-neuron networks invokethe weak coupling limit and find that electrical and chemi-cal synapses combine their effect in a linear manner [17, 18].Moreover, the shape of the spikes determines whether elec-trical synapses cooperate of compete with inhibition for syn-chrony [18]. In contrast, works investigating large heteroge-neous networks suggest that electrical and inhibitory synapsesplay distinct roles: While strong electrical coupling leads tocollective synchrony [22–24], strong inhibition typically leadsto the suppression of firing and destroys synchrony [25–28].Rather than investigating the dynamics of large popula-tions of spiking neurons, an alternative and widespread the-oretical approach is to use so-called mean-field models (alsocalled firing-rate or neural-mass models) [29]. Such modelsare simplified, low-dimensional mathematical descriptions ofthe mean activity of the population but they only characterizepopulations with chemical synapses. Yet, an important the-oretical achievement linking individual and global dynamics has been recently accomplished with the advent of the Ott-Antonsen theory [30], and with its application to populationsof θ -neurons [31] and Quadratic Integrate-and-Fire (QIF) neu-rons [32]. Notably, these novel theoretical approaches readilyapply to networks with electrical synapses [33–36].Unfortunately, the mean field theories proposed in [31, 32]rely on an assumption of critical importance for electricallycoupled networks: Spikes need to have a very particular sym-metric shape that is unrealistic (see Fig. 1), and this maylargely alter the network dynamics. As we show below thisis the reason why the low-dimensional firing-rate equations(FREs) for electrically coupled networks originally derivedin [33, 34] fail to elucidate how chemical and electrical cou-plings might add their effects linearly, as found in [17, 18].Therefore, though the mean field theories in [31, 32] havebeen successfully applied to investigate networks with chem-ical synapses [28, 37–57], their suitability to investigate syn-chrony in networks with gap junctions remains questionable.In this Letter we show that the mean-field theory origi- -1000100 0 40 80 V p = 100 − V r = − (a) a = V p V r = 1 V V p = 100 − V r = − (b) a = 4 V time (ms) FIG. 1. Time series of an oscillatory QIF neuron, Eq. (1), with (a)symmetric, and (b) non-symmetric resetting rule. η = 1 , τ = 10 ms. a r X i v : . [ q - b i o . N C ] S e p nally proposed in [32] can be generalized to account for non-symmetric spikes. Remarkably, this leads to an extension ofthe FREs in [34], which provides the first precise mechanis-tic explanation of the dual role of electrical synapses in col-lective synchrony. This explanation reconciles previous dis-parate theoretical and numerical results on the intricate effectsof electrical coupling in synchrony, and consolidates the FREsfor QIF neurons as a valid and singular low-dimensional de-scription for populations of electrically coupled neurons. Quadratic integrate-and-fire neuron.–
First of all we in-troduce the QIF model, and explain the main idea to gener-alize the mean-field theory developed in [32]. The model de-scribes the evolution of the membrane potential variable V viathe first-order ordinary differential equation [58–60] τ ˙ V = V + η ; if V > V p then V ← − V r . (1)The overdot denotes derivation with respect to time, τ is themembrane’s time constant, and η represents an external cur-rent. Due to the quadratic nonlinearity in Eq. (1), V mayescape to infinity in a finite time. To prevent this, the QIFmodel incorporates a resetting rule: each time the neuron’smembrane potential V reaches the peak value V p , the neuronemits a spike and the voltage is reset to − V r (here we consider V p , V r > ), see Fig. 1.It is convenient to define a positive, real parameter charac-terizing spike asymmetry as the ratio a ≡ V p V r . (2)When a = 1 , the resetting rule of the QIF model is sym-metric, see Fig. 1(a). This symmetry is implicitly assumed inthe transformation between the QIF and the θ -neuron models,which can be formally performed in the limit V r = V p →∞ [61].The mean field theory in [32] assumes symmetric resetting, a = 1 , and then adopts the limit V r → ∞ . We next show thatthe assumption of symmetric resetting can be relaxed in thistheory, leading to a novel low-dimensional firing-rate modelwith arbitrary spike asymmetry. Population model.–
We investigate a population of QIFneurons τ ˙ V j = V j + η j + Jτ r ( t ) + g [ v ( t ) − V j ] , (3)with the resetting rule of Eq. (1). The index j = 1 , . . . , N labels all the neurons in the population, and the parameters η j represent external currents that are taken from some pre-scribed probability distribution G ( η ) .In Eq. (3) electrical coupling, of strength g ≥ , diffu-sively couples each neuron’s membrane potential with themean membrane potential v ( t ) = N (cid:80) Nj =1 V j ( t ) . Addition-ally, each neuron in the population is connected to all theother neurons via chemical synapses of strength J . Electri-cal synapses typically connect inhibitory neurons and hence J should be thought of as a negative parameter. Chemi-cal coupling is mediated by the average firing rate r ( t ) = N (cid:80) Nj =1 (cid:80) k \ t kj Next, we perform the thermo-dynamic limit, drop the indices of Eqs. (3), and definethe time-dependent conditional density ρ ( V | η, t ) , such that ρ ( V | η, t ) dV is the fraction of neurons with membrane volt-age between V and V + dV , and heterogeneity parameter η .The density ρ necessarily obeys the continuity equation τ ∂ t ρ + ∂ V (cid:8)(cid:2) V + η + Jτ r + g ( v − V ) (cid:3) ρ (cid:9) = 0 . (4)This equation can be solved resorting to the Lorentzian ansatz(related to the Ott-Antonsen ansatz [30] through a conformalmapping [32]) ρ ( V | η, t ) = 1 π x ( η, t )[ V − y ( η, t )] + x ( η, t ) , (5)which is the asymptotic shape of the density in the long timelimit assuming that V spans over the whole real line [32, 34].We note that in numerical simulations V ∈ ( − V r , V p ) , andtherefore the Lorentzian ansatz Eq. (5) is an approximationthat works progressively better as V p , V r → ∞ .Inserting Eq. (5) into Eq. (4), we get the evolution equa-tions for x ( η, t ) and y ( η, t ) . Nevertheless, before doing so, itis convenient to express the global quantities r ( t ) and v ( t ) inEq. (4) in terms of x ( η, t ) and y ( η, t ) . Regarding the mean fir-ing rate r , it is related with the width x ( η, t ) of the Lorentzianansatz [32]. Indeed, the firing rate for neurons with a given η value is the probability flux at V p (taken at infinity), whichgives the identity: r ( η, t ) = x ( η, t ) / ( πτ ) . Then, the meanfiring rate r ( t ) is given by the integral over all ηr ( t ) = 1 πτ (cid:90) ∞−∞ x ( η, t ) G ( η ) dη. (6) Mean membrane potential for general resetting rule.– Next we show that, for general resetting, the mean membranepotential v depends on the central voltages y ( η, t ) as well ason the widths x ( η, t ) in Eq. (5). This combined dependenceunderlies the dual effects of electrical synapses.The average voltage for neurons with a certain η value canbe calculated taking the following limit [62] v ( η, t ) = lim V r →∞ (cid:90) V p = aV r − V r ρ ( V | η, t ) V dV. (7)The mean field theory originally proposed in [32] considers a = 1 . However, here we relax this assumption and split theintegral in Eq. (7) into two parts: one integral with symmetricintegration limits (Cauchy principal value), and another onewith the remaining integration interval v ( η, t ) = (cid:34) p . v . (cid:90) ∞−∞ + lim V r →∞ (cid:90) aV r V r (cid:35) ρ ( V | η, t ) V dV. (8)The first integral simply yields the center of the distributionof membrane potentials, y ( η, t ) , see Eq. (5). The second inte-gral is the contribution to the mean membrane voltage due toasymmetric spike resetting, and it can be evaluated in closedform. Then, after taking the limit, we find v ( η, t ) = y ( η, t ) + ln aπ x ( η, t ) . (9)This identity pinpoints the deviation of the mean membranepotential, v ( η, t ) , with respect to the center of the distributionof voltages, y ( η, t ) , due to asymmetric spike resetting, a (cid:54) = 1 .Remarkably, this deviation is both proportional to the firingrate through x ( η, t ) , and to ln a . Therefore, for spikes with a > , the mean membrane voltage is above the center of thedistribution, whereas for spikes with < a < the meanmembrane voltage is below y ( η, t ) .For symmetric resetting, the mean membrane potential v s ( t ) of the entire population is obtained integrating y ( η, t ) over all η values: v s ( t ) ≡ (cid:90) ∞−∞ y ( η, t ) G ( η ) dη. (10)In the general case of asymmetric spike resetting, the meanmembrane voltage is obtained integrating Eq. (9) over η .Hence, using Eqs. (6) and (10), we find v ( t ) = v s ( t ) + ( τ ln a ) r ( t ) . (11) Dynamics in the “Lorentzian manifold”.– The evolutionequations for x ( η, t ) and y ( η, t ) , obtained substituting Eq. (5)into Eq. (4), can be condensed into one via the complex vari-able w ( η, t ) ≡ x ( η, t ) + iy ( η, t ) : τ ∂ t w ( η, t ) = i (cid:2) η + ( J + g ln a ) τ r − w (cid:3) + g ( iv s − w ) , (12)where r and v s are the integrals in Eqs. (6) and (10). Note thatEq. (12) is an infinite dimensional system, since each w ( η, t ) is coupled to all the other functions w ( η (cid:48) , t ) .Three results follow from Eq. (12), which are valid irre-spective of the specific form of the distribution G ( η ) : First :Only in the presence of electrical coupling ( g > ), the effectsof asymmetric spike resetting influence the collective dynam-ics of the network. Second : The term ln a weights the ‘virtualchemical coupling’ that linearly adds to the actual chemicalcoupling constant J . Moreover, the sign of ln a determines ifthis contribution is excitatory or inhibitory. Thus, in the pres-ence of electrical coupling, the effective chemical couplingconstant becomes J eff ( J, g, a ) = J + g ln a. (13) Third : Terms with electrical coupling g appear at two differentplaces in Eq. (12), indicating a dual role of gap junctions in thedynamics of w . a = 4 a = 1 a = 1 / TB CuspSNSL Sync g / √ ∆ ¯ η/ ∆ HopfHomoclinicSN/SNIC FIG. 2. Phase diagram of the FREs (14) for electrical coupling only( J = 0 ) and for three different values of the asymmetry parameter a . The region of synchronization (‘Sync’), located in the upper-rightpart of the diagram, enlarges as a grows. Symbol × indicates theparameters value (1 , . , used in Fig. 3. To lighten de diagram,homoclinic bifurcations and codimension-2 points —Cusp, Saddle-Node-Separatrix-Loop (SNSL), and Takens-Bogdanov (TB)—, areomitted in the cases a = 1 and a = 1 / . Firing-rate equations.– To gain further insight on the ef-fects of electrical coupling we adopt hereafter Lorentzian dis-tributed currents with center at ¯ η and half-width ∆ : G ( η ) =(∆ /π ) / [( η − ¯ η ) +∆ ] . In this case a maximal dimensionalityreduction is achieved, since the integrals in Eqs. (6) and (10)can be evaluated applying the residue theorem. We analyt-ically extend w ( η, t ) into complex η [32], close the integralsby an arc at infinity in the complex half-plane Im ( η ) < [63],and obtain: πτ r ( t ) + iv s = w (¯ η − i ∆ , t ) . Then, we evaluateEq. (12) at η = ¯ η − i ∆ , obtaining the FREs τ ˙ r = ∆ τ π + 2 rv s − gr, (14a) τ ˙ v s = v s + ¯ η − ( πτ r ) + ( J + g ln a ) τ r. (14b)This system of two ordinary differential equations for themean firing rate r and for the auxiliary variable v s describesthe dynamics of the ensemble exactly in the limit of infinite N , V p and V r . [64]. The dual role of the electrical coupling g is clearly described by the FREs: First, electrical couplingenters in Eq. (14a), reducing the firing rate, and hence reduc-ing the width of the distribution of membrane potentials. Thiscontribution homogenizes membrane potentials and, for largeenough g , favors the emergence of collective synchronization[65]. Second, electrical coupling enters in Eq. (14b), con-tributing as a virtual chemical coupling. Depending on thevalue of a , electrical and chemical synapses cooperate or com-pete for synchrony.The relevant parameters of Eqs. (14) are ¯ η , the couplingconstants J and g , and the spike asymmetry a (one can assume ∆ = 1 and τ = 10 ms without loss of generality). Addition-ally, parameters a and J are both acting in the same term ofEq. (14) and hence they can not produce qualitatively different (a) a = 1 / , V p = 100 r ( H z ) Firing rate modelQIF network (b) a = 1 , V p = 100 0250500 0 60 120 (c) a = 4 , V p = 1000250 0 60 120 (d) a = 1 / , V p = 1000 r ( H z ) time (ms) (e) a = 1 , V p = 1000 time (ms) (f) a = 4 , V p = 1000 time (ms) FIG. 3. Time series of r ( t ) for FREs (14) (red) and network Eqs. (3) with N = 10 (black). Peak values are (a,b,c) V p = 100 , and (d,e,f) V p = 1000 . Three values of the spike asymmetry parameter are used: (a, d) a = 1 / , (b,e) a = 1 , (c,f) a = 4 . Parameters correspond to thesymbol × in Fig. 2: J = 0 , g = 2 . , ¯ η = ∆ = 1 , and τ = 10 ms. We used the explicit Euler scheme with dt = 10 − ms and τ s = 10 − ms. dynamical behaviors. Accordingly, we consider J = 0 here-after and borrow the results of the bifurcation analysis in [34],which can be directly applied here replacing the chemical cou-pling constant J by the effective chemical coupling Eq. (13)[66]. The phase diagram in Fig. 2 shows the bifurcation lociof the FREs for three values of a . Three bifurcations, Hopf(red), Saddle-Node (black) and Homoclinic (green), tangen-tially meet at the Takens-Bogdanov (TB) point, (¯ η/ ∆) TB = − π ln a, (15)around which the phase diagram organizes. The homocliniccurve (shown only for a = 4 ) moves parallel to the Hopfline for a while, and then tangentially meets the upper branchof the Saddle-Node (SN) bifurcating line at a saddle-node-separatrix-loop (SNSL) point. At this point, the SN bound-ary becomes a SN on the Invariant Circle (SNIC) boundary(black) that together with the Hopf and homoclinic lines en-closes the region of synchronization (Sync) where collectiveoscillations occur. Note that two small regions of bistabilityare located in the region limited by the three codimension-twopoints: Cusp, TB, and SNSL.In Fig. 2 the SNIC bifurcation lines have a vertical asymp-tote precisely at the same ¯ η value than the TB point Eq. (15).For < a < , the asymptote is located at positive values of ¯ η , and shifts to the left as a increases. Hence, synchroniza-tion is hindered for spikes with < a < , and favored forspikes with a > [67]. In Fig. 3 we depict the time seriesof the firing rate in numerical simulations of both the networkEqs. (3) and the FREs Eqs. (14), for three values of the pa-rameter a = { / , , } . In the first row, we used spikes with V p = 100 in Eqs. (3). Though the agreement between themodels is not perfect, the mean field Eqs. (14) accurately pre-dict the emergence of oscillations in Eqs. (3) as a is increased.Finally, for V p = 1000 (second row), the agreement greatlyimproves showing that the dynamics of the network Eqs. (3) converges to that of the FREs for large values of V p . Conclusions.– We showed that the spike resetting rule ofthe QIF model can be incorporated in the mean-field theoryproposed in [32]. This extension of the theory reveals a non-trivial dependence of the mean membrane voltage on the fir-ing rate —see Eq. (11)— which turns out to be crucial fordeciphering the dual role of electrical synapses in synchrony,and to reconcile previous theoretical and numerical work. Asthe dynamics in the “Lorentzian manifold” Eq. (12) as wellas the FREs (14) are valid for arbitrary coupling strengths,we conclude that electrical and chemical coupling do not sim-ply add linearly for weak coupling, as suggested by previousworks, see e.g. [23]. In view that traditional mean-field theo-ries for neural networks only encompass chemically coupledensembles, we regard the mean-field theory presented here asa unique tool to investigate the interplay between chemicaland electrical synapses.We acknowledge support by the Agencia Estatal de Inves-tigaci´on and Fondo Europeo de Desarrollo Regional underProject No. FIS2016-74957-P (AEI/FEDER, EU). [1] M. Galarreta and S. Hestrin, Nature Reviews Neuroscience ,425 (2001).[2] M. V. Bennett and R. Zukin, Neuron , 495 (2004).[3] B. W. Connors and M. A. Long, Annual Review of Neuro-science , 393 (2004).[4] J. I. Nagy, A. E. Pereda, and J. E. Rash, Biochimica et Bio-physica Acta (BBA) - Biomembranes , 102 (2018).[5] R. D. Traub, M. A. Whittington, R. Guti´errez, and A. Draguhn,Cell and Tissue Research , 671 (2018).[6] A. Draguhn, R. Traub, D. Schmitz, and J. Jefferys, Nature ,189 (1998).[7] P. Mann-Metzer and Y. Yarom, Journal of Neuroscience ,3298 (1999). [8] J. L. Perez Velazquez and P. L. Carlen, Trends in Neurosciences , 68 (2000).[9] G. Tam´as, E. H. Buhl, A. L¨orincz, and P. Somogyi, Natureneuroscience , 366 (2000).[10] M. R. Deans, J. R. Gibson, C. Sellitto, B. W. Connors, andD. L. Paul, Neuron , 477 (2001).[11] S. G. Hormuzdi, I. Pais, F. E. LeBeau, S. K. Towers, A. Rozov,E. H. Buhl, M. A. Whittington, and H. Monyer, Neuron ,487 (2001).[12] S. Stagkourakis, C. T. P´erez, A. Hellysaz, R. Ammari, andC. Broberger, Elife , e33144 (2018).[13] C. Bou-Flores and A. J. Berger, Journal of Neurophysiology ,1543 (2001), pMID: 11287478.[14] K. Vervaeke, A. Lrincz, P. Gleeson, M. Farinella, Z. Nusser,and R. A. Silver, Neuron , 435 (2010).[15] B. W. Connors, Developmental Neurobiology , 610 (2017).[16] C. C. Chow and N. Kopell, Neural Computation , 1643(2000).[17] T. J. Lewis and J. Rinzel, Journal of Computational Neuro-science , 283 (2003).[18] B. Pfeuty, G. Mato, D. Golomb, and D. Hansel, Neural Com-putation , 633 (2005).[19] B. Pfeuty, G. Mato, D. Golomb, and D. Hansel, Journal ofNeuroscience , 6280 (2003).[20] T. J. Lewis and F. K. Skinner, “Understanding activity inelectrically coupled networks using PRCs and the theory ofweakly coupled oscillators,” in Phase Response Curves in Neu-roscience: Theory, Experiment, and Analysis , edited by N. W.Schultheiss, A. A. Prinz, and R. J. Butera (Springer New York,New York, NY, 2012) pp. 329–359.[21] C. B¨orgers, An introduction to modeling neuronal dynamics ,Vol. 66 (Springer, 2017).[22] S. Ostojic, N. Brunel, and V. Hakim, Journal of ComputationalNeuroscience , 369 (2009).[23] N. Kopell and B. Ermentrout, Proceedings of the NationalAcademy of Sciences , 15482 (2004).[24] B. Pfeuty, D. Golomb, G. Mato, and D. Hansel, Frontiers inComputational Neuroscience , 8 (2007).[25] X.-J. Wang and G. Buzs´aki, The Journal of Neuroscience ,6402 (1996).[26] J. A. White, C. C. Chow, J. Rit, C. Soto-Trevi˜no, and N. Kopell,Journal of Computational Neuroscience , 5 (1998).[27] P. Tiesinga and J. V. Jos´e, Network: Computation in NeuralSystems , 1 (2000).[28] F. Devalle, A. Roxin, and E. Montbri´o, PLoS ComputationalBiology , 1 (2017).[29] H. R. Wilson and J. D. Cowan, Biophys. J. , 1 (1972).[30] E. Ott and T. M. Antonsen, Chaos , 037113 (2008).[31] T. B. Luke, E. Barreto, and P. So, Neural Comput. , 3207(2013).[32] E. Montbri´o, D. Paz´o, and A. Roxin, Phys. Rev. X , 021028(2015).[33] C. R. Laing, SIAM Journal on Applied Dynamical Systems ,1899 (2015).[34] B. Pietras, F. Devalle, A. Roxin, A. Daffertshofer, and E. Mont-bri´o, Phys. Rev. E , 042412 (2019).[35] A. Daffertshofer and B. Pietras, “Phase synchronization in neu-ral systems,” in Encyclopedia of Complexity and Systems Sci-ence , edited by R. A. Meyers (Springer Berlin Heidelberg,Berlin, Heidelberg, 2020) pp. 1–14.[36] ´A. Byrne, J. Ross, R. Nicks, and S. Coombes, bioRxiv (2020).[37] P. So, T. B. Luke, and E. Barreto, Physica D , 16 (2014).[38] C. R. Laing, Phys. Rev. E , 010901 (2014).[39] J. Roulet and G. B. Mindlin, Chaos: An Interdisciplinary Jour- nal of Nonlinear Science , 093104 (2016).[40] D. Paz´o and E. Montbri´o, Phys. Rev. Lett. , 238101 (2016).[41] I. Ratas and K. Pyragas, Phys. Rev. E , 032215 (2016).[42] J. M. Esnaola-Acebes, A. Roxin, D. Avitabile, and E. Mont-bri´o, Phys. Rev. E , 052407 (2017).[43] G. Dumont, G. B. Ermentrout, and B. Gutkin, Phys. Rev. E ,042311 (2017).[44] F. Devalle, E. Montbri´o, and D. Paz´o, Phys. Rev. E , 042214(2018).[45] H. Schmidt, D. Avitabile, E. Montbri´o, and A. Roxin, PLoSComputational Biology , 1 (2018).[46] M. di Volo and A. Torcini, Phys. Rev. Lett. , 128301 (2018).[47] I. Ratas and K. Pyragas, Phys. Rev. E , 052224 (2018).[48] S. Coombes and ´A. Byrne, in Nonlinear Dynamics in Com-putational Neuroscience , edited by F. Corinto and A. Torcini(Springer International Publishing, Cham, 2019) pp. 1–16.[49] A. Byrne, D. Avitabile, and S. Coombes, Phys. Rev. E ,012313 (2019).[50] G. Dumont and B. Gutkin, PLOS Computational Biology , 1(2019).[51] S. Keeley, ´A. Byrne, A. Fenton, and J. Rinzel, Journal of Neu-rophysiology , 2181 (2019).[52] H. Schmidt and D. Avitabile, Chaos: An Interdisciplinary Jour-nal of Nonlinear Science , 033133 (2020).[53] R. Gast, H. Schmidt, and T. R. Knoesche, Neural Computation , 1615 (2020).[54] M. Segneri, H. Bi, S. Olmi, and A. Torcini, Frontiers in Com-putational Neuroscience , 47 (2020).[55] A. Ceni, S. Olmi, A. Torcini, and D. Angulo-Garcia, Chaos:An Interdisciplinary Journal of Nonlinear Science , 053121(2020).[56] H. Bi, M. Segneri, M. di Volo, and A. Torcini, Phys. Rev. Re-search , 013042 (2020).[57] C. Bick, M. Goodfellow, C. R. Laing, and E. A. Martens, TheJournal of Mathematical Neuroscience , 1 (2020).[58] E. M. Izhikevich, Dynamical Systems in Neuroscience (TheMIT Press, Cambridge, Massachusetts, 2007).[59] P. Latham, B. Richmond, P. Nelson, and S. Nirenberg, Journalof Neurophysiology , 808 (2000).[60] D. Hansel and G. Mato, Neural Computation , 1 (2003).[61] B. Ermentrout and N. Kopell, SIAM J. Appl. Math. , 233(1986).[62] The integration limits cannot be simply replaced by ±∞ , dueto the slow decay of the Lorentzian distribution ρ ( V | η, t ) V ∼ V − as | V | → ∞ .[63] The integrals along the arc at infinity are zero, since accordingto Eq. (12) we have: | w | ∼ (cid:112) | η | . Accordingly, e.g. for Eq. (6),we have (cid:82) − π Re( w ) G ( η ) | η | ie iφ dφ ∼ | η | / / | η | → ; like-wise for Eq. (10).[64] It is also possible to write Eqs. (14) using the variables r and v .Here we use the auxiliary variable v s instead of v since in thisform the FREs are simpler to interpret.[65] For symmetric resetting, a = 1 , the mean membrane voltage is v = v s and the FREs derived in [34] are recovered. In this casechemical and electrical coupling act independently, in differentterms of Eqs. (14).[66] The weakly nonlinear analysis performed in [34] also applieshere, and confirms that —in contrast to [22]— the Hopf bifur-cation remains always supercritical for any value of a .[67] The asymptote of the SNIC boundary is located at ¯ η < for a > . Therefore, even if a majority of the neurons are intrin-sically quiescent, electrical coupling is able to trigger a self-sustained collective oscillation. Similar effects have been re- ported for pairs of electrically coupled neurons [68, 69].[68] Y. Manor, J. Rinzel, I. Segev, and Y. Yarom, Journal of Neuro- physiology , 2736 (1997).[69] T. Kepler, E. Marder, and L. Abbott, Science248