Bumps and Oscillons in Networks of Spiking Neurons
BBumps and Oscillons in Spiking Networks
Bumps and Oscillons in Networks of Spiking Neurons
Helmut Schmidt a) and Daniele Avitabile Max Planck Institute for Human Cognitive and Brain Sciences, Stephanstrasse 1a, 04103 Leipzig,Germany Department of Mathematics, Vrije Universiteit (VU University Amsterdam), Faculteit der Exacte Wetenschappen,De Boelelaan 1081a, 1081 HV Amsterdam,Mathneuro Team, Inria Sophia Antipolis, 2004 Rue des Lucioles, 06902 Sophia Antipolis,Cedex. (Dated: 23 March 2020)
We study localized patterns in an exact mean-field description of a spatially-extended network of quadraticintegrate-and-fire (QIF) neurons. We investigate conditions for the existence and stability of localized so-lutions, so-called bumps , and give an analytic estimate for the parameter range where these solutions existin parameter space, when one or more microscopic network parameters are varied. We develop Galerkinmethods for the model equations, which enable numerical bifurcation analysis of stationary and time-periodicspatially-extended solutions. We study the emergence of patterns composed of multiple bumps, which arearranged in a snake-and-ladder bifurcation structure if a homogeneous or heterogeneous synaptic kernel issuitably chosen. Furthermore, we examine time-periodic, spatially-localized solutions ( oscillons ) in the pres-ence of external forcing, and in autonomous, recurrently coupled excitatory and inhibitory networks. In bothcases we observe period doubling cascades leading to chaotic oscillations.
Spatially extended networks of spiking modelneurons are capable of producing spatio-temporalpatterns that are observed experimentally in neu-ronal tissues. An important tool in investigat-ing such patterns are low-dimensional neural fieldmodels, which describe the macroscopic dynam-ics of such networks. Many neural field modelsare derived heuristically, and do not fully de-scribe the dynamics of the underlying networkof spiking neurons. We utilize a recently derivedmean-field description for networks of quadraticintegrate-and-fire neurons, which yields an accu-rate description of the mean firing rate and themean membrane potential of the network. Con-trary to other neural field models, this model onlycontains nonlinearities of the mean field variablesup to quadratic order, which are amenable to thedevelopment of Galerkin methods for the numeri-cal approximation of the problem. This allows usto study the bifurcation structure of stationarylocalized solutions, and of time-varying localizedsolutions up to the emergence of chaos.
I. INTRODUCTION
Localized states in neuronal networks, so-calledbumps, are related to working memory and featureselectivity , whereby neurons encoding similar stimuli orfeatures show an increased firing rate for the duration ofthe related cognitive task. Neural fields are well-known a) [email protected]; corresponding author coarse-grained models of spatio-temporal neuronal activ-ity , capable of reproducing dynamic phenomena foundexperimentally, such as traveling waves, temporal oscilla-tions, and spatially localized states . A challenge facedin the derivation of neural field models is to establish anaccurate mean-field description of the spiking dynamicsof the underlying microscopic neural network. Classi-cal neural field models recover the microscopic dynamicsonly in the limit of slow synapses , and the derivationof neural mass or neural field models from networks ofspiking neurons is still an active area of research . Inaddition, in neural fields the network firing rate is not anemergent quantity, but rather the result of a modellingchoice.Some limitations can be overcome if the microscopicmodel is a heterogeneous network of synaptically coupled θ or QIF neurons, subject to random, Cauchy-distributedbackground currents. Recently, it has been shown thatheterogeneous networks of θ - and QIF neurons admit anexact mean field description , which has later been ex-panded to spatially-extended networks . In the ther-modynamic limit, the network admits an exact mean fielddescription in terms of the network mean firing rate andvoltage , or in terms of a complex-valued order param-eter .We study a network of n quadratic integrate-and-fireneurons: ˙ V i = V i + η i + Js i , i = 1 , . . . , n, (1)where V i is the membrane potential of the i th neu-ron, η i an intrinsic current, s i the synaptic input and J > a global coupling parameter. The i th neu-ron emits a spike when V i reaches the firing thresh-old V θ , and V i is reset immediately to V r . FollowingReference , we distribute { η i : i = 1 , . . . , n } accord-ing to a Lorentzian distribution using the formula η i = η +∆ tan [ π/ j − n − / ( n + 1)] , where η is the center a r X i v : . [ n li n . PS ] M a r umps and Oscillons in Spiking Networks 2and ∆ is the half-width of the Lorentzian distribution, re-spectively. An important difference in the model consid-ered in the present paper is that neurons are distributedin space, in a domain Ω = ( − L/ , L/ , with L (cid:29) , atevenly spaced positions { x i = iδx − L/ i ∈ , . . . , n } with δx = L/n . We associate with each lattice point x i a random component of the vector { η j } , without rep-etitions. The synaptic current received by a neuron isdetermined by the synaptic footprint as follows s i ( t ) = 1 n n (cid:88) j =1 w ( x i , x j ) (cid:88) k : t kj Stationary states of Equation (3) are determined bythe conditions ∂ t r = ∂ t v = 0 . Bounded solutions with r ( x ) > satisfy π r + η + Jw ⊗ r − π r , v = − ∆2 πr . (6)The model supports both uniform and non-uniformsteady states, which we discuss below in further detail. A. Spatially uniform states Solutions to (6) depend in general on x . Spatially-uniform solutions, for which Jw ⊗ r = Jr , satisfy thequartic equation r − Jπ r − ηπ r − ∆ π = 0 , (7)which has the following four solutions, r , = J π + 12 √ S ± (cid:113) − p − S − q/ √ S,r , = J π − √ S ± (cid:113) − p − S + 2 q/ √ S, (8)where p , q , and S are given by p = − ηπ − J π , q = − J π − Jηπ , (9) S = − p + 13 ( Q + R /Q ) , (10)with Q = (cid:18) (cid:18) R + (cid:113) R − R (cid:19)(cid:19) / ,R = η π − π ,R = − η π − J ∆ π − η ∆ π , (11)respectively. Physically-relevant solutions are positiveand real, and an inspection of the equations above re-veals that r must be discarded, and the system admitseither or homogeneous steady states. At sufficientlysmall (large) η only one stable fixed point exists, repre-sented by r ( r ); also, there exists an interval in param-eter space where the stable solutions r , r coexist with r , which is unstable. The conclusions presented abovejustify the bifurcation diagram found in References ,and reported in Figure 2a.Loci of saddle-node bifurcations in the ( η, J ) -plane canbe found by setting d η/ d r = 0 in the first equation in (6) which, combined with (7) yields a parameterization in rη sn = − π r − π r ,J sn = 2 π r + ∆ π r , (12)or, more explicitly, J sn = √ π (cid:113) − η sn ± (cid:112) η sn − + √ π ∆ (cid:16) − η sn ± (cid:112) η sn − (cid:17) / , (13)where ± denote two bifurcation branches of saddle-nodebifurcation which collide at a cusp ( η c , J c ) = (cid:18) −√ , π (cid:113) √ (cid:19) . (14) B. Turing bifurcations A first step towards the construction of heterogeneoussteady states is the determination of Turing bifurcations,which mark points in parameter space where a spatiallyuniform solution becomes unstable to spatially periodicpatterns. We remark that it is known that spatially-extended networks of QIF or θ neurons display this in-stability , and here we present an analytic determi-nation of the loci of such bifurcation in parameter space.Turing bifurcations of a homogeneous steady state ( r, v ) can be identified by linear stability analysis of the modelequations in Fourier space, which results in the followingeigenvalue problem: λ ( k ) (cid:18) ˜ r ˜ v (cid:19) = (cid:18) v rJ ˆ w ( k ) − π r v (cid:19) (cid:18) ˜ r ˜ v (cid:19) := ˆ A ( k ) (cid:18) ˜ r ˜ v (cid:19) , where ˆ w ( k ) is the Fourier transform of the connectivitykernel. A sufficient condition for a Turing bifurcation isthe existence of a critical wavenumber k c > for which det A ( k c ) = 0 , which yields r T = 1 π (cid:118)(cid:117)(cid:117)(cid:116) − η T ˆ w ( k c )2(2 − ˆ w ( k c ) ± (cid:115) η T ˆ w ( k c ) (2 − ˆ w ( k c )) − (2+ ˆ w ( k c ))∆ − ˆ w ( k c ) , (15)and J T = 1ˆ w ( k c ) (cid:18) ∆ π r T + 2 π r T (cid:19) . (16)Combining (15) and (16) results in an equation for theloci of the Turing bifurcation in the ( η, J ) -plane. As k c → the resulting equation recovers (13), since ˆ w ( k c ) → . This analytic result agrees well with the numericalcalculations of these loci, which will be presented furtherbelow.umps and Oscillons in Spiking Networks 4 C. Spatial dynamical system After studying uniform and spatially-periodic steadystates, we move to the construction of localized steadystates supported by the models. One strategy to study lo-calized stationary states in nonlinear models posed on R is to construct solutions to boundary-value problems de-rived from the model’s steady state equations . Withthis approach, localized steady states correspond to ho-moclinic orbits of a dynamical system in which x playsthe role of time (hence the term spatial dynamics ).In this section we make some preliminary considera-tions on the spatial dynamics of steady state solutionsto (3), although we do not explicitly study the associ-ated spatial-dynamical system, as we will construct oursolutions numerically in the following sections. Using thepositions u ( r ) = − ∆ π r − η + π r = − v ( r ) − η + π r , (17) f ( u ( r )) = Jr, (18)the steady state equation (6) is recast as − u + w ⊗ f ( u ) . (19)We note that Eq. (19) is formally equivalent to the Amaristeady state equation ; in Amari’s theory u representsthe voltage, whereas in this case u combines the steadystate’s voltage and rate, and scales as u ∼ − v − η forsmall r and u ∼ π r − η for large r , respectively.Importantly, the identification with the Amari equa-tion allows us to use spatial dynamics to characterizelocalized steady state solutions . The Fourier trans-form of w is of the form ˆ w ( k ) = Π ( k ) / Π ( k ) , with Π i being a polynomial of order i , if w is the biexponentialkernel (5). Hence, the integral kernel can be regarded asGreen’s function of a differential operator. In particu-lar, the bi-exponential kernel (5) leads to the differentialequation u (cid:48)(cid:48)(cid:48)(cid:48) − u (cid:48)(cid:48) + 14 u − f ( u ) + 74 [ f ( u )] (cid:48)(cid:48) = 0 , (20)where a prime denotes differentiation with respect to x .The equation above can be cast as a 4D, first-order spa-tial dynamical system in the vector ( u, u (cid:48) , u (cid:48)(cid:48) , u (cid:48)(cid:48)(cid:48) ) , whichwe omit here for brevity. To construct localized solutionsto (3) we proceed in the same spirit as : to eachhomogeneous steady state of (3) corresponds one value r j in (8), and hence one value u j in (17), and one con-stant solution ( u j , , , to (20); in addition, there existsa region in parameter space where r and r coexist andare stable (see also Figure 2a). A localized steady stateof (3) is identified with a bounded, sufficiently regularfunction u : R → R which satisfies (20) with boundaryconditions lim x →−∞ (cid:0) u ( x ) , u (cid:48) ( x ) , u (cid:48)(cid:48) ( x ) , u (cid:48)(cid:48)(cid:48) ( x ) (cid:1) = ( u , , , , lim x → + ∞ (cid:0) u ( x ) , u (cid:48) ( x ) , u (cid:48)(cid:48) ( x ) , u (cid:48)(cid:48)(cid:48) ( x ) (cid:1) = ( u , , , . Furthermore, we note that the quantity H ( u, u (cid:48) , u (cid:48)(cid:48) , u (cid:48)(cid:48)(cid:48) , x ) = u (cid:48)(cid:48)(cid:48) u (cid:48) − 12 ( u (cid:48)(cid:48) ) − 58 ( u (cid:48) ) + 18 u − (cid:90) u f ( z ) d z + 74 (cid:90) x [ f ( u )] (cid:48)(cid:48) ( z ) u (cid:48) ( z ) d z, is conserved in the sense that, if (20) holds, thendd x H ( u ( x ) , u (cid:48) ( x ) , u (cid:48)(cid:48) ( x ) , u (cid:48)(cid:48)(cid:48) ( x ) , x ) = 0 . Therefore, we expect to construct a localized sta-tionary state in a region of parameter space where H ( u , , , , 0) = H ( u , , , , . With a slight abuseof notation, we write this condition in terms of the vari-able r , as H ( r ) = H ( r ) , where H is given by H ( r ) = − ηr − Jr + 43 r . (21)In analogy with the literature mentioned above, we called Maxwell points the values on the ( η, J ) -plane where thecondition H ( r ) = H ( r ) is met. We display the Maxwellpoint for our standard parameter set in Figure 2b, andwe plot the locus of Maxwell points and the bistabilityregion in Figure 2c. III. NUMERICAL SCHEMES As anticipated in the previous sections, stationarystates beyond onset are computed numerically, hence wepresent in this section several numerical schemes used inthe upcoming computations. In preparation for present-ing the schemes, we rewrite the model as an ODE on afunction space. To simplify the notation we apply in thissection the scaling r (cid:55)→ r/π , J (cid:55)→ πJ to (3), and obtain ˙ r = ∆ + 2 rv, ˙ v = η + v − r + W r, (22)where W is the integral operator defined as ( W r )( x ) = J ( w ⊗ r )( x ) = J (cid:82) Ω w ( x, y ) r ( y ) d y . In the system abovewe assume r, v : R → L ρ (Ω , C ) (also denoted by L ρ (Ω) ),that is, at each time t , r ( t ) and v ( t ) belong to a weightedLebesgue space of complex-valued functions defined on Ω , with inner product (cid:104) f, g (cid:105) ρ = (cid:90) Ω f ( x ) g ∗ ( x ) ρ ( x ) d x, and norm (cid:107) f (cid:107) ρ = (cid:104) f, f (cid:105) / . Note that the subscript ρ will be omitted when ρ ( x ) ≡ . We assume that, oncecomplemented with initial conditions, system (22) definesa well-posed Cauchy problem on L ρ (Ω) × L ρ (Ω) .umps and Oscillons in Spiking Networks 5 FIG. 2. (a) Solution branches of uniform solutions around the bistable regime with saddle-node bifurcations (S). (b) Theassociated conserved quantity H , as defined in (21), showing a Maxwell point (M) at η M ≈ − . . (c) A plot of the bistableregion of homogeneous states (shaded) on the ( J, η ) -plane, delimited by saddle-node bifurcations (S), emanating from a cusp(C). We also show the locus of Maxwell points (M) on the plane. Parameters: J = 15 √ , ∆ = 2 . A. Galerkin Schemes Galerkin schemes are derived by introducing a com-plete orthogonal basis { ϕ i : i ∈ N } for the weightedspace L ρ (Ω) , and seeking an approximation in the n -dimensional subspace spanned by { ϕ i : i ∈ Λ n } , where Λ n is an index set with n elements, as follows r n ( x, t ) = (cid:88) i ∈ Λ n R i ( t ) ϕ i ( x ) , v n ( x, t ) = (cid:88) i ∈ Λ n V i ( t ) ϕ i ( x ) . A Galerkin scheme for (22) is then given by (cid:104) ϕ i , − ˙ r n + ∆ + 2 r n v n (cid:105) ρ = 0 , (cid:104) ϕ i , − ˙ v n + η + v n − r n + W r n (cid:105) ρ = 0 , i ∈ Λ n , that is, ˙ R i = α i ∆ + 2 (cid:88) j,k ∈ Λ n γ ijk R j V k , ˙ V i = α i η + (cid:88) j ∈ Λ n β ij R j + 2 (cid:88) j,k ∈ Λ n γ ijk ( V j V k − R j R k ) , for i ∈ Λ n , with coefficients given by α i = (cid:104) ϕ i , (cid:105) ρ , β ij = (cid:104) ϕ i , W ϕ j (cid:105) ρ , γ ijk = (cid:104) ϕ i , ϕ j ϕ k (cid:105) ρ . 1. Fourier-Galerkin Scheme When Ω = ( − L/ , L/ ∼ = S , the functions r ( t ) and v ( t ) are L -periodic. Therefore we choose the Fourier ba-sis ϕ j ( x ) = exp(i j πx/L ) , j ∈ Z , which is a completeorthogonal basis for L (Ω) . The index set for this case is Λ n = {− n/ , . . . , n/ − } with n even. Exploiting thetrigonometric properties of the Fourier basis, we obtain α i = (cid:40) L if i = 0 , otherwise, γ ijk = (cid:40) L if i + j + k = 0 , otherwise. In passing we note that β ij can also be expressed com-pactly, in terms of the Fourier coefficients of the kernel w ,if the operator W is convolutional. In addition, requiring r and v to be real-valued implies ( R i , V i ) = ( R ∗− i , V ∗− i ) .We call this method the Fourier-Galerkin scheme . 2. Hermite-Galerkin Scheme When Ω = R , a natural basis for the Galerkin schemeis given by the Hermite polynomials ϕ j ( x ) = H j ( x ) = ( − j exp( x ) d j d x j exp( − x ) , which are a complete orthogonal set for L ρ (Ω , R ) withweight ρ ( x ) = exp( − x ) . For this scheme Λ n = { , . . . , n − } . To avoid problems with the numericalevaluations of ϕ j for large | x | , we derive an alternativescheme, which uses inner products with weight ρ ( x ) ≡ ,as the Fourier Galerkin scheme. We seek a solution to(22) in the form r = R + ˜ r, v = V + ˜ v, with R , V constant in x , and ˜ r, ˜ v ∈ L (Ω , R ) . Thisleads to the system ˙ R = ∆ + 2 R V , ˙ V = η + V − R + JR , ˙˜ r = 2 R ˜ v + 2 V ˜ r + 2˜ r ˜ v, ˙˜ v = 2 V ˜ v − R ˜ r + ˜ v − ˜ r + W ˜ r, in which the homogeneous background dynamics for ( R , V ) is decoupled from (˜ r, ˜ v ) , and follows thespatially-clamped QIF mean field . Since the Hermitefunctions ϕ j ( x ) = exp( − x / H j − ( x ) , j ∈ N > , umps and Oscillons in Spiking Networks 6are an orthogonal set for L (Ω , R ) , an approximation to ˜ r, ˜ v is sought in the space spanned by ϕ j , with j ∈ Λ n = { , . . . , n } , giving the scheme ˙ R = ∆ + 2 R V , ˙ V = η + V − R + JR , ˙ R i = 2 (cid:88) j ∈ Λ n ( R V j + V R j ) + 2 (cid:88) j,k ∈ Λ n γ ijk R j V k , ˙ V i = (cid:88) j ∈ Λ n [2 V V j + ( β ij − R ) R j ]+ (cid:88) j,k ∈ Λ n γ ijk ( V j V k − R j R k ) , for i ∈ Λ n . We call this method the Hermite–Galerkinscheme . B. Fourier Collocation Scheme A Fourier collocation scheme can be derived in the case Ω = ( − L/ , L/ ∼ = S . This method, which has been usedin the past for Amari neural field models and the QIFneural field model , represents ( r n , v n ) by its values atthe gridpoints x j = − L + 2 Lj/n , j ∈ Λ n = { , . . . , n } , ˙ R i = ∆ + 2 R i V i , ˙ V i = η + V i + R i + ( W r n ) i , and evaluates ( W r n ) i either with a quadrature rule or,more efficiently, with a pseudospectral evaluation if W isconvolutional. C. Numerical considerations To the best of our knowledge, the methods presentedabove are novel, and we leave the analysis of the numeri-cal properties of these schemes to a separate publication.The calculations presented here have been tested againstevent-driven simulations of large network of spiking neu-rons. We employ our schemes as follows: the Fouriercollocation scheme with n = 5000 is generally used fortime simulations, to obtain accurate initial guesses for thecontinuation. However, we observed that time-periodicorbits are reproduced with a similar accuracy by theHermite–Galerkin scheme with just n = 50 modes, hencewe select this scheme to continue periodic orbits. Fi-nally, we use the Fourier–Galerkin scheme with n = 200 for bifurcation analysis of steady states on large domains,when solutions are non-localized.We compare the results of the QIF neural field modelwith the dynamics of the underlying network of spik-ing neurons. We integrate equation (1) using the Eulermethod, with time step ∆ t = 10 − . The domain is cho-sen to be Ω = ( − L/ , L/ with periodic boundary con-ditions. We choose L = 50 and × model neurons, which ensures a good correspondence to the neural fieldmodel ( δx = 10 − ); with L = 50 the relative error of thenormalization of the synaptic kernel, (cid:82) Ω w ( y ) d y = 1 , isless than − . We note here that the microscopic de-scription only matches the mean field description if δx ismuch smaller than the characteristic length scale of thesynaptic kernel w . The synaptic input to each neuronis computed with equation (2). We follow reference incomputing the synaptic integration across a time win-dow with a τ ( t ) = Θ( τ − t ) /τ , τ = 10 − , and in setting V θ = − V r = 100 with a refractory period of /V i once the i th neuron has exceeded V θ . The latter approximates thelimit V θ = − V r → ∞ . The refractory period is roundedto the nearest multiple of ∆ t , and after the refractoryperiod V i is set to − V i . Rasterplots are generated usinga subset of randomly chosen model neurons. IV. STATIONARY LOCALIZED SOLUTIONS We use the numerical schemes presented in the pre-vious section to study the bifurcation structure of sta-tionary localized solution to the QIF mean field model.We initially study the model with our default excitatory-inhibitory kernel (5), and then show that a snaking bifur-cation scenario is supported when the kernel is switchedto a homogeneous oscillatory kernel, or to a kernel withharmonic heterogeneities, similarly to what is found forAmari neural field models. A. Local excitation, lateral inhibition kernel We set w as in (5), generate a stationary localized so-lution by numerically integrating the model equations intime, and then implement the Fourier–Galerkin schemeto continue the localized solutions in η , using AUTO.In Figure 3 we show the bifurcation diagram of local-ized solutions. Across a range of parameters, these occuras a pair of one wide, stable solution and one narrow,unstable solution. The solution branch connects to thebranch of uniform solutions at points where Turing bifur-cations occur, which also give rise to a branch of periodicsolutions. Using the Fourier basis, it can be shown thatthe stable solution branch approaches the Maxwell pointasymptotically, and solutions grow wider, which resembletwo (stationary) interacting wave fronts. Because of theperiodic boundary conditions, the solution branch growslarger again and forms another stable/unstable solutionpair of locally low activity (not shown). The latter couldbe regarded as stationary versions of traveling anti-pulses reported in refs. .Because stable solutions are of particular interest, wepresent a two-parameter bifurcation diagram (Figure 4a)of the saddle-node bifurcations that delimit the branchof stable solutions. As expected, the locus of saddle-node bifurcations of localized states enclose the Maxwellpoint. In addition, the loci of saddle-node bifurcations ofumps and Oscillons in Spiking Networks 7 FIG. 3. (a) Bifurcation diagram in η of localized solu-tions (blue), periodic solutions (green), and uniform solutions(black). (b) Exemplary profiles of unstable narrow (1), stablenarrow (2), unstable wide (3), and stable periodic (4) local-ized solutions, close to the Maxwell point (vertical line in (a)).Parameters: ∆ = 2 , J = 15 √ ∆ . localized and uniform steady states meet at two separatecusps, as shown in Figure 4b.The bifurcation behavior of localized solutions de-scribed above is robust to changes in coupling param-eters but, as we shall see below, it is strongly affected bychanges in the kernel. B. Snaking with homogeneous kernel Homoclinic snaking is a phenomenon that describesthe formation of multiple, coexisting localized solutionsin spatially-extended models. Steady states are ar-ranged in branches of intertwined snaking bifurcation di-agrams, connected via ladders . Adopting thespatial-dynamics approach outlined above, localized so-lutions are interpreted as homoclinic orbits to a fixedpoint. Snaking solution branches correspond to sym-metries of the problems, which are broken along theladder branches . This scenario is not limited toPDEs, but have also been studied in the non-local Swift-Hohenberg equation , as well as in neural field mod-els .In the simplest setting, localized snaking solutionsare found in regions of parameter space where there isbistability between a stationary homogeneous state anda periodic state. In nonlocal neural fields, homoclinicsnaking has been observed with the following homoge-neous damped-oscillatory kernel w ( x ) = 1 + b b e − b | x | ( b sin | x | + cos x ) , (23)which we now adopt also for the QIF neural field model.This kernel leads to a sub-critical Turing bifurcationof the lower stable branch of uniform solutions, fromwhich an unstable branch of spatially-periodic solutionsemerges. This branch undergoes a saddle-node bifur-cation, where spatially-periodic solutions become sta-ble. Eventually, the branch connects to the upper stablebranch of uniform solutions, see Figure 5. FIG. 4. (a) The parameter space in which stable localizedsolutions exist is delimited by loci of saddle-node bifurcations.They can be approximated by the saddle-node bifurcations ofspatially uniform solutions, and the Maxwell point. (b) Insetof (a), showing additionally the loci of Turing bifurcations andsaddle-node bifurcations of periodic solutions. The saddle-node bifurcations of the bump solutions form a cusp where theTuring bifurcation changes from supercritical to subcritical.Parameters: ∆ = 2 .FIG. 5. (a) The damped-oscillatory kernel (23). (b) Fouriertransform of the kernel, which shows a maximum at a non-zero wave-number k c . (c) As η is varied, the branch of homo-geneous steady states (black) undergoes a Turing bifurcation,from which a branch of periodic solutions with wavenumber k c emerges (green). Parameters: J = 39 , ∆ = 1 , b = 0 . . umps and Oscillons in Spiking Networks 8 FIG. 6. (a) Snakes-and-ladders bifurcation scenario. (b)Representative solutions for branches of symmetric (1,3) andasymmetric solutions (2). Parameters: J = 39 , ∆ = 1 , b =0 . . As anticipated, spatially localized snaking solutions arefound in this region of parameter space, and they are ar-ranged in a typical snakes-and-ladders bifurcation struc-ture, which is displayed in Figure 6. C. Snaking with heterogeneous kernel It is known that snaking bifurcation scenarios canbe triggered by heterogeneities in the underlying evo-lution equations. Examples discussed in the literatureinclude the Swift–Hohenberg , Amari , and Ginzburg-Landau equations. In neural field models, hetero-geneities are naturally introduced via harmonic pertur-bations of a homogeneous (distance-dependent) kernel,which break the translational invariance of the prob-lem . In Reference we have shown that the fol-lowing kernel leads to snaking in the Amari model w ( x, y ) = 12 e −| x − y | (1 + a cos( ky )) , (24)and we therefore investigate the effect of this kernel onthe QIF neural field model.In the absence of spatial forcing ( a = 0 ), a system withexponential connectivity does not yield stable localizedsolutions (see Figure 7). In the presence of modulation,we find snaking branches that oscillate around the branchobtained for a = 0 (see Figure 8). Furthermore, for smallvalues of a , the snaking width increases proportionally tothe value of a (not shown). These findings indicate thatthe snaking phenomenon in the QIF neural field model isentirely determined by the kernel choice, as in the Amaricase. V. OSCILLONS Various nonlinear models including chemical, fluid-dynamical, and particle systems, support time-periodic, FIG. 7. (a) Bifurcation diagram of spatially uniform solu-tions, and of localized solutions generated with the exponen-tial kernel (24) with a = 0 . The lack of lateral inhibitionresults in the entire branch of solutions being unstable. (b)Representative solutions. Parameters: J = 15 √ ∆ , ∆ = 2 .FIG. 8. Snaking induced by spatial periodic modulationof the exponential connectivity kernel (24). (a) Bifurcationdiagram obtained for a = 0 . (orange and purple branches,ladders not shown) compared with the one obtained for a = 0 (blue). (b) Representative solutions. Parameters: J = 15 √ ∆ , ∆ = 2 . spatially-localized states termed oscillons (see Refer-ence and references therein). A comprehensive the-ory for the existence and bifurcation structure of suchsolutions is the subject of experimental, numerical, andanalytical investigations. We study oscillons in the QIFneural field model in the two main settings where theyare observed in other media: (i) a non-autonomous set-ting, whereby oscillons emerge as the medium is subjectto a homogeneous, exogenous, time-periodic forcing; (ii)an autonomous setting, whereby oscillons emerge spon-taneously as one of the model parameters is varied. A. Oscillons induced by harmonic forcing We setup the QIF neural field model subject to atime-dependent, homogeneous, sinusoidal forcing withfrequency ω , ∂ t r = ∆ π + 2 rv,∂ t v = v + Jw ⊗ r − π r + η + A sin( ωt ) , umps and Oscillons in Spiking Networks 9and cast it in the following, equivalent autonomous modelformulation to perform numerical bifurcation analysis ∂ t r = ∆ π + 2 rv,∂ t v = v + Jw ⊗ r − π r + η + Aξ, ˙ ξ = ξ + ωζ − ( ξ + ζ ) ξ, ˙ ζ = ζ − ωξ − ( ξ + ζ ) ζ. (25)Note that the numerical framework proposed here is ap-plicable also if the forcing is heterogeneous.In this setting we expect oscillons to emerge without bi-furcation from a localized steady state of the QIF neuralfield model with A = 0 , upon imposing a small-amplitudeforcing, A (cid:28) . We therefore select the default kernel (5),set η = − , for which the model with A = 0 supportsone stable (wide) and one unstable (narrow) bump (seeFigure 3), and continue time-periodic solutions to (25)in A > for ω = 4 , close to the network’s resonant fre-quency .One stable and one unstable branch of oscillons emergefrom A = 0 , as shown in Figure 9, and connect at asaddle-node bifurcation. The stable branch undergoesa sequence of period-doubling bifurcations leading tochaos, and examples of a period-doubled solution anda chaotic solution are shown in Figure 9, demonstratingthe correspondence between the QIF neural field modeland the spiking network model.In a recent study we have investigated the effect of pe-riodic forcing on a population of excitatory spiking neu-rons , whose solutions correspond to the spatially uni-form states of the present model. In that context it wasshown that a sufficiently large forcing amplitude is ableto suppress homogeneous oscillations. Here we reportthat the same statement holds true for forced oscillons:no localized time-periodic solution is found to the rightof the saddle-node bifurcation in Figure 9, where the at-tractor is a spatially-homogeneous, time-periodic state,which can be found by continuing in A the low-activityuniform steady state (not shown). B. Spontaneous oscillons in coupled networks of excitatoryand inhibitory neurons In the second scenario, oscillons occur in autonomoussystems. Direct numerical simulations of reaction dif-fusion systems display oscillons in the proximity ofcodimension-two Turing–Hopf bifurcations of the homo-geneous steady state . Oscillons in these systems havetypically been observed as large-amplitude structures,hence they are conjectured to form via a subcritical Hopfbifurcation of a heterogeneous, spatially-localized steadystate. This conjecture, however, has not yet been con-firmed by numerical bifurcation analysis which, in con-trast to direct numerical simulations, allows to track bothstable and unstable states. Here we employ the Hermite–Galerkin scheme to studythe formation of oscillons in the QIF neural field model.As mentioned above, a necessary ingredient for oscil-lons is the presence of oscillatory bifurcations. Thesebifurcations are precluded in one-populations networksof QIF neurons but, as we shall see, are possible in two-population models, therefore we turn our attention to thefollowing network of coupled excitatory and inhibitorypopulations ˙ r e = ∆ π + 2 r e v e , ˙ v e = v e + η e + J e w e ⊗ r e − J i τ i w i ⊗ r i − π r e ,τ i ˙ r i = ∆ π + 2 τ i r i v i ,τ i ˙ v i = v i + η i + J e w e ⊗ r e − J i τ i w i ⊗ r i − π τ i r i . The subscripts e , i indicate whether a variable or param-eter refers to the excitatory or inhibitory population, re-spectively: the two populations have, for simplicity, thesame heterogeneity parameter ∆ , but they have possi-bly different membrane time constants and average back-ground currents. In single-population mean fields, exci-tation and inhibition are artificially lumped into a singleexcitatory-inhibitory kernel (see for instance (5), (23),and (24)), whereas in the new, more realistic model thekernels are separate w e ( x ) = e −| x | , w i ( x ) = 14 e −| x | / . (26)The connectivity parameters are chosen to be J e = J i = J to recover a similar setting used in the lumped model.In Figure 10a we show the bifurcation diagram of local-ized solutions using η e as bifurcation parameter. The bi-furcation structure is similar to the lumped model, withthe exception that the range of parameters for which sta-ble solutions exist is narrower. This computation con-firms that stationary bumps are supported by the two-population network. In order to hunt for oscillons, wecontinue the solution for η e = − in the parameter τ i :the bump becomes unstable at a subcritical Hopf bifur-cation at τ i ≈ . , restabilizes at a saddle-node bifur-cation, and undergoes a sequence of saddle-node bifur-cations leading to a torus bifurcation (i.e. generalizedHopf bifurcation). The branch eventually restabilizes ata further saddle-node, leading to a period-doubling cas-cade which initiates around τ i ≈ . , and to chaos at τ i > . (Figure 10b).In Figure 10 we also show numerical examples of astable period-doubled solution at τ i = 1 . and a chaoticsolution at τ i = 1 . . We do not observe oscillons be-yond τ i = 1 . . Chaotic solutions can also be reproducedin the spiking network model, see Figure 11. VI. DISCUSSION We introduced a framework to study localized solu-tions in a neural field model that was recently derivedumps and Oscillons in Spiking Networks 10 FIG. 9. Maximum values of R ( t ) plotted against the amplitude of sinusoidal forcing. At A ≈ . a period-doubling bifurcationoccurs, which is the starting point of a period-doubling cascade leading to chaos at A (cid:46) . . A period-doubled solution ( A = 4 )and a chaotic solution ( A = 4 . ) is shown. Parameters: ∆ = 2 , J = 15 √ ∆ , η = − , ω = 4 .FIG. 10. (a) Bifurcation diagram of bump solutions in the E-I network (blue), and spatially uniform solutions (black) for τ i = 1 . (b) Bifurcation diagram of emerging limit cycles (showing maxima of R ( t ) ) using τ i as bifurcation parameter. H: Hopfbifurcation, S: saddle node bifurcation, T: torus bifurcation, P: period doubling bifurcation. Parameters: η e = η i = − , J e = J i = 15 √ ∆ , ∆ = 2 . as an exact representation of the mean field dynamics ofnetworks of spiking neurons. Although this model doesnot permit closed-form solutions such as the Amari modelwith Heaviside firing rates, we show that it is possible togive an analytical estimate for the range of model param-eters for which stable localized solutions exist. The struc-ture of the QIF neural field model permits the straight- forward use of Galerkin methods, which unlike the Amarimodel has a linear nonlocal term.We have demonstrated that stationary equations canbe transformed into a formulation that is equivalent tothe stationary Amari model, provided an effective firing-rate function is defined. The significance of such a firingrate is chiefly mathematical: the neural field possesses aumps and Oscillons in Spiking Networks 11 t x r e r i tt xr r i r e (a) (b)(c) (d) FIG. 11. Chaotic solution in spiking network. (a) Raster-gram of excitatory population. (b) Rastergram of inhibitorypopulation. (c) Excitatory and inhibitory spike rates aver-aged on interval − < x < and sliding window in t (width − ). (d) Phase portrait of spike rates in (c). Parameters: η e = η i = − , J e = J i = 15 √ ∆ , ∆ = 2 , τ i = 1 . . Theseparameters correspond to point (4) in Figure 10. rate variable, which is combined with the voltage vari-able in the effective firing rate; however, this transfor-mation allows to map out patterned steady states of theQIF neural field model using the same toolkit availablefor the Amari formulation. In both models localized so-lutions emerge subcritically from a branch of homoge-neous steady states, which then restabilize at a saddlenode bifurcation. In the Amari model, this behavioris parametrized by a firing threshold, whereas here weuse the average excitability of the network to map outsolutions. However, there is a correspondence betweenthe excitability of the model used here and the firingthreshold in the Amari model, in the sense that an in-crease in the firing threshold in the latter corresponds toa decrease in the excitability in the former. In addition,techniques developed for piecewise-linear firing rate func-tions in the Amari model , could be adapted to work forsteady states in the QIF neural field model, using the cor-respondence described above. Furthermore, all branchesof stationary solutions computed in this paper, includingthe snaking branches, also occur in standard rate-basedmodels. The crucial difference lies in the transient dy-namics of the two models, which makes the model con-sidered here dynamically richer and more realistic.The development of a Galerkin method opened up thepossibility to study oscillons using numerical bifurcationanalysis. We focused here on sinusoidal forcing of bumpsolutions, which is a proxy of oscillations ubiquitous inneuronal systems. In previous work , the neural massversion of this model was studied in terms of its response to oscillatory forcing in various frequency bands, and thepresent paper makes this exploration feasible also in thespatially-extended model. We leave this exploration to afuture publication.In coupled networks of excitatory and inhibitory pop-ulations, a small change in the inhibitory membrane timescale can have a significant effect on the existence and dy-namics of bump solutions, and can elicit oscillons. Thiswas demonstrated for instantaneous synapses, and it re-mains to be seen how the dynamics changes when synap-tic delays are introduced to the model. Interestingly, os-cillatory solutions which undergo torus bifurcations havebeen observed in spatially extended networks of excita-tory and inhibitory neurons with conductance-based dy-namics . Another natural extension would be to exam-ine coupled multi-layer neural field models , which areknown to give rise to localized bump solutions when nei-ther layer does in isolation .The Galerkin numerical methods derived in this pa-per can be applied directly to more general spatially-extended models of QIF networks, such as the ones men-tioned above. For instance, adding a synaptic variablecan be accounted for with an additional Galerkin expan-sion, and n scalar variables per additional evolution equa-tion.Single population, QIF neural mass models with chem-ical as well as electrical synapses have recently been de-veloped , and it was found that oscillations originateat Hopf bifurcations. Spatially-extended versions of thismodel would then have the possibility of forming oscil-lons with a single population, although it is not clearwhether Hopf bifurcations of bumps will occur near Hopfbifurcations of homogeneous states, which are the onesmapped in Reference .Understanding how slow-fast temporal scales are gen-erated by the discrete network is an open question, whichhas recently been addressed in networks of sparsely-coupled networks of QIF neurons . Employing our nu-merical methodology to these macroscopic mean fields isalso possible, and one could study how such slow-fastphenomena occur in more realistic, spatially-extendednetworks. Acknowledgments HS acknowledges financial support from the SpanishMinistry of Economics and Competitiveness through theMaría de Maeztu Programme for Units of Excellencein R&D (MDM-2014-0445) and grant MTM2015-71509-C2-1-R, and from the German Research Council (DFG(KN 588/7-1) within priority programme ‘ComputationalConnectomics’ (SPP 2041) ).umps and Oscillons in Spiking Networks 12 References A. Compte, N. Brunel, P. S. Goldman-Rakic, and X.-J. Wang,“Synaptic mechanisms and network dynamics underlying spatialworking memory in a cortical network model.” Cereb. Cortex ,910 – 923 (2000). K. Wimmer, D. Q. Nykamp, C. Constantinidis, and A. Compte,“Bump attractor dynamics in prefrontal cortex explains behav-ioral precision in spatial working memory.” Nat. Neurosc. , 431– 439 (2014). S. Kim, H. Rouault, S. Druckmann, and V. Jayaraman, “Ringattractor dynamics in the drosophila brain.” Science , 849 –853 (2017). H. R. Wilson and J. D. Cowan, “A mathematical theory of thefunctional dynamics of cortical and thalamic nervous tissue.” Ky-bernetik , 55 – 80 (1973). P. L. Nunez, “The brain wave equation: A model for the EEG,”Mathematical Biosciences , 279 – 297 (1974). S. Amari, “Dynamics of pattern formation in lateral-inhibitiontype neural fields,” Biol. Cybern. , 77–87 (1977). S. Coombes, “Waves, bumps, and patterns in neural field theo-ries.” Biol. Cybern. , 91 – 108 (2005). P. C. Bressloff, “Spatiotemporal dynamics of continuum neuralfields.” J. Phys. A , 033001 (2012). G. B. Ermentrout, “Reduction of conductance-based models withslow synapses to neural nets.” Neural Computation , 679 – 695(1994). S. Ostojic and N. Brunel, “From spiking neuron models to linear-nonlinear models.” PLOS Comp. Biol. , e1001056 (2011). E. S. Schaffer, S. Ostojic, and L. F. Abbott, “A complex-valuedfiring-rate model that approximates the dynamics of spiking net-works,” PLOS Comp. Biol. , e1003301 (2013). M. A. Buice and C. C. Chow, “Dynamic finite size effects inspiking neural networks.” PLOS Comp. Biol. , e1002872 (2013). S. Visser and S. A. V. Gils, “Lumping Izhikevich neurons,” EPJNonlinear Biomed , 6 (2014). M. Mattia, “Low-dimensional firing rate dynamics of spiking neu-ron networks.” arXiv , 160908855 (2016). T. Schwalger, M. Deger, and W. Gerstner, “Towards a theory ofcortical columns: from spiking neurons to interacting neural pop-ulations of finite size.” PLOS Comp. Biol. , e1005507 (2017). M. Augustin, J. Ladenbauer, F. Baumann, and K. Obermayer,“Low-dimensional spike rate models derived from networks ofadaptive integrate-and-fire neurons: comparison and implemen-tation.” PLOS Comp. Biol. , e1005545 (2017). Y. Park and G. B. Ermentrout, “A multiple timescales approachto bridging spiking- and population-level dynamics,” Chaos ,083123 (2018). S.-W. Qiu and C. C. Chow, “Finite-size effects for spiking neuralnetworks with spatially dependent coupling,” Phys. Rev. E ,062414 (2018). T. B. Luke, E. Barreto, and P. So, “Complete classification ofthe macroscopic behavior of a heterogeneous network of thetaneurons.” Neural Comput. , 3207 – 3234 (2013). E. Montbrió, D. Pazó, and A. Roxin, “Macroscopic descriptionfor networks of spiking neurons,” Phys. Rev. X (2015). C. R. Laing, “Derivation of a neural field model from a networkof theta neurons,” Phys. Rev. E , 010901 (2014). C. R. Laing, “Exact neural fields incorporating gap junctions,”SIAM Journal on Applied Dynamical Systems , 1899–1929(2015). J. M. Esnaola-Acebes, A. Roxin, D. Avitabile, and E. Montbrió,“Synchrony-induced modes of oscillation of a neural field model,”Phys. Rev. E (2017). A. Byrne, D. Avitabile, and S. Coombes, “Next-generation neu-ral field model: The evolution of synchrony within patterns andwaves,” Phys. Rev. E , 012313 (2019). S. Coombes and À. Byrne, “Nonlinear dynamics in computationalneuroscience,” (Springer, 2019) Chap. Next generation neural mass models. F. Devalle, A. Roxin, and E. Montbrió, “Firing rate equa-tions require a spike synchrony mechanism to correctly describefast oscillations in inhibitory networks,” PLOS Comp. Biol. ,e1005881 (2017). I. Ratas and K. Pyragas, “Macroscopic oscillations of a quadraticintegrate-and-fire neuron network with global distributed-delaycoupling,” Phys. Rev. E , 052224 (2018). R. Zillmer, R. Livi, A. Politi, and A. Torcini, “Stability of thesplay state in pulse-coupled networks,” Phys. Rev. E , 046102(2007). E. Ott and T. M. Antonsen, “Low dimensional behavior of largesystems of globally coupled oscillators,” Chaos , 037113 (2008). H. Schmidt, D. Avitabile, E. Montbrió, and A. Roxin, “Networkmechanisms underlying the role of oscillations in cognitive tasks,”PLOS Comp. Biol. , e1006430 (2018). Y. Kawamura, “From the Kuramoto-Sakaguchi model to theKuramoto-Sivashinsky equation,” Phys. Rev. E , 010901(2014). A. R. Champneys, “Homoclinic orbits in reversible systems andtheir applications in mechanics, fluids and optics,” Physica D:Nonlinear Phenomena , 158 – 186 (1998), proceedings of theWorkshop on Time-Reversal Symmetry in Dynamical Systems. J. Burke and E. Knobloch, “Homoclinic snaking: Structure andstability,” Chaos , 037102 (2007). J. Burke and E. Knobloch, “Snakes and ladders: localized statesin the Swift–Hohenberg equation,” Physics Letters A , 681 –688 (2007). J. Knobloch and T. Wagenknecht, “Homoclinic snaking near aheteroclinic cycle in reversible systems,” Physica D , 82 – 93(2005). M. Beck, J. Knobloch, D. J. B. Lloyd, B. Sandstede, and T. Wa-genknecht, “Snakes, ladders, and isolas of localised patterns.”SIAM J. Math. Anal. , 936 – 972 (2009). C. R. Laing, W. C. Troy, B. Gutkin, and G. B. Ermentrout,“Multiple bumps in a neuronal model of working memory,” SIAMJ. Appl. Math. , 62 – 97 (2002). S. Coombes, G. J. Lord, and M. R. Owen, “Waves and bumpsin neuronal networks with axo-dendritic synaptic interactions,”Physica D , 219 – 241 (2003). A. J. Elvin, C. R. Laing, R. I. McLachlan, and M. G. Roberts,“Exploiting the Hamiltonian structure of a neural field model,”Physica D , 537 – 546 (2010). D. Avitabile and H. Schmidt, “Snakes and ladders in an inhomo-geneous neural field model,” Physica D , 24 – 36 (2015). J. Rankin, D. Avitabile, J. Baladron, G. Faye, and D. J. B.Lloyd, “Continuation of localised coherent structures in nonlocalneural field equations,” SIAM J. Sci. Comput. , B70 – B93(2013). C. R. Laing and S. Coombes, “The importance of different tim-ings of excitatory and inhibitory pathways in neural field models,”Network , 151 – 172 (2006). H. G. Meijer and S. Coombes, “Travelling waves in models ofneural tissue: from localised structures to periodic waves,” EPJNonlinear Biomed Phys , 3 (2014). A. Champneys, “Homoclinic orbits in reversible systems and theirapplications in mechanics, fluids and optics,” Physica D: Nonlin-ear Phenomena , 158 – 186 (1998), Proceedings of the Work-shop on Time-Reversal Symmetry in Dynamical Systems. D. J. B. Lloyd, B. Sandstede, D. Avitabile, and A. R. Champ-neys, “Localized hexagon patterns of the planar Swift–Hohenbergequation,” SIAM J. Appl. Dyn. Syst. , 1049 – 1100 (2008). D. Avitabile, D. J. B. Lloyd, J. Burke, E. Knobloch, and B. Sand-stede, “To snake or not to snake in the planar Swift-Hohenbergequation,” SIAM J. Appl. Dyn. Syst. , 704 – 733 (2010). D. Morgan and J. H. P. Dawes, “The Swift–Hohenberg equationwith a nonlocal nonlinearity,” Physica D , 60 – 80 (2014). C. R. Laing and W. C. Troy, “PDE methods for nonlocal models,”SIAM J. Appl. Dyn. Syst. , 487 – 516 (2003). umps and Oscillons in Spiking Networks 13 G. Faye, J. Rankin, and P. Chossat, “Localized states in an un-bounded neural field equation with smooth firing rate function: amulti-parameter analysis,” J. Math. Biol. , 1303 – 1338 (2013). G. Faye, J. Rankin, and D. J. B. Lloyd, “Localized radial bumpsof a neural field equation on the euclidean plane and the poincarédisk,” Nonlinearity , 437 (2013). H.-C. Kao, C. Beaume, and E. Knobloch, “Spatial localizationin heterogeneous systems,” Physical Review E , 012903 (2014). B. C. Ponedel and E. Knobloch, “Forced snaking: Localizedstructures in the real Ginzburg-Landau equation with spatiallyperiodic parametric forcing,” Eur. Phys. J. Spec. Top. , 2549– 2561 (2016). P. C. Bressloff, “Traveling fronts and wave propagation failurein an inhomogeneous neural network,” Physica D , 83 – 100(2001). Z. P. Kilpatrick, S. E. Folias, and P. C. Bressloff, “Travelingpulses and wave propagation failure in inhomogeneous neural me-dia,” SIAM J. Appl. Dyn. Syst. , 161 – 185 (2008). H. Schmidt, A. Hutt, and L. Schimansky-Geier, “Wave fronts ininhomogeneous neural field models,” Physica D , 1101 – 1112(2009). C. Coombes and C. R. Laing, “Pulsating fronts in periodicallymodulated neural field models,” Phys. Rev. E , 011912 (2011). E. Knobloch, “Spatially localized structures in dissipative sys-tems: open problems,” Nonlinearity , T45 – T60 (2008). V. Vanag and I. Epstein, “Stationary and Oscillatory LocalizedPatterns, and Subcritical Bifurcations,” Physical Review Letters , 128301 (2004). V. K. Vanag and I. R. Epstein, “Localized patterns in reaction-diffusion systems,” Chaos: An Interdisciplinary Journal of Non-linear Science , 037110 (2007). S. Coombes and H. Schmidt, “Neural fields with sigmoidal firingrates: Approximate solutions.” Discrete and Continuous Dynam-ical Systems - A , 1369 – 1379 (2010). S. E. Folias and G. B. Ermentrout, “Spatially localized syn-chronous oscillations in synaptically coupled neuronal networks:conductance-based models and discrete maps.” SIAM J. Appl.Dyn. Syst. , 1019–1060 (2010). P. C. Bressloff and S. R. Carroll, “Laminar neural field modelof laterally propagating waves of orientation selectivity,” PLOSComp. Biol. , e1004545 (2015). S. E. Folias and G. B. Ermentrout, “New patterns of activity ina pair of interacting excitatory-inhibitory neural fields,” Phys.Rev. Lett. , 228103 (2011). B. Pietras, F. Devalle, A. Roxin, A. Daffertshofer, and E. Mont-brió, “Exact firing rate model reveals the differential effects ofchemical versus electrical synapses in spiking networks,” Physi-cal Review E , 042412 (2019).65