An exact firing rate model reveals the differential effects of chemical versus electrical synapses in spiking networks
Bastian Pietras, Federico Devalle, Alex Roxin, Andreas Daffertshofer, Ernest Montbrió
AAn exact firing rate model reveals the differential effects of chemical versus electrical synapses inspiking networks
Bastian Pietras,
1, 2, 3, 4
Federico Devalle,
5, 2
Alex Roxin,
6, 7
Andreas Daffertshofer, and Ernest Montbri´o Faculty of Behavioural and Movement Sciences, Amsterdam Movement Sciences & Institute of Brain and Behavior Amsterdam,Vrije Universiteit Amsterdam, van der Boechorststraat 9, Amsterdam 1081 BT, The Netherlands. Department of Physics, Lancaster University, Lancaster LA1 4YB, United Kingdom. Institute of Mathematics, Technical University Berlin, 10623 Berlin, Germany. Bernstein Center for Computational Neuroscience Berlin, 10115 Berlin, Germany. Department of Information and Communication Technologies, Universitat Pompeu Fabra, 08003 Barcelona, Spain. Centre de Recerca Matem`atica, Campus de Bellaterra, Edifici C, 08193 Bellaterra (Barcelona), Spain. Barcelona Graduate School of Mathematics, Barcelona, Spain. (Dated: September 23, 2019)Chemical and electrical synapses shape the dynamics of neuronal networks. Numerous theoretical studieshave investigated how each of these types of synapses contributes to the generation of neuronal oscillations,but their combined effect is less understood. This limitation is further magnified by the impossibility of tradi-tional neuronal mean-field models —also known as firing rate models, or firing rate equations— to account forelectrical synapses. Here we introduce a novel firing rate model that exactly describes the mean field dynamicsof heterogeneous populations of quadratic integrate-and-fire (QIF) neurons with both chemical and electricalsynapses. The mathematical analysis of the firing rate model reveals a well-established bifurcation scenario fornetworks with chemical synapses, characterized by a codimension-2 Cusp point and persistent states for strongrecurrent excitatory coupling. The inclusion of electrical coupling generally implies neuronal synchrony byvirtue of a supercritical Hopf bifurcation. This transforms the Cusp scenario into a bifurcation scenario char-acterized by three codimension-2 points (Cusp, Takens-Bogdanov, and Saddle-Node Separatrix Loop), whichgreatly reduces the possibility for persistent states. This is generic for heterogeneous QIF networks with bothchemical and electrical coupling. Our results agree with several numerical studies on the dynamics of largenetworks of heterogeneous spiking neurons with electrical and chemical coupling.
PACS numbers: 05.45.Xt
I. INTRODUCTION
Collective oscillations and synchrony are prominent fea-tures of neuronal circuits, and are fundamental for thewell-timed coordination of neuronal activity. Such oscil-lations are profoundly shaped by the presence of chemicalsynapses [1]. An increasing number of experimental stud-ies indicate both the prevalence and functional importance ofelectrical synapses (formed by gap junctions between neu-rons) in many diverse regions of central nervous systems, es-pecially in inhibitory interneurons [2–4]. Electrical synapsesparticipate in mediating synchronization of neuronal networkactivity [5, 6], suggesting that electrical interaction may beinterrelated with the generation of oscillations via chemicaltransmission.The mechanisms by which chemical synapses mediatelarge-scale synchronous activity have been extensively inves-tigated, see e.g. [1, 7]. However, only a few studies ad-dressed the synchronization of large networks in which neu-rons are not only interacting via excitation and/or inhibition,but also via electrical synapses [8–20]. This limited theo-retical progress for networks of electrically coupled neurons,compared to chemically coupled networks, is magnified dueto the technical challenges faced when developing simplifiedmean field models —often called firing rate models, or fir-ing rate equations (FRE)— for networks involving electricalsynapses. While firing rate models turned out to be very use-ful to explain key aspects of the dynamics of spiking neuron networks with chemical synapses [21–33], it remains an openquestion whether there are similar simplified mean field theo-ries for networks involving electrical interactions.Recently, a novel method has been found to exactly deriveFRE for populations of heterogeneous quadratic integrate-and-fire (QIF) neurons with chemical coupling [34]. Themethod, related to the so-called Ott-Antonsen ansatz [35–41],allows to obtain exact, low-dimensional firing rate equationsfor ensembles of QIF neurons, see also [42–44]. The FRE forQIF neurons have been used to investigate numerous problemsregarding the dynamics of networks of chemically coupledQIF neurons [45–64]. Remarkably, previous work has alsosought to apply this approach to networks with both chemicaland electrical coupling [65]. However, in [65], the electricalcoupling has been treated by making use of an approximationwhich renders the resulting FRE analytically intractable. Webuild on this previous work and derive a set of FRE for net-works with chemical and electrical coupling, but without theneed for any approximation. The resulting system is not onlyanalytically tractable but also allows, in a unified framework,for carrying out a complete analysis of the possible dynamicsand bifurcations of networks with mixed chemical and elec-trical synapses. In
Appendix B we show that our exact FREare recovered by appropriately relaxing the approximation in-voked in [65].The structure of the paper is as follows: In Section II, wedescribe the spiking neuron network under investigation, andbriefly illustrate the impact of electrical coupling in the dy-namics of two nonidentical QIF neurons. In Section III, we in- a r X i v : . [ n li n . AO ] S e p troduce the FRE corresponding to the thermodynamic limit ofthe QIF network. The detailed derivation is performed in Ap-pendix A . In Section IV, we perform a comparative analysis ofthe fixed points and their bifurcations in networks with electri-cal coupling vs. networks with chemical coupling. Finally, weinvestigate the dynamics of a QIF network with both electri-cal and chemical synapses and demonstrate that the presenceof electrical coupling critically determines the bifurcation sce-nario of the neuronal network. Finally, we discuss our resultsin Section V.
II. QUADRATIC INTEGRATE-AND-FIRE NEURONSWITH ELECTRICAL AND CHEMICAL SYNAPSES
We consider a large population of globally electrically andchemically-coupled QIF neurons, with membrane potentials { V j } j =1 ,...,N and N (cid:29) . Their dynamics reads τ ˙ V j = V j + η j + g ( v − V j ) + Jτ s, (1)where τ denotes the cells’ common membrane time constant,and parameter η j represents an external input current flowinginto cell j . To model the action potential, the continuous dy-namics Eq. (1) is supplemented by a discrete resetting rule.Here, we assume that if V j reaches infinity, neuron j emitsa spike and its membrane potential is reset to minus infin-ity [79]. The mean membrane voltage v = 1 N N (cid:88) k =1 V k , to which all cells are diffusively coupled with strength g ≥ ,mediates the electrical coupling. The constant J quantifiesthe coupling strength of chemical synapses. The coupling viachemical synapses is mediated by the mean synaptic activa-tion function s ( t ) = 1 N N (cid:88) j =1 τ s t (cid:90) t − τ s (cid:88) k δ (cid:0) t (cid:48) − t kj (cid:1) dt (cid:48) , (2)where t kj denotes the time of the k -th spike of the j -th neuron, δ ( t ) is the Dirac delta function, and τ s is a synaptic time con-stant [80]. The synaptic weight J can be positive or negativedepending on whether the chemical synapses are excitatory orinhibitory, respectivelyIn the absence of coupling, J = g = 0 , the QIF neuronsare either quiescent ( η i < ), or oscillatory ( η i > ) withfrequency f i = 1 τ π √ η i . (3)These two dynamical regimes of individual neurons are con-nected by a saddle-node on the invariant circle (SNIC) bifur-cation, which occurs when η i = 0 , with f i = 0 .Electrical coupling tends to equalize the membrane poten-tials of the neurons they connect and may favor synchrony. − (a) V − − (b) V time (ms) FIG. 1: Strong electrical coupling suppresses oscillatory activity fornegative mean currents, ¯ η < . The panels show the time seriesof the membrane voltages of N = 2 electrically coupled QIF neu-rons. (a) Self-oscillatory neuron with η = π ; (b) Quiescent neuronwith η = − π . The two neurons are either uncoupled (black thincurves, g = 0 ), weakly coupled (red curves, g = 1 ), or stronglycoupled (blue thick curves, g = 6 ). We used τ = 10 ms and J = 0 . Yet, if a large fraction of cells in the network is quiescent,gap junctions may suppress oscillations and neural synchrony.Next we illustrate this phenomenon for two nonidentical QIFneurons that are coupled via a gap junction [81]. The resultsof this analysis will later be useful to understand some aspectsof the dynamics of a large network of electrically coupled QIFneurons.
A. Strong coupling limit of two electrically coupled QIFneurons
We consider a network of N = 2 nonidentical QIF neu-rons with dynamics Eq. (1). The neurons are coupled via gapjunctions only, i.e. J = 0 but g > . We are interested in thestrong coupling limit g (cid:29) when η > and η < . InFig. 1 we depict the corresponding time series of cell 1 (panela) and cell 2 (panel b). Black thin curves correspond to thedynamics of the uncoupled ( g = 0 ) cells: cell 1 fires period-ically, while cell 2 remains quiescent. When the neurons areelectrically coupled (red curves), the membrane voltage of cell2 displays a series of so-called ‘spikelets’ [82]. Moreover, theelectrical interaction brings cell 1 closer to its firing thresholdand, hence, its frequency f is reduced. When g is increasedfurther, cell 1 becomes quiescent (blue thick curves).Although analyzing the dynamics of the two cells for arbi-trary coupling strength g is a challenge, there exists a simpleand general result valid in the large g limit, and of relevancefor the large- N analysis carried out below. Indeed, for large g , the dynamics of the N = 2 network simply depends on thesign of the mean current [66] ¯ η = η + η . For ¯ η > , the quiescent cell eventually becomes self-oscillatory as g is increased from zero. By contrast, for ¯ η < ,the oscillatory cell eventually turns quiescent in the strongcoupling limit; see the blue lines in Fig. 1 [83]. III. FIRING RATE MODEL
In the following, we introduce the FRE corresponding tothe thermodynamic limit of Eqs. (1). The detailed derivationof the model closely follows the lines of [34] and is given in
Appendix A .For N → ∞ , one can drop the indices in Eq. (1) and de-fine a density function ρ such that ρ ( V | η, t ) dV denotes thefraction of neurons with membrane potentials between V and V + dV and parameter η at time t . In the limit of instanta-neous synaptic processing, i.e. for τ s → , Eq. (2) reduces to s ( t ) = r ( t ) with r ( t ) being the population-mean firing rate. Ifthe external currents are distributed according to a Lorentziandistribution centered around η = ¯ η with half-width ∆ , L ∆ , ¯ η ( η ) = 1 π ∆( η − ¯ η ) + ∆ , (4)we find that the asymptotic mean-field dynamics evolves ac-cording to the following FRE [84] τ ˙ r = ∆ τπ + 2 rv − gr, (5a) τ ˙ v = v + ¯ η − ( πτ r ) + Jτ r. (5b)The variables r and v are the mean firing rate and mean mem-brane potential, respectively. They determine the total volt-age density for the network Eq. (1), which turns out to bea Lorentzian distribution centered at v ( t ) and of half-width πτ r ( t ) , ρ ( V, t ) = 1 π πτ r ( t )[ V − v ( t )] + [ πτ r ( t )] . (6)The structure of the FRE Eqs. (5) reveals an interesting fea-ture: Electrical coupling is solely mediated by the firing ratethrough the negative feedback term − gr in the r -dynamicsEq. (5a), and not by membrane potential differences [85].That is, electrical coupling leads to a narrowing of the voltagedistribution Eq. (6), i.e. a decrease in firing rate. This confirmsour initial sketch that electrical coupling tends to equalize theneurons’ membrane potentials and, under suitable conditions,this may promote synchrony. By contrast, chemical couplingshifts the center of the distribution Eq. (6) of voltages via thefeedback term Jr in the v -dynamics Eq. (5b). The followingphase plane and bifurcation analysis of the FRE (5) allows forunderstanding the collective dynamics of the QIF network. IV. ANALYSIS OF THE FIRING RATE EQUATIONSA. Electrical vs. chemical coupling
In the absence of chemical coupling, our previous discus-sion of the case N = 2 hints at two distinct dynamical regimesfor positive and negative values of ¯ η . With respect to the − (a) g/ √ ∆ = 0 r ( π τ / √ ∆ ) v/ √ ∆ − (b) v/ √ ∆ FIG. 2: The sign of the mean current ¯ η determines the behavior ofthe fixed points of the FRE (5) with electrical coupling only, J = 0 .The panels show the nullclines of the FRE (5) with only electricalcoupling ( J = 0 ) for negative ( ¯ η/ ∆ = − ) and positive ( ¯ η/ ∆ =0 . ) values of ¯ η (panel a and b, respectively), and g/ √ ∆ = 0 , , .The black points correspond to the intersections of r -nullclines ( ˙ r =0 , red) and v -nullclines ( ˙ v = 0 , gray) and are fixed points of Eqs. (5). fixed points ( r ∗ , v ∗ ) of the FRE (5) for J = 0 , we find the v -nullcline to be πτ r = (cid:112) v + ¯ η. Note that if ¯ η is negative, there exists a range of ‘forbid-den’ values of v . Fig. 2(a) shows the nullclines for ¯ η < and for different values of the ratio g/ √ ∆ . Since the ma-jority of the neurons are quiescent, an increase in couplingstrength g causes active neurons to reduce firing, which leadsto a progressive decrease of the firing rate r ∗ . By contrast, inFig. 2(b) the majority of the cells are self-oscillatory, ¯ η > ,and strong electrical coupling forces quiescent neurons to fire.This yields an increase of v ∗ . Interestingly, the firing rate r ∗ isa non-monotonic function of g/ √ ∆ : While v ∗ remains nega-tive, the voltages are pushed to subthreshold values, decreas-ing the firing rate. This behavior is reverted when v ∗ becomespositive and all voltages are pushed towards values above thefiring threshold. The different behaviors of Eqs. (5) with elec-trical coupling for positive and negative values of ¯ η are clearlyrevealed in the corresponding bifurcation diagrams shown inFigs. 3(a,c).The case of networks with only chemical coupling, g = 0 ,is simpler [34]. The bifurcation diagram depicted in Fig. 3(d)shows that v ∗ remains always negative and converges asymp-totically to zero as J increases. The firing rate r ∗ , depicted inFig. 3(b), also increases with J . For ¯ η < and strong recur-rent excitatory coupling, the system undergoes a cusp bifurca-tion and two saddle-node (SN) bifurcations are created. Thisimplies the existence of a parameter regime where a persistent,high-activity state (stable focus) coexists with a low-activitystate (stable node) —see Fig. 7(a), and [34]. This coexistencebetween persistent and low-activity states also occurs in net-works with electrical synapses, but it is located in a very smallregion of parameters as we show below, see Fig. 6(b).We next explore the linear stability of the fixed points ofEqs. (5), see also Appendix C . We find that a Hopf bifurcation
02 0 4¯ η/ ∆ = 2¯ η/ ∆ = 0 . η/ ∆ = − Electrical coupling(a) r ∗ ( π τ / √ ∆ ) − Chemical coupling(b) −
202 0 4¯ η/ ∆ = 2¯ η/ ∆ = 0 . η/ ∆ = − (c) v ∗ / √ ∆ g/ √ ∆ − − (d) J/ ( π √ ∆) FIG. 3: The bifurcation diagrams of the FRE (5) for networks with(left) electrical and (right) chemical coupling are qualitatively differ-ent. Top panels show the scaled firing rate r ∗ / ( πτ √ ∆) versus theratio of (a) electrical g/ √ ∆ and (b) chemical J/ ( π √ ∆) synapticstrengths. Panels (c,d) show the same bifurcation diagrams for thescaled mean membrane potential v ∗ / √ ∆ . occurs along the boundary (cid:16) ¯ η ∆ (cid:17) H = 4∆ g − g − Jπg , (7)with frequency f H = 1 πτ (cid:115) ¯ η + ∆ π Jg . (8)The Hopf boundary Eq. (7) is depicted in red in the phasediagrams of Figs. 6, 7. Note that ¯ η/ ∆ → + ∞ as g → according to Eq. (7), which indicates that electrical couplingis a necessary ingredient for the Hopf bifurcation to exist [86].To confirm the presence of collective oscillations in theoriginal network of electrically coupled QIF neurons withdynamics Eq. (1), we carried out numerical simulations andcompared them with those of the FRE (5). Fig. 4 showsthe time series of the firing rate in the full and in the re-duced system, which display a very good agreement. Inpanel (a) we considered a network with electrical couplingonly. The frequency of the oscillations, f ≈ . Hz, isclose to the theoretical value at criticality, given by Eq. (8): f H = 100 /π ≈ . Hz. Therefore, in absence of chemi-cal coupling and near the Hopf bifurcation, the frequency ofthe oscillations is almost independent of the coupling strength g and closely follows Eq. (8). To further test the validity ofEq. (8) far from criticality, we numerically evaluated the fre-quency of the limit cycle of the FRE (5) (black solid line,Fig. 5) as the the coupling strength g is increased from theHopf bifurcation (at g H ≈ . ). The black dotted line corre-sponds to the Hopf frequency Eq. (8). We find that the fre-quency of the limit cycle remains close to this for a broadrange of g values. (a) r ( H z ) (b) r ( H z ) time (ms) FIG. 4: Electrical coupling promotes collective synchrony. The in-clusion of inhibitory coupling degrades synchrony and slows downoscillations. The figure shows the time series of the mean firingrate r of a network of N = 10 QIF neurons with dynamics (1)(black) and of the FRE (5) (red). Panel (a) shows collective oscilla-tions (frequency f ≈ . Hz) of a network with gap junctions only( J = 0 ). Panel (b) corresponds to a network with both gap junctionsand inhibitory chemical coupling ( J = − π ), which both slows down( f ≈ . Hz) and reduces the amplitude of collective oscillations.Parameters: g = 3 , ¯ η = 1 , τ = 10 ms and ∆ = 1 . Hopf instability in networks of electrically coupled QIFneurons occurs like the transition to synchronization in theKuramoto model of coupled phase oscillators [67]. Consider-ing J = 0 , we find the main features of the Kuramoto tran-sition to collective synchronization: ( i ) In the limit of weakelectrical coupling g , the Hopf boundary Eq. (7) can be writ-ten as g H ≈ √ ¯ η . (9)For ¯ η = 1 , Eq. (9) coincides with Kuramoto’s critical cou-pling for synchrony. ( ii ) As previously discussed, macro-scopic oscillations emerge with a frequency determined by themost likely value of the natural frequencies in the network, seeEq. (3). For the case of the Lorentzian distribution of currents,Eq. (4), the most likely value of the frequency is ¯ f = √ ¯ ηπτ . (10)( iii ) The Hopf bifurcation is always supercritical; cf. Ap-pendix D . Taken together, for ¯ η > and given a certainlevel of heterogeneity ∆ , synchronization occurs —at a criti-cal coupling approximately given by Eq. (9)— with the nucle-ation of a small cluster of oscillators with natural frequenciesEq. (3) near ¯ f . As electrical coupling g is further increased,more and more oscillators become entrained to the frequency ¯ f , resulting in a continuous and monotonous increase in theamplitude of the oscillations. This transition is in contrast tothat of networks with inhibitory coupling and synaptic kinet-ics and/or delays, where synchrony is only achieved for weakheterogeneity and weak coupling, see, e.g., [48, 49].The phase diagram depicted in Fig. 6 characterizes the dy-namics of the firing rate model Eq. (5) with only electrical J = π J = 0 J = − π f ( H z ) g FIG. 5: In the absence of chemical coupling (black), electrical cou-pling g has little effect on the frequency of the oscillations Excita-tory (red)/Inhibitory (blue) coupling speeds-up/slows-down collec-tive oscillations. These effects tend to disappear for strong electricalcoupling. The figure shows the frequency of the oscillations f asa function of the strength of electrical coupling g in networks withExcitation, J = π , Inhibition, J = − π , and without chemical cou-pling, J = 0 . Symbols ( ◦ ) are frequencies obtained from numericalsimulations of a network of N = 10 QIF neurons Eqs. (1). Solidlines are numerically obtained frequencies from the FRE (5). Dottedlines correspond to the Hopf frequency given by Eq. (8). Parameters: ¯ η = 1 , τ = 10 ms, and ∆ = 1 . coupling. The red curve corresponds to the Hopf bifurcationline given by Eq. (7). According to Eq. (8), the frequency ofthe collective oscillations approaches zero as ¯ η → . This in-dicates that the Hopf line ends in a Takens-Bogdanov (TB) bi-furcation at ¯ η = 0 , see Fig. 6(b). At this codimension-2 point,the Hopf boundary tangentially meets a SN bifurcation and ahomoclinic bifurcation. The homoclinic line moves parallelto the Hopf line for a while, it makes a sharp backward turnand then tangentially joins onto the upper branch of the SN bi-furcation curve (two branches of SN bifurcations are createdat the Cusp point), at a saddle-node-separatrix-loop (SNSL)point. At this point the SN boundary becomes a SNIC bound-ary that, together with the Hopf and homoclinic lines, enclosesthe region of synchronization (Sync) featuring collective os-cillations. Note that in Fig. 6(b) we encounter a very small re-gion of bistability between a Low-Activity State (LAS, node)and a persistent state (focus). Electrical coupling destabilizesthe persistent state almost immediately after the SN line, lead-ing to another small region of bistability between LAS and asmall amplitude limit cycle (Sync) —which disappears in thehomoclinic bifurcation.Finally, the SNIC curve asymptotically approaches ¯ η = 0 as g/ √ ∆ → ∞ (as suggested by the N = 2 analysis inSection II A). In this limit, all neurons are strongly coupled( g → ∞ ) and/or are nearly identical ( ∆ → ) so that theybehave as a single QIF neuron with input current ¯ η [87].
048 -2 0 2
Async Sync S N I C H o p f State N o d e F o c u s (a) g / √ ∆ ¯ η/ ∆ TB CuspSNSL h o m LAS + Persistent StateLAS + Sync (b) ¯ η/ ∆ FIG. 6: The phase diagram of the FRE (5) for electrical couplingonly ( J = 0 ) is characterized by the presence of three codimension-2 bifurcation points (Takens-Bogdanov, TB; Saddle-Node SeparatixLoop, SNSL; Cusp), all located at ¯ η ≥ . The region of syn-chronization (Sync) is limited by supercritical Hopf (red), SNIC(black), and homoclinic (green) bifurcations. Dashed line: Focus-Node boundary of the Asynchronous State. Panel (b) Enlargementof the region near the three codimension-2 points. There are twosmall regions of bistability between Asynchronous, low-activity-states (LAS) and Asynchronous Persistent States, and between LASand Sync. Two Saddle-Node (SN) bifurcations are created at a Cusppoint, at (1 / (3 √ , √ / / ) ≈ (0 . , . . The upper SNline meets the homoclinic (hom) bifurcation in a SNSL point. At thispoint the upper SN becomes a SNIC bifurcation. The other SN bifur-cation tangentially meets the homoclinic and the Hopf lines at a TBpoint, at (0 , √ ≈ (0 , . . The Hopf boundary correspondsto Eq. (7). SN/SNIC and Focus/Node boundaries are obtained inparametric and explicit form, respectively, in Appendix C . The ho-moclinic boundary has been obtained numerically. The symbol × indicates the parameter value considered in Fig. 4(a). B. Networks with both chemical and electrical coupling
We finally analyze the dynamics of a population ofQIF neurons with mixed, chemical and electrical synapses.Fig. 7(a) presents the possible dynamical regimes of a pop-ulation with chemical synapses only, g = 0 . In contrast tonetworks with pure electrical coupling, where the bifurcationscenario is determined by the presence of three codimension-2points, cf. Fig. 6, here there is only a Cusp point, see also [34].This entails the presence of a persistent state (focus) coexist-ing with an asynchronous, LAS (node) within the cusp-shapedregion in the top-left corner of Fig. 7(a). Additionally, thedashed line indicates that the asynchronous state is of focustype in a vast region of parameters for excitatory coupling,and always for inhibitory coupling.Including electrical coupling, g > , yields the Hopf bi-furcation given by Eq. (7), which joins onto the lower branchof the SN bifurcation curve at a TB point, see Figs. 7(b,c).Hence, the bifurcation scenario for networks with electricaland chemical synapses matches that for networks with elec-trical synapses only: Similar to Fig. 6, the Hopf line cutsthrough the cusp-shaped region —the TB bifurcation demar-cates the point where the Hopf boundary and the lower SN -606 -6 0 6 LAS +Persistent State
FocusNode (a) J / ( π √ ∆ ) -606 -6 0 6(b) Sync TB NodeFocus -606 -6 0 6(c) Sync TB N o d e F o c u s J / ( π √ ∆ ) ¯ η/ ∆ -606 -6 0 6(d) Sync N o d e F o c u s ¯ η/ ∆ FIG. 7: The phase diagram for networks with only chemical cou-pling, panel (a), is characterized by the presence of a Cusp bifurca-tion point. The inclusion of electrical coupling, panels (b-d), trans-forms the Cusp bifurcation scenario into that of Fig. 6, characterizedby the presence of three codimension-2 bifurcation points. The pan-els show the phase diagrams of the FRE (5) with chemical coupling(
J > : Excitatory, J < : Inhibitory) for (a) g/ √ ∆ = 0 , (b) g/ √ ∆ = 1 , (c) g/ √ ∆ = 3 , (d) g/ √ ∆ = 5 . The Hopf bound-aries (red lines) are straight lines given by Eq. (7). SN/SNIC (blacklines) and Focus/Node (dashed) boundaries are obtained in paramet-ric form in Appendix C . Hopf and SN boundaries meet at a TB bifur-cation point. To lighten the diagrams, Cusp, SNSL and homoclinicbifurcations are not shown. Symbols × and + indicate the parametervalues considered in Fig. 4. line intersect. Then, due to the presence of electrical cou-pling, the persistent state becomes only stable in a small pa-rameter region confined between the Hopf and the lower SNline, see Fig. 7(b). As electrical coupling is increased, theTB point approaches the Cusp bifurcation, which results inan even smaller range of parameters for which the persistentstate is stable. This agrees with numerical results using largenetworks of noisy, conductance-based and QIF neurons, andhas been hypothesized to be a possible reason why electricalsynapses are rarely found between excitatory neurons [10].Returning to the analysis of the FRE (5), we find for lowvalues of g that synchronization emerges predominantly forexcitatory coupling, J > , see Fig. 7(b). As electrical cou-pling is increased, the Sync region extends to the inhibitoryregion, J < , and to larger values of ¯ η —note that in this cou-pling regime the emergence of collective oscillations mainlyoccurs via a SNIC bifurcation for excitation and via a Hopfbifurcation for inhibition, see Fig. 7(c). For even larger elec-trical coupling, the TB point moves further into the inhibitoryregion. That is, for strong electrical coupling the J -coordinateof the TB bifurcation rapidly decreases towards minus infin-ity whereas the other coordinate stays relatively close to the ¯ η = 0 axis [88]. The SNIC bifurcation tilts towards a ver-tical line close to the ¯ η = 0 axis, because strong electricalcoupling coerces all neurons to behave as a single QIF neuronwith common input η = ¯ η . Then, the SNIC bifurcation be-comes the only transition between the two possible dynamicalregimes, asynchrony or synchrony, see Fig. 7(d). Fig. 4 shows how the addition of inhibitory coupling into anetwork with only electrical synapses degrades synchrony —parameters used in Fig. 4 correspond to the symbols shown inFig. 7(c). The presence of inhibition clearly slows down theoscillations, as predicted by Eq. (8).Although Eq. (8) is strictly valid only at the Hopf bifur-cation, it is a good estimate of the frequency of the oscilla-tions of the FRE (5). Fig. 5 depicts the comparison betweenEq. (8) as a function of g (dotted lines), with the actual fre-quencies numerically obtained using the FRE (5) (solid lines)and the QIF network Eq. (1) (symbols). In excitatory net-works, the oscillations already emerge for weak values of g . Incontrast, synchronizing inhibitory networks requires a muchlarger value of g , i.e. inhibition does not promote synchro-nization. Remarkably, only the presence of chemical couplingallows the frequency of the oscillations to deviate from ¯ f , seeEqs. (8,10). Oscillations emerge with f > ¯ f for excitationand with f < ¯ f for inhibition, while they remain f ≈ ¯ f fornetworks with only electrical coupling. As g increases, theeffects of chemical coupling are gradually washed out sincean increasing number of neurons are entrained by electricalcoupling to the most-likely frequency of the uncoupled net-work, ¯ f . This dependence is well described by Eq. (8). Fi-nally, since the level of heterogeneity ∆ degrades synchrony,in Eq. (8) this term favors the deviation of the frequency from ¯ f , and compensates for the homogenizing effect of electricalcoupling. In the limit of identical neurons ∆ → , the effectsof instantaneous chemical coupling on the frequency vanish,since neurons synchronize in-phase and, at the instant of fir-ing, all neurons become refractory. V. CONCLUSION AND DISCUSSION
Firing rate models are very useful tools for investigatingthe dynamics of large networks of spiking neurons that inter-act via excitation and inhibition, see e.g. [21–33]. Remark-ably, using a recently proposed approach to derive exact fir-ing rate equations for networks of excitatory and/or inhibitoryQIF neurons [34, 42], Laing has found that electrical synapsescan also be incorporated in the framework of firing rate mod-els [65]. Yet, the FRE in [65] are not exact, and their mathe-matical form makes the analysis intractable.Here we showed that the FRE corresponding to a networkof QIF neurons with both chemical and electrical synapses canbe exactly obtained, without the need for any approximation.Much in the spirit of firing rate models, the resulting FRE (5)are simple in form and highly amenable to analysis. More-over, in
Appendix B , we demonstrate that relaxing the approx-imation invoked in [65] the FRE derived by Laing simplify toour Eqs. (5).At first glance, the mathematical form of Eqs. (5) alreadyunveils two interesting features of electrical and chemical cou-pling, see also Eq. (6): ( i ) Chemical coupling tends to shiftthe center of the distribution of membrane potentials (givenby v ), while electrical coupling tends to reduce the width ofthe distribution (given by r ), potentially promoting the emer-gence of synchronization; ( ii ) While in the original networkof QIF neurons electrical coupling is mediated by membranepotential differences, at the mean-field level the electrical in-teraction is solely mediated by the mean firing rate r .The mathematical analysis of the FRE (5) unravels howchemical and electrical coupling shape the dynamics of glob-ally coupled populations of QIF neurons with Lorentzian het-erogeneity. Some of our results were already reported inprevious work, and confirm the value and the validity ofthe FRE (5). An important conclusion of our study is thatthe presence of electrical coupling, g (cid:54) = 0 , generally im-plies the appearance of a supercritical Hopf bifurcation, seeEq. (7). This Hopf bifurcation meets a SN bifurcating line ina codimension-2 TB point, causing a drastic reduction of theregion of bistability between low-activity and persistent asyn-chronous states, see Figs. 6 and 7. The Hopf bifurcation desta-bilizes the persistent state producing synchronous oscillations,which then are abolished via a homoclinic bifurcation. Previ-ous studies of networks of excitatory and inhibitory neuronsshowed that synchrony often destroys persistent states [68–72]. Moreover, the generality of the bifurcation scenario ofFigs. 6 and 7 —characterized by three codimension-2 points,TB, Cusp and SNSL—, is confirmed in previous studies ana-lyzing closely related systems [73–76].Networks of spiking neurons with strong excitatory cou-pling display robust persistent states. These states emerge at aCusp bifurcation, see Fig. 7(a). Of particular relevance to ourstudy is the work by Ermentrout [10]. He found that electri-cal coupling tends to synchronize neurons, and that this anihi-lates persistent states via the bifurcation scenario described inFigs. 6 and 7. Persistent activity may underlie important cog-nitive functions such as working memory, and has been sug-gested as a possible reason for the lack of electrical couplingbetween excitatory neurons [10]. According to [10], ‘the mainrole for gap junctions is to encourage synchronization duringrhythmic behavior. Synchrony, because it leads to a sharedrefractory period between neurons can lead to the extinctionof persistent activity’ .The Hopf bifurcation is always supercritical. In [17], Os-tojic et al. analyzed the super- or sub-critical character of theHopf bifurcation in networks of electrically coupled leaky in-tegrate and fire (LIF) neurons. At variance with QIF neurons,LIF neurons do not have spikes and, hence, modeling electri-cal coupling requires an additional parameter [15]. This pa-rameter enables one to adjust the shape of the spikelet elicitedin the postsynaptic cell due to an action potential in the presy-naptic cell. Ostojic et al. [17] found that the Hopf bifurca-tion is supercritical when the spikelets are effectively excita-tory, while inhibitory spikelets lead to subcritical Hopf bifur-cations. For the QIF model, the spikelet elicited in a postsy-naptic cell by the transmission of a presynaptic spike has anet excitatory effect —see Fig. 1(b)—, and hence our resultthat the Hopf bifurcation is supercritical is in agreement withthe results in [17]. Yet, we note that our result also includesnetworks with chemical synapses, and not only networks withelectrical synapses, as in [17].Another important result by Ostojic et al. [17] is that elec-trical coupling can lead to oscillations even in the presenceof strong heterogeneity. Our Eq. (7) is consistent with this. Kopell and Ermentrout [8] also investigated the robustness ofsynchrony against current heterogeneities in networks withboth electrical and inhibitory synapses. They found thata small amount of electrical coupling, added to an alreadysignificant inhibitory coupling, can increase synchronizationmore than a very large increase in the inhibitory coupling. InFig. 4 we show that increasing inhibition reduces the ampli-tude of the oscillations in a network with g (cid:54) = 0 . In addition,Fig. 7 shows that, for a given value of g , increasing inhibitionleads to asynchrony. This level of inhibition increases withelectrical coupling, in line with the results in [8]. Two stud-ies [11, 17] also investigated the frequency of the emergingoscillations in networks with electrical synapses. This fre-quency remains tied to the mean firing rate f i in the network(i.e. near ¯ f ), as our Eq. (8) suggests. In Fig. 5 we confirmthat, in networks with only electrical synapses, the frequencyof the oscillations remains near the most likely f i -value: ¯ f .The result that the Hopf bifurcation is always supercritical,and that the frequency of the emerging oscillations is given by ¯ f evoke the paradigmatic synchronization transition in the Ku-ramoto model [67]. For weakly electrically-coupled networks,we find that the onset of oscillations occurs at the Kuramoto’scritical coupling for synchrony, Eq. (9). When consideringchemical coupling, the frequency of the oscillations deviatesfrom ¯ f , increasing/decreasing for excitatory/inhibitory cou-pling. The intensity of this deviation depends on the ratio ofchemical to electrical coupling, as Eq. (8) suggests. Fig. 5confirms that strong electrical coupling overcomes the effectof excitation/inhibition onto the frequency of the oscillations,approximately as dictated by Eq. (8).Together with the firing rate model derived in [65], theFRE (5) constitute a unique example of a firing rate modelwith both electrical and chemical coupling. The numericalsimulations of the original QIF network Eq. (1) are in agree-ment with the FRE (5) —see Figs. 4, 5—, underlining the va-lidity of the reduction method applied. Interestingly, the fixedpoints of the FRE (5) with chemical synapses ( g = 0 , J (cid:54) = 0 )can be cast in the form of a traditional firing rate model r ∗ = Φ(¯ η + Jτ r ∗ ) , (11)where Φ( x ) = (cid:112) x + √ x + ∆ / ( √ πτ ) is the so-calledtransfer function of the heterogeneous QIF network [48–50].The FRE (5) with electrical synapses ( g (cid:54) = 0 ), however, can-not be written in the form of Eq. (11). Therefore, the linkpointed by Eq. (11) between traditional firing rate models andEqs. (5) is lost when electrical coupling is considered. Acknowledgements
This work was supported by ITN COSMOS funded by theEU Horizon 2020 Research and Innovation programme underthe Marie Skłodowska-Curie Grant Agreement No. 642563.EM acknowledges support by the Spanish Ministry of Econ-omy and Competitiveness under Project No. PSI2016-75688-P. [1] X.-J. Wang, Physiological Reviews , 1195 (2010).[2] J. I. Nagy, A. E. Pereda, and J. E. Rash, Biochimica et Biophys-ica Acta (BBA) - Biomembranes , 102 (2018).[3] B. W. Connors, Developmental Neurobiology , 610 (2017).[4] R. D. Traub, M. A. Whittington, R. Guti´errez, and A. Draguhn,Cell and Tissue Research , 671 (2018).[5] M. V. Bennett and R. Zukin, Neuron , 495 (2004).[6] B. W. Connors and M. A. Long, Annual Review of Neuro-science , 393 (2004).[7] M. Whittington, R. Traub, N. Kopell, B. Ermentrout, andE. Buhl, Int. Journal of Psychophysiol. , 315 (2000), ISSN0167-8760.[8] N. Kopell and B. Ermentrout, Proceedings of the NationalAcademy of Sciences , 15482 (2004).[9] B. Pfeuty, D. Golomb, G. Mato, and D. Hansel, Frontiers inComputational Neuroscience , 8 (2007), ISSN 1662-5188.[10] B. Ermentrout, Phys. Rev. E , 031918 (2006).[11] A. Viriyopase, R.-M. Memmesheimer, and S. Gielen, Journalof Neurophysiology , 232 (2016).[12] A. Holzbecher and R. Kempter, European Journal of Neuro-science , 3446 (2018).[13] D. Guo, Q. Wang, and M. c. v. Perc, Phys. Rev. E , 061905(2012).[14] T. Tchumatchenko and C. Clopath, Nature Communications ,5512 (2014).[15] T. J. Lewis and J. Rinzel, Journal of Computational Neuro-science , 283 (2003).[16] C. C. Chow and N. Kopell, Neural Computation , 1643(2000).[17] S. Ostojic, N. Brunel, and V. Hakim, Journal of ComputationalNeuroscience , 369 (2009), ISSN 1573-6873.[18] B. Pfeuty, G. Mato, D. Golomb, and D. Hansel, Journal of Neu-roscience , 6280 (2003).[19] S. Coombes, SIAM Journal on Applied Dynamical Systems ,1101 (2008).[20] J. G. Mancilla, T. J. Lewis, D. J. Pinto, J. Rinzel, and B. W.Connors, Journal of Neuroscience , 2058 (2007).[21] H. R. Wilson and J. D. Cowan, Biophys. J. , 1 (1972).[22] P. Dayan and L. F. Abbott, Theoretical neuroscience (Cam-bridge, MA: MIT Press, 2001).[23] G. B. Ermentrout and D. H. Terman,
Mathematical foundationsof neuroscience , vol. 64 (Springer, 2010).[24] J. J. Hopfield, Proceedings of the national academy of sciences , 3088 (1984).[25] G. Mongillo, O. Barak, and M. Tsodyks, Science , 1543(2008).[26] R. Ben-Yishai, R. L. Bar-Or, and H. Sompolinsky, Proc. Nat.Acad. Sci. , 3844 (1995).[27] D. Hansel and H. Sompolinsky, in Methods in Neuronal Mod-elling: From Ions to Networks , edited by C. Koch and I. Segev(MIT Press, Cambridge, 1998), pp. 499–567.[28] J. Tabak, W. Senn, M. J. ODonovan, and J. Rinzel, J. Neurosci. , 3041 (2000).[29] J. Rankin, E. Sussman, and J. Rinzel, PLoS Computational Bi-ology , 1 (2015).[30] D. Mart´ı and J. Rinzel, Neural computation , 1 (2013).[31] M. H. Tsodyks M., Pawelzik K., Neural Comput. , 821(1998).[32] A. Roxin, N. Brunel, and D. Hansel, Phys. Rev. Lett. ,238103 (2005).[33] A. Roxin and E. Montbri´o, Physica D , 323 (2011). [34] E. Montbri´o, D. Paz´o, and A. Roxin, Phys. Rev. X , 021028(2015).[35] E. Ott and T. M. Antonsen, Chaos , 037113 (2008).[36] E. Ott and T. M. Antonsen, Chaos , 023117 (2009).[37] E. Ott, B. R. Hunt, and T. M. Antonsen, Chaos , 025112(2011).[38] A. Pikovsky and M. Rosenblum, Phys. Rev. Lett. , 264103(2008).[39] S. A. Marvel, R. E. Mirollo, and S. H. Strogatz, Chaos: An In-terdisciplinary Journal of Nonlinear Science , 043104 (2009).[40] A. Pikovsky and M. Rosenblum, Physica D , 872 (2011).[41] B. Pietras and A. Daffertshofer, Chaos: An InterdisciplinaryJournal of Nonlinear Science , 103101 (2016).[42] T. B. Luke, E. Barreto, and P. So, Neural Comput. , 3207(2013).[43] P. So, T. B. Luke, and E. Barreto, Physica D , 16 (2014).[44] C. R. Laing, Phys. Rev. E , 010901 (2014).[45] D. Paz´o and E. Montbri´o, Phys. Rev. Lett. , 238101 (2016).[46] I. Ratas and K. Pyragas, Phys. Rev. E , 032215 (2016).[47] J. Roulet and G. B. Mindlin, Chaos: An Interdisciplinary Jour-nal of Nonlinear Science , 093104 (2016).[48] F. Devalle, E. Montbri´o, and D. Paz´o, Phys. Rev. E , 042214(2018).[49] F. Devalle, A. Roxin, and E. Montbri´o, PLoS ComputationalBiology (2017).[50] J. M. Esnaola-Acebes, A. Roxin, D. Avitabile, and E. Montbri´o,Phys. Rev. E , 052407 (2017).[51] I. Ratas and K. Pyragas, Phys. Rev. E , 042212 (2017).[52] G. Dumont, G. B. Ermentrout, and B. Gutkin, Phys. Rev. E ,042311 (2017).[53] ´A. Byrne, M. J. Brookes, and S. Coombes, Journal of Compu-tational Neuroscience , 143 (2017).[54] C. R. Laing, The Journal of Mathematical Neuroscience , 4(2018), ISSN 2190-8567.[55] H. Schmidt, D. Avitabile, E. Montbri´o, and A. Roxin, PLoSComputational Biology (2018).[56] I. Ratas and K. Pyragas, Phys. Rev. E , 052224 (2018).[57] M. di Volo and A. Torcini, Phys. Rev. Lett. , 128301 (2018).[58] G. Dumont and B. Gutkin, PLOS Computational Biology , 1(2019).[59] A. Akao, S. Shirasaka, Y. Jimbo, B. Ermentrout, and K. Kotani,arXiv preprint arXiv:1903.12155 (2019).[60] H. Bi, M. Segneri, M. di Volo, and A. Torcini, bioRxiv (2019).[61] 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.[62] S. Keeley, A. Byrne, A. Fenton, and J. Rinzel, Journal of Neu-rophysiology , 2181 (2019).[63] A. Byrne, D. Avitabile, and S. Coombes, Phys. Rev. E ,012313 (2019).[64] S. Boari, G. Uribarri, A. Amador, and G. B. Mindlin, Mathe-matical and Computational Applications (2019).[65] C. R. Laing, SIAM Journal on Applied Dynamical Systems ,1899 (2015).[66] D. Paz´o and E. Montbri´o, Phys. Rev. E , 055202 (2006).[67] Y. Kuramoto, in International Symposium on MathematicalProblems in Theoretical Physics , edited by H. Araki (Springer,Berlin, 1975), vol. 39 of
Lecture Notes in Physics , pp. 420–422.[68] D. Hansel and G. Mato, Phys. Rev. Lett. , 4175 (2001).[69] D. Hansel and G. Mato, Neural Computation , 1 (2003).[70] B. S. Gutkin, C. R. Laing, C. L. Colby, C. C. Chow, and G. B. Ermentrout, Journal of Computational Neuroscience , 121(2001).[71] C. R. Laing and C. C. Chow, Neural Computation , 1473(2001).[72] T. Kanamaru and M. Sekine, Phys. Rev. E , 031916 (2003).[73] H. Sakaguchi, S. Shinomoto, and Y. Kuramoto, Progress ofTheoretical Physics , 600 (1988).[74] M. A. Zaks, A. B. Neiman, S. Feistel, and L. Schimansky-Geier,Phys. Rev. E , 066206 (2003).[75] L. M.Childs and S. H. Strogatz, Chaos , 043128 (2008).[76] L. F. Lafuerza, P. Colet, and R. Toral, Phys. Rev. Lett. ,084101 (2010).[77] H. Daido and K. Nakanishi, Phys. Rev. Lett. , 104101 (2004).[78] Y. Kuramoto, Chemical Oscillations, Waves, and Turbulence (Springer-Verlag, Berlin, 1984).[79] In our numerical simulations (Euler scheme, δt = 10 − ), theresetting rule was applied as follows: When V j ≥ , themembrane voltage is held at V j for a time interval τ /V j . Then,a spike is emitted, and the voltage is reset and kept at − V j for asubsequent interval τ /V j . In Figs. 4 and 5, to numerically eval-uate the mean membrane potential v , the population average iscomputed discarding those neurons with | V j | ≥ .[80] We used τ s = 10 − ms, and τ = 10 ms. The instantaneousfiring rates shown in Fig. 4 were obtained by binning time andcounting spikes within a sliding time window of size δt = 2 . × − ms.[81] This phenomenon has been termed ‘Aging transition’ in the lit-erature [66, 77].[82] Note that these spikelets depend on the shape of the presynapticspike, and thus on the particular neuron model considered.[83] The value ¯ η = 0 determines a boundary separating networkoscillations and quiescence, cf. Eq (6) with p ∞ c = 1 / in [66].[84] See Appendix B for the comparison of the firing rate modelEqs.(5) with the FRE derived in [65].[85] This might be understood as follows: The evolution equationfor the mean membrane potential v is obtained summing upthe N differential equations Eq. (1): N (cid:80) Ni =1 ˙ V i = ˙ v = N (cid:80) Ni =1 ( V i + η i ) + Js + g ( v − N (cid:80) Ni =1 V i ) . Although thisis not a closed equation for v and r , one finds that the last term—corresponding to diffusive coupling— cancels to zero.[86] For non-instantaneous inhibitory synapses, oscillations emerg-ing through a Hopf bifurcation are also encountered for weakcoupling and weak heterogeneity [48, 49]. These oscillationsare often referred to as ‘Interneuronal Gamma Oscillations’,ING [7]. To keep our analysis simple, here we do not considerING oscillations, since Eqs. (5) become higher-dimensional andphase plane analysis is no longer possible in this case.[87] Exactly the same scenario is found in systems of globally, sinecoupled ‘active rotators’. Such systems are, in fact, closely re-lated to the dynamics (1) with J = 0 [73–76].[88] It can be shown that the TB point behaves proportional to (˜ η, ˜ J ) TB ∝ ( g , − g ) in the limit g → ∞ . Appendix A: Derivation of the Firing Rate Equations
In the thermodynamic limit, N → ∞ , we drop the in-dices for the individual neuronal dynamics Eq. (1), and de-note ρ ( V | η, t ) dV as the fraction of neurons with membranepotentials between V and V + dV , and parameter η at time t .Accordingly, the parameter η becomes a continuous randomvariable that is distributed according to a probability distri- bution function, which here is considered to be a Lorentzian L ∆ , ¯ η ( η ) of half-width ∆ and centered at ¯ η , see Eq. (4). Theconservation of the number of neurons leads to the continuityequation τ ∂ t ρ + ∂ V (cid:2) ( V + η + g ( v − V ) + Jτ r ) ρ (cid:3) = 0 , (A1)where we explicitly included the velocity given by the contin-uous equivalent of Eq. (1). We also defined the mean value ofthe membrane potential as v ( t ) = (cid:90) ∞−∞ (cid:90) ∞−∞ ρ ( V | η, t ) V L ∆ , ¯ η ( η ) dV dη. (A2)Next, we consider the family of conditional density func-tions [34] ρ ( V | η, t ) = 1 π x ( η, t )[ V − y ( η, t )] + x ( η, t ) , (A3)which are Lorentzian functions with time-dependent half-width x ( η, t ) , centered at y ( η, t ) . Substituting (A3) into thecontinuity equation (A1), we find that, for each value of η ,variables x and y must obey two coupled equations, τ ˙ x ( η, t ) = 2 x ( η, t ) y ( η, t ) − gx ( η, t ) , (A4a) τ ˙ y ( η, t ) = η − x ( η, t ) + y ( η, t ) (A4b) − g (cid:2) y ( η, t ) − v (cid:3) + Jτ r, that can be written in complex form as τ ∂ t w ( η, t ) = i (cid:2) η − w ( η, t ) + Jτ r (cid:3) + g (cid:2) iv − w ( η, t ) (cid:3) (A5)where w ( η, t ) ≡ x ( η, t ) + iy ( η, t ) . For a particular value of η , the firing rate r of the population of QIF neurons is relatedto the width x of the Lorentzian ansatz (A3). Specifically, thefiring rate r ( η, t ) for each η value at time t is the probabilityflux at infinity: r ( η, t ) = ρ ( V → ∞| η, t ) ˙ V ( V → ∞| η, t ) ,which yields the identity x ( η, t ) = πτ r ( η, t ) . (A6)Hence, integrating this quantity over the distributions of cur-rents Eq. (4) provides the mean firing rate r ( t ) = 1 τ π (cid:90) ∞−∞ x ( η, t ) L ∆ , ¯ η ( η ) dη. (A7)Likewise, we can link the center y ( η, t ) of the Lorentzianansatz Eq. (A3) with the mean of the (conditional) membranepotential via y ( η, t ) = p.v. (cid:90) ∞−∞ ρ ( V | η, t ) V dV. (A8)Note that the Lorentzian distribution does not have finite mo-ments so that the integral in Eq. (A8) needs to be takenas the Cauchy principal value (i.e. p.v. (cid:82) ∞−∞ ρV dV =lim R →∞ (cid:82) R − R ρV dV ). Then, Eq. (A2) becomes v ( t ) = (cid:90) ∞−∞ y ( η, t ) L ∆ , ¯ η ( η ) dη. (A9)0The integrals in (A7,A9) can be evaluated closing the in-tegral contour in the complex η -plane and using Cauchy’sresidue theorem. The integrals must however be performedcarefully, so that the variable x ( η, t ) remains non-negative.To make the analytic continuation of w ( η, t ) from real tocomplex-valued η , we define η ≡ η r + iη i . This contin-uation is possible into the lower half-plane η i < , sincethis guarantees the half-width x ( η, t ) to remain non-negative: ∂ t x ( η, t ) = − η i > , at x = 0 . Therefore, we perform con-tour integration in Eq. (A7) and Eq. (A9) along the arc | η | e iϑ with | η | → ∞ and ϑ ∈ ( − π, . This contour encloses onepole of the Lorentzian distribution Eq. (4). Then, we find thatthe firing rate and the mean membrane potential depend onlyon the value of w at the pole of L ∆ , ¯ η ( η ) in the lower half η -plane: πτ r ( t ) + iv ( t ) = w (¯ η − i ∆ , t ) , As a result, we only need to evaluate Eq. (A5) at η = ¯ η − i ∆ , and obtain a system of FRE composed of two ordinarydifferential equations as given in Eq. (5), τ ˙ r = ∆ τ π + 2 rv − gr,τ ˙ v = v + ¯ η − ( πτ r ) + Jτ r, in terms of the population-mean firing rate r and thepopulation-mean membrane potential v . Multiplying theLorentzian ansatz Eq. (A3) by L ∆ , ¯ η ( η ) and integrating over η , we finally obtain the total density of neurons Eq. (6) as ρ ( V, t ) = 1 π πτ r ( t )[ V − v ( t )] + π τ r ( t ) , where we again applied Cauchy’s residue theorem by usingthat the ansatz Eq. (A3) is analytic in the lower η -complexplane. Hence, the total density of the population of QIF neu-rons is a Lorentzian distribution centered at v ( t ) and half-width πτ r ( t ) , which evolves according to the FRE (5). Appendix B: Connection between the FRE in [65] and Eq. (5)
The derivation of the FRE (5) is exact in the thermodynamiclimit, and does not rely on any approximation. Here we showthat the Eqs. (2.35&2.36) in [65] reduce to our Eq. (5) afteradopting a limit in which the derivation performed in [65] be-comes exact.In contrast to our Eq. (5b), note that Eq. (2.36) in [65] con-tains a diffusive term, g [ Q ( t ) − v ( t )] , (B1)where the function Q ( t ) is defined as Q ( t ) = i ∞ (cid:88) m =1 ρ m +1 − ρ m − ρ + 1 + (cid:15) [ z m − ¯ z m ] (B2)with < (cid:15) (cid:28) , and ρ = √ (cid:15) + (cid:15) − − (cid:15) . The variable z in Eq. (B2) is the complex Kuramoto order parameter (the bar denotes complex conjugation), which is related to the vari-ables r and v in the FRE (5) via the change of variables [34] πr + iv = 1 − ¯ z z . (B3)The parameter (cid:15) , defined in Eq. (2.7) in [65], was used to ap-proximate the mean voltage v , see also [10]. In the limit (cid:15) → this approximation becomes exact, but this limit was not con-sidered in [65]. In consequence, to use the Eqs. (2.35&2.36)in [65], the infinite series Eq. (B2) was truncated after 100terms, and the bifurcation analysis of the mean-field modelcould only be performed numerically.Using the geometric series formula ( | z | < ) and the trans-formation of variables Eq. (B3), we find lim (cid:15) → Q ( t ) = i ∞ (cid:88) m =1 ( − m [ z m − ¯ z m ] = 2 Im ( z )(1 + z )(1 + ¯ z ) = v. Hence we have showed that the diffusive term Eq. (B1) iden-tically vanishes when the mean-field reduction becomes exact(i.e. in the limit (cid:15) → ), and the FRE in [65] reduce to Eqs. (5). Appendix C: Bifurcation analysis of the Firing Rate Equations
The FRE (5) have five free parameters. The number ofeffective parameters can be reduced to three through non-dimensionalization, defining ˜ η = ¯ η/ ∆ , ˜ g = g/ √ ∆ , ˜ J = J/ ( π √ ∆) , and rescaling variables as ˜ r = τ πr/ √ ∆ , ˜ v = v/ √ ∆ , ˜ t = √ ∆ t/τ. Then, the firing rate model becomes d ˜ rd ˜ t = 1 + 2˜ r ˜ v − ˜ g ˜ r, (C1a) d ˜ vd ˜ t = ˜ v + ˜ η − ˜ r + ˜ J ˜ r, (C1b)The fixed points (˜ r ∗ , ˜ v ∗ ) of Eq. (C1) satisfy ˜ v ∗ = ˜ g − r ∗ . (C2)Linearization about the fixed points Eq. (C2) gives the eigen-values λ ± = 12 (cid:18) v ∗ − ˜ g ± (cid:113) ˜ g + 8˜ r ∗ ( ˜ J − r ∗ ) (cid:19) . (C3)For networks with only chemical synapses (i.e. g = 0 ), thereal part of the eigenvalues remains always negative (since v ∗ < ), and a Hopf bifurcation is not possible. However,chemical coupling has a direct influence on the real part of theeigenvalues Eq. (C3), and may produce oscillatory instabili-ties if the argument of the square root is a real number.1 A. Hopf boundaries and Takens-Bogdanov point
The Hopf boundaries can be obtained when imposingRe( λ ± ) = 0 in Eq. (C3), which gives ˜ g = 4˜ v ∗ . Then, us-ing Eq. (C2), we find ˜ g H = 2 / ˜ r ∗ . (C4)Substituting Eq. (C4) in the v -fixed point equation Eq. (C1b),and solving for ˜ η , we obtain ˜ η H = r ∗ − ˜ J ˜ r ∗ − r ∗ . (C5)Solving Eq. (C4) for ˜ r ∗ and substituting it into Eq. (C5) weobtain the Hopf boundaries in explicit form ˜ η H = − J ˜ g + 4˜ g − ˜ g . (C6)The frequency of the oscillations is given by the imaginarypart of the eigenvalues Eq. (C3) at criticality that, using thefixed points of Eqs. (C1) and Eq. (C4), reduces to the explicitformula ˜ f H = 1 π (cid:115) ˜ η + ˜ J ˜ g . (C7)The frequency becomes zero at a Takens-Bogdanov (TB)point, when ˜ η = − ˜ J/ ˜ g . Inserting this condition into Eq. (C6)we obtain the coordinates of the TB point (cid:16) ˜ η, ˜ J (cid:17) TB = (cid:18) ˜ g − g , g − ˜ g (cid:19) , (C8)see also Fig. 7. For ˜ J = 0 , the TB point is located at (˜ η, ˜ g ) TB = (cid:16) , √ (cid:17) (C9)in the phase diagram Fig. 6. B. Saddle-node boundaries
The boundaries of the saddle-node bifurcations are ob-tained by setting λ ± = 0 in Eq. (C3), using Eq. (C2), andsolving for ˜ g : ˜ g sn = 1˜ r ∗ − J ˜ r ∗ + 4˜ r ∗ . (C10)Substituting (C10) in the v -fixed point Eq. (C1b), and solvingfor ˜ η , we obtain ˜ η sn = ˜ r ∗ − r ∗ + ˜ J ˜ r ∗ (cid:0) r ∗ − J ˜ r ∗ − (cid:1) . (C11)The saddle node boundaries are plotted in the (˜ η, ˜ g ) phase di-agram in Fig. 6. The same boundaries can be represented inthe (˜ η, ˜ J ) phase diagram when solving Eq. (C10) for ˜ J ˜ J sn = 12˜ r ∗ − ˜ g r ∗ + 2˜ r ∗ (C12) and replacing ˜ J by Eq. (C12) in Eq. (C11) gives ˜ η sn = ˜ g ˜ r ∗ − ˜ g − r ∗ − r ∗ . These saddle-node boundaries are shown in black for differentvalues of ˜ g in Fig. 7. C. Focus-Node boundaries
The boundaries in the phase diagram Fig. 7 in which thestable asynchronous state changes from Focus to Node canbe obtained in parametric form equating the square root inEq. (C3) to zero. This gives ˜ J F N = 16˜ r ∗ − ˜ g r ∗ . Substituting ˜ J F N into the v -fixed point Eq. (C1b), and usingEq. (C2) we find ˜ η F N = 18˜ r ∗ (cid:0) − g ˜ r ∗ − ˜ g ˜ r ∗ − r ∗ (cid:1) . For networks without chemical coupling, ˜ J = 0 , the Focus-Node boundary can be obtained in explicit form ˜ η F N = 2 − g − g . This is the dashed boundary depicted in Fig. 6.
Appendix D: Small-amplitude equation near the Hopfbifurcation
In this Appendix we derive the small amplitude equationnear the Hopf bifurcation, and show that the Hopf bifurca-tion is always supercritical. The derivation is performed usingmultiple-scales analysis, see e.g. [78]. We first expand the so-lution (cid:18) ˜ r ˜ v (cid:19) = (cid:18) r v (cid:19) + (cid:15) (cid:18) r v (cid:19) + (cid:15) (cid:18) r v (cid:19) + . . . (D1)in powers of a small parameter (cid:15) (cid:28) , about a fixed point ( r , v ) of Eqs. (C1) at the Hopf bifurcation. In addition, weintroduce the deviation from the Hopf bifurcation Eq. (C6) ofparameter ˜ η as ˜ η − ˜ η H = χ(cid:15) , (D2)where χ determines the sign of the deviation. Finally, we de-fine the slow time T = (cid:15) t. (D3)Then, the time differentiation is transformed as ddt → ∂ t + (cid:15) ∂ T . (D4)2Plugging Eqs. (D1,D2,D4) into Eq. (C1) gives (cid:2) ∂ t + (cid:15) ∂ T − L (cid:3) (cid:20) (cid:15) (cid:18) r v (cid:19) + (cid:15) (cid:18) r v (cid:19) + . . . (cid:21) − (cid:15) (cid:18) χ (cid:19) == (cid:15) N + (cid:15) N + . . . , (D5)where L = (cid:18) v − ˜ g r ˜ J − r v (cid:19) . (D6)and N = (cid:18) r v v − r (cid:19) , N = (cid:18) r v + 2 r v v v − r r (cid:19) . (D7) Critical eigenvectors
Using Eqs. (C3,C4), we define the critical frequency as ω = 12˜ g (cid:113) − ˜ g −
16 ˜ J ˜ g, (D8)so that the matrix Eq. (D6) can be written as L = (cid:18) − ˜ g/ / ˜ g − ˜ g/ g + 4 ω ) ˜ g/ (cid:19) . (D9)The right-eigenvector of Eq. (D9) is u R = (cid:18) / ˜ g ˜ g/ iω (cid:19) . (D10)Imposing the condition u L u R = 1 , the left-eigenvector of L is u L = 12 ω (cid:18) gω + i ˜ g , − i (cid:19) . (D11) Analysis of multiple scales
At order (cid:15) , Eq. (D5) is (cid:18) ˙ r ˙ v (cid:19) − L (cid:18) r v (cid:19) = (cid:18) (cid:19) . (D12) This system of differential equations has a general solution (cid:18) r v (cid:19) = Ae iω ˜ t u R + c.c , (D13)which is the so-called neutral solution.At order (cid:15) , Eq. (D5) is (cid:18) ˙ r ˙ v (cid:19) − L (cid:18) r v (cid:19) − (cid:18) χ (cid:19) = N . (D14)Substituting the neutral solution Eq. (D13) into Eq. (D7) wefind N = (cid:18) g + 4˜ g ω − / (2˜ g ) (cid:19) | A | + (cid:18) iω / ˜ g ˜ g / i ˜ gω − ω − / ˜ g (cid:19) A e iω ˜ t + c.c.Next we use the following ansatz (cid:18) r v (cid:19) = (cid:18) r v (cid:19) + (cid:18) r v (cid:19) e iω ˜ t + c.c. (D15)and substitute it into Eq. (D14). We find r = 2˜ g ω (cid:2) g χ − (64 + ˜ g − g ω ) | A | (cid:3) ,v = 14˜ gω (cid:2) g χ − (64 + ˜ g + 4˜ g ω ) | A | (cid:3) ,r = A g ω (cid:2)