Macroscopic description for networks of spiking neurons

A major goal of neuroscience, statistical physics and nonlinear dynamics is to understand how brain function arises from the collective dynamics of networks of spiking neurons. This challenge has been chiefly addressed through large-scale numerical simulations. Alternatively, researchers have formulated mean-field theories to gain insight into macroscopic states of large neuronal networks in terms of the collective firing activity of the neurons, or the firing rate. However, these theories have not succeeded in establishing an exact correspondence between the firing rate of the network and the underlying microscopic state of the spiking neurons. This has largely constrained the range of applicability of such macroscopic descriptions, particularly when trying to describe neuronal synchronization. Here we provide the derivation of a set of exact macroscopic equations for a network of spiking neurons. Our results reveal that the spike generation mechanism of individual neurons introduces an effective coupling between two biophysically relevant macroscopic quantities, the firing rate and the mean membrane potential, which together govern the evolution of the neuronal network. The resulting equations exactly describe all possible macroscopic dynamical states of the network, including states of synchronous spiking activity. Finally we show that the firing rate description is related, via a conformal map, with a low-dimensional description in terms of the Kuramoto order parameter, called Ott-Antonsen theory. We anticipate our results will be an important tool in investigating how large networks of spiking neurons self-organize in time to process and encode information in the brain.

Processing and coding of information in the brain necessarily imply the coordinated activity of large ensembles of neurons.Within sensory regions of the cortex, many cells show similar responses to a given stimulus, indicating a high degree of neuronal redundancy at the local level.This suggests that information is encoded in the population response and hence can be captured via macroscopic measures of the network activity [1].Moreover, the collective behavior of large neuronal networks is particularly relevant given that current brain measurement techniques, such as electroencephalography (EEG) or functional magnetic resonance imaging (fMRI), provide data which is necessarily averaged over the activity of a large number of neurons.
Despite these efforts, to date there is no exact theory linking the dynamics of a large network of spiking neurons with that of the firing rate.Specifically, current macroscopic descriptions do not offer a precise correspondence between the microscopic dynamics of the individual neurons, e.g.individual membrane potentials, and the firing rate dynamics of the neuronal network.
Indeed, FREs are generally derived through heuristic arguments which rely on the underlying spiking activity of the neurons being asynchronous and hence uncorrelated.As such, firing rate descriptions are not sufficient to describe network states involving some degree of spike synchronization.Synchronization is, however, an ubiquitous phenomenon in the brain, and its potential role in neuronal computation is the subject of intense research [14,[30][31][32][33][34][35][36].Hence, the lack of firing rate descriptions for synchronous states limits the range of applicability of mean-field theories to investigate neuronal dynamics.
Here we propose a method to derive the FREs for networks of heterogeneous, all-to-all coupled quadratic integrate-andfire (QIF) neurons, which is exact in the thermodynamic limit, i.e. for large numbers of neurons.We consider an ansatz for the distribution of the neurons' membrane potentials that we denominate the Lorentzian ansatz (LA).The LA solves the corresponding continuity equation exactly, making the system amenable to theoretical analysis.Specifically, for particular distributions of the heterogeneity, the LA yields a system of two ordinary differential equations for the firing rate and mean membrane potential of the neuronal population.These equations fully describe the macroscopic states of the networkincluding synchronized states-, and represent the first example of an exact firing-rate description for a network of recurrently connected spiking neurons.We finally show how the LA transforms, via a conformal mapping, into the so-called Ott-Antonsen ansatz (OA) that is extensively used to investigate the low-dimensional dynamics of large populations of phase oscillators in terms of the Kuramoto order parameter [37].

I. MODEL DESCRIPTION
Hodgkin-Huxley-type neuronal models can be broadly classified into two classes, according to the nature of their transition to spiking in response to an injected current [38,39].Neuronal models with so-called Class I excitability, generate action potentials with arbitrarily low frequency, depending on the strength of the applied current.This occurs when a resting state disappears through a saddle-node bifurcation.In contrast, in neurons with Class II excitability the action potentials are generated with a finite frequency.This occurs when the resting state loses stability via a Hopf bifurcation.The QIF neuron is the canonical model for Class I neurons, and thus generically describes their dynamics near the spiking threshold [5,40,41].Our aim here is to derive the FREs corresponding to a heterogeneous all-to-all coupled population of N QIF neurons.The correspondence is exact in the thermodynamic limit, i.e. when N → ∞ (this convergence has been recently studied in [42]).
The microscopic state of the population of QIF neurons is given by the membrane potentials {V j } j=1,...,N , which obey the following ordinary differential equations [5]: Here the overdot denotes the time derivative, and I j represents an input current.Each time a neuron's membrane potential V j reaches the peak value V p , the neuron emits a spike and its voltage is reset to the value V r .In our analysis we consider the limit V p = −V r → ∞.This resetting rule captures the spike reset as well as the refractoriness of the neuron.Without loss of generality we have rescaled the time and the voltage in (1) to absorb any coefficients which would have appeared in the first two terms.The form for the input currents is: where the external input has a heterogeneous, quenched component η j as well as a common time-varying component I(t), and the recurrent input is the synaptic weight J times the mean synaptic activation s(t), which is written: Here, t k j is the time of the kth spike of jth neuron, δ(t) is the Dirac delta function, and a τ (t) is the normalized synaptic activation caused by a single pre-synaptic spike with time scale τ , e.g. a τ (t) = e −t/τ /τ .

A. Continuous formulation
In the thermodynamic limit N → ∞, we drop the indices in Eqs. ( 1), (2) and denote ρ(V |η, t)dV as the fraction of neurons with membrane potentials between V and V + dV , and parameter η at time t.Accordingly, parameter η becomes now a continuous random variable distributed according to a probability distribution function g(η).The total voltage density at time t is then given by ∞ −∞ ρ(V |η, t) g(η) dη.The conservation of the number of neurons leads to the following continuity equation: where we have explicitly included the velocity given by ( 1) and (2).

II. RESULTS
The continuity equation ( 4) without temporal forcing I(t) = 0 has a trivial stationary solution.For each value of η, this solution has the form of a Lorentzian function: Physically, the Lorentzian density means that firing neurons with the same η value will be scattered with a density inversely proportional to their speed (1), i.e. they will accumulate at slow places and thin out at fast places on the V axis.In addition, for those η values corresponding to quiescent neurons, the density ρ 0 collapses at the rest state in the form of a Dirac delta function.
Next we assume that, independently of the initial condition, solutions of (4) generically converge to a Lorentzian-shaped function, so that all relevant dynamics occur inside that lower dimensional space.This fact is mathematically expressed by the following Lorentzian Ansatz (LA) for the conditional density functions: which is a Lorentzian function with time-dependent halfwidth x(η, t) and center at y(η, t).In the following we assume that the LA (5) completely describes the macroscopic dynamics of the network of QIF neurons and postpone the mathematical justification of its validity to Sec.II E.

A. Macroscopic observables: Firing rate and mean membrane potential
The half-width x(η, t) of the LA has a particularly simple relation with the firing rate of the neuronal population (i.e. the number of spikes per unit time).Indeed, the firing rate for each η value at time t, r(η, t), can be computed by noting that neurons fire at a rate given by the probability flux at infinity: r(η, t) = ρ(V → ∞|η, t) V (V → ∞|η, t).The limit V → ∞ on the right hand side of this equation can be evaluated within the LA, and gives the simple identity x(η, t) = πr(η, t).
The (total) firing rate r(t) is then r/∆ FIG. 1. (color online).Analysis of the steady states of FREs (12).(a) Phase diagram: In the wedge-shaped cyan-shaded region there is bistability between a high and a low activity state.The boundary of the bistability region is the locus of a saddle-node bifurcation which is exactly obtained in parametric form: To the right of the dashed line, defined by ηf = −[J/(2π)] 2 − (π∆/J) 2 , there is a stable focus (shaded regions).(b) r-η and (c) v-η bifurcation diagrams for J/∆ 1/2 = 15.Square symbols: Results obtained from numerical simulations of QIF neurons (see Appendix A for details).(d) Phase portrait of the system in the bistable region (η/∆ = −5, J/∆ 1/2 = 15, triangle in panel (a)) with three fixed points: a stable focus (with its basin of attraction shaded), a stable node, and a saddle point.
Additionally, the quantity y(η, t) is, for each value of η, the mean of the membrane potential: Here we take the Cauchy principal value, defined as P.V.
, to avoid the otherwise ill-defined integral.The mean membrane potential is then

B. Firing-rate equations
Substituting the LA (5) into the continuity equation ( 4), we find that, for each value of η, variables x and y must obey two coupled equations which can be written in complex form as where w(η, t) ≡ x(η, t) + iy(η, t).Closing this equation requires expressing the mean synaptic activation s(t) as a function of w(η, t).The simplest choice is to take the limit of infinitely fast synapses, τ → 0 in (3), so that we get an equality with the firing rate: s(t) = r(t).This allows for the system of QIF neurons (1)-(3) to be exactly described by Eqs. ( 10) and (7); an infinite set of integro-differential equations if g(η) is a continuous distribution.Equation ( 10) is useful for general distributions g(η) (see Appendix B), but a particularly sharp reduction in dimensionality is achieved if η is distributed according to a Lorentzian distribution of half-width ∆ centered at η: Note that this distribution accounts for the quenched variability in the external inputs.The fact that it is Lorentzian is unrelated to the LA for the density of membrane potentials.Assuming (11) the integrals in ( 7) and ( 9) can be evaluated closing the integral contour in the complex η-plane and using the residue theorem [43].Notably, the firing rate and the mean membrane potential depend only on the value of w at the pole of g(η) in the lower half η-plane: As a result, we only need to evaluate (10) at η = η − i∆, and thereby obtain a system of FREs composed of two ordinary differential equations (12b) This nonlinear system describes the macroscopic dynamics of the population of QIF neurons in terms of the population firing rate r and mean membrane potential v.
It is enlightening to compare mean-field model ( 12) with the equations of the spiking neurons.Note that Eq. (12b) resembles equations ( 1) and ( 2) for the individual QIF neuron, but without spike resetting.Indeed, the macroscopic firingrate variable r enters as a negative feedback term in Eq. (12b), and impedes the explosive growth of the mean membrane potential v.
This negative feedback, combined with the coupling term on the right hand side of (12a), describe the effective interaction between the firing rate and mean membrane potential at the network level.Therefore, the FREs (12) describe the effect of the single-cell spike generation and reset mechanism at the network level.
In the following we examine the dynamics described by Eq. ( 12), and show that they fully reproduce the macroscopic dynamics of the network of QIF neurons, even during episodes of strong spike synchrony.

C. Analysis of the firing-rate equations
To begin with the analysis of Eq. ( 12), we first note that in the absence of forcing, I(t) = 0, the only attractors of ( 12) are fixed points.Figure 1(a) shows a phase diagram of the system as a function of the mean external drive η and synaptic weight J, both normalized by the width of the input distribution [44].There are three qualitatively distinct regions of the phase diagram: 1 -A single stable node corresponding to a low-activity state (white), 2 -A single stable focus (spiral) generally corresponding to a high-activity state (gray), and 3 -A region of bistability between low and high firing rate (cyan; see a phase portrait of this region in Fig. 1(d)).Comparison of a sample bifurcation diagram of the fixed points from numerical simulation of networks of QIF neurons with that obtained from the FREs (12) shows an excellent correspondence, see Fig. 1(b,c).
A similar phase diagram can be readily reproduced by traditional heuristic firing-rate models, with one significant qualitative difference: the presence of a stable focus -and hence damped oscillations.Specifically, in the gray region of the phase diagram in Fig. 1(a), the system undergoes oscillatory decay to the stable fixed point.This oscillatory decay occurs as well for the high-activity state over a large extent of the region of bistability (cyan), see e.g.Fig. 1(d).
The presence of damped oscillations at the macroscopic level reflects the transitory synchronous firing of a fraction of the neurons in the ensemble.This behavior is common in spiking neuron models with weak noise, and is not captured by traditional firing rate models (see e.g.[22]).

D. Analysis of the firing-rate equations: Non-stationary inputs
To show that the FREs (12) fully describe the macroscopic response of the population of QIF neurons to time-varying stimuli (up to finite-size effects), we consider two types of stimulus I(t): 1 -a step function and 2 -a sinusoidal forcing.In both cases we simulate the full system of QIF neurons and the FREs (12).
Figure 2 shows the system's response to the two different inputs.In both cases the system is initially (t < 0) in a bistable regime and set in the low activity state, with parameters corresponding to those of Fig. 1(d).Left panels of Fig. 2 show the response of the system to a step current, applied at t = 0.The applied current is such that the system abandons the bistable region -see Fig. 1(a)-and approaches the high activity state, which is a stable focus.This is clearly reflected in the time series r(t), v(t), where the rate equations (12) exactly predict the damped oscillations exhibited by the meanfield of the QIF neurons.The raster plot in panel (e) shows the presence of the oscillations, which is due to the transitory synchronous firing of a large fraction of neurons in the population.Finally, at t = 30, the current is removed and the system converges -again, showing damped oscillations-to the new location of the (focus) fixed point, which clearly coexists with the stable node where it was originally placed (t < 0).
The right hand panels of Fig. 2 show the response of the model to a periodic current, which drives the system from one side of the bistable region to the other.As a result, we observe periodic bursting behavior when the system visits the stable focus region of the phase diagram Fig. 1(a).
To further illustrate the potential of the FREs (12) to predict and investigate complex dynamics in ensembles of spiking neurons, we present here a simple situation where the system of QIF neurons exhibits macroscopic chaos.This is observed by increasing the frequency ω of the sinusoidal driving, so that the system cannot trivially follow the stable fixed point at each cycle of the applied current.
Figure 3(a) shows a phase diagram obtained using Eq. ( 12).The shaded regions indicate parameter values where the rate Finally, to illustrate this chaotic state at the microscopic level, Fig. 3(e) shows the raster plot for 300 randomly chosen neurons, corresponding to the time series in Fig. 3(d).The irregular firing of neurons in Fig. 3(e) has some underlying structure that may be visualized ordering the same set of neurons according to their intrinsic currents as η k < η k+1 , with k ∈ [1, 300], see Fig. 3(f).Clearly, the maxima of the firing rate coincide with the synchronous firings of clusters of neurons with similar η values.The size of these clusters is highly irregular in time, in concomitance with the chaotic behavior.

E. Validity of the Lorentzian ansatz
Thus far, we have shown that the LA (5) solves the continuity equation ( 4), and confirmed that these solutions agree with the numerical simulations of the original system of QIF neurons (1)-( 2).Here we further clarify why the LA holds for ensembles of QIF neurons.
Transforming the voltage of the QIF neuron into a phase via V j = tan(θ j /2), the system (1)-( 2) becomes an ensemble of 'theta neurons' [40]: In the new phase variable θ ∈ [0, 2π), the LA (5) becomes: where the function α(η, t) is related with w(η, t) as Equations ( 5) and ( 14) are two representations of the socalled Poisson kernel on the half-plane and on the unit disk, respectively.These representations are related via equation (15), that establishes a conformal mapping from the half-plane Re(w) ≥ 0 onto the unit disk |α| ≤ 1.In the next section we show that variables r and v can be related, via the same conformal map, with the Kuramoto order parameter, which is a macroscopic measure of phase coherence [45,46].
The key observation supporting the applicability of the LA is the fact that equation ( 14) turns out to be the ansatz discovered in 2008 by Ott and Antonsen [37].According to the Ott-Antonsen (OA) theory, in the thermodynamic limit the dynamics of the class of systems generally converges to the OA manifold (14).In our case, for equation ( 13), we have Ω(η, t) = 1+η+Js+I and H(η, t) = i(−1 + η + Js + I).Thus far, the convergence of ( 16) to the OA manifold has been proven only for H(η, t) = H(t).This includes the well-known Kuramoto and Winfree models [47,48], but not system (13).However, there are theoretical arguments [49] that strongly suggest that the OA manifold is also attracting for η-dependent H, as numerically confirmed by a number of recent papers using theta neurons [50][51][52] and other phase-oscillator models [53][54][55][56][57].

F. Firing rate and Kuramoto order parameter
There exists a mapping between the macroscopic variables r and v and the so-called Kuramoto order parameter Z. Equation (15) relates, for each value of η, the firing rate and the mean membrane potential, both contained in w, with the uniformity of the phase density, measured by α -note that α = 0 in ( 14) yields a perfectly uniform density.The Kuramoto order parameter is obtained by integrating α over the whole population as where we assumed that the density ρ is the OA ansatz Eq. ( 14).
The quantity Z can be seen as the center of mass of a population of phases distributed across the unit circle e iθ .For a Lorentzian distribution of currents, g(η), we may derive the exact formula (see Appendix B) relating the Kuramoto order parameter Z with the firing-rate quantity W ≡ πr + iv. Figure 4 illustrates this conformal mapping of the right half-plane (r ≥ 0) onto the unit disk (|Z| ≤ 1).

III. CONCLUSIONS
We have presented a method for deriving firing rate equations for a network of heterogeneous QIF neurons, which is exact in the thermodynamic limit.To our knowledge, the resulting system of ordinary differential equations (12) represents the first exact FREs of a network of spiking neurons.
We would like to emphasize that the derivation of the FREs does not rely on assumptions of weak coupling, separation of time-scales, averaging, or any other approximation.Rather, the assumptions underlying the validity of equation ( 12) are that i) the QIF neurons be all-to-all connected, ii) heterogeneity be quenched, and iii) inputs be distributed according to a Lorentzian distribution.
The last two assumptions may seem particularly restrictive.Concerning (iii), it must be stressed that for arbitrary distributions of quenched heterogeneity the LA (5) remains generally valid, and therefore equation ( 10) can be used to discern the stability of the network states (see Appendix C).Furthermore, in Appendix C we also show numerical simulations using uniform-and Gaussian-distributed inputs, that reveal very similar macroscopic dynamics in response to time varying inputs.Even relaxing assumption (ii), numerical simulations with identical QIF neurons driven by external independent Gaussian noise sources show qualitative agreement with the FREs (12) (see Appendix D).In sum, the choice of a quenched Lorentzian distribution g(η) is thus a mere mathematical convenience, whereas the insights gained from the resulting firing-rate description are valid more generally.
Our derivation represents a sharp departure from, and a major advance compared to previous studies in several regards.Firstly, in the past it has only been possible to calculate the approximate firing rate of networks of spiking neurons for stationary states, or for weak deviations of such states [9-12, 15, 20, 22].The equations which result from these derivations are difficult to solve, and typically require special numerical methods.In contrast, FREs ( 12) can be easily analyzed and simulated, and exactly reproduce the behavior of the spiking network far from any fixed point and for arbitrary external currents.Secondly, firing-rate descriptions traditionally assume that the activity of a population of neurons is equivalent to a set of uncorrelated stochastic processes with a given rate, see e.g.[13,14].However, in simulations of spiking networks it is well known that the response of the network to non-stationary inputs generally involves some degree of spike synchrony.Recent theoretical work has sought to improve on classical rate models, which act as a low-pass filter of inputs, by fitting them with linear filters extracted from the corresponding Fokker-Planck equation for the network [20,22].These filters tend to generate damped oscillations when the external noise level is not too high, reflecting the presence of spike synchrony [15,22].The FREs (12) also capture this phenomenon, and reveal that the underlying mechanism is an interplay between firing rate and subthreshold voltage.Furthermore, because these equations are exact, we can explore the full nonlinear response of the network, such as the generation of chaotic states of synchronous bursting as shown in Fig. 3.
Recently it has been shown that, not only Kuramoto-like models, but a much wider class of phase-models, evolve in the OA manifold defined by (14).Specifically, networks of pulse-coupled oscillators [48] and theta-neurons [50][51][52] allow for an exact, low-dimensional description in terms of the Kuramoto order parameter (17).Although these later works use a finite-width phase coupling function that differs from the synaptic coupling (3), the obtained low-dimensional description in terms of the Kuramoto order parameter is analogous to ours, but in a different space [58].Indeed, we showed that the OA ansatz ( 14) is related, via the nonlinear transformation of variables (15), to the LA ansatz (5).Remarkably, this transformation establishes an exact correspondence between the Kuramoto order parameter and a novel, biophysically maningful macroscopic observable which describes the firing rate and mean membrane potential of the neuronal network.Interestingly, the low-dimensional description in terms of firing rates seems to be a more natural description for networks of spiking neurons, compared to that in terms of the Kuramoto order parameter.The firing rate equations (12) take a surprisingly simple form in the LA manifold, what makes them a valuable tool to explore and and understand the mechanisms governing the macroscopic dynamics of neuronal networks.
Finally, since the OA ansatz is the asymptotic solution for systems of the form in equation ( 16), applying the change of variables V = tan(θ/2) automatically implies that the LA should hold for populations governed by: where A, B, and C are related with Ω and Notably, equation ( 19) defines a wide family of ensembles of QIF neurons.
Therefore, the LA is actually valid for more general networks beyond the one investigated here -η j in Eq. ( 19) may also be a vector containing several forms of disorder.As particularly relevant cases, in Appendix E we provide the FREs governing the dynamics of an excitatory network with both, heterogeneous inputs η and synaptic weights J, as well as pair of interacting excitatory and inhibitory populations of QIF neurons.In addition, according to Eq. ( 19), the LA is also valid if synapses are modeled as conductances, in which case the reversal potentials may be distributed as well.Moreover, the role of gap-junctions or synaptic kinetics may be considered within the same framework.
The time it takes for the membrane potential of a QIF neuron (with I j > 0) to reach infinity from a given positive value of the membrane potential is arctan( I j /V p )/ I j .For I j ≪ V p , this expression can be approximated as: In simulations we considered V p = −V r = 100, and used the previous approximation.Thus, the time for the neurons to reach infinity from V p is 1/V p ≈ 10 −2 .Additionally, the time taken from minus infinity to V r is 1/V r ≈ 10 −2 .Numerically, once a neuron's membrane potential V j satisfies V j ≥ V p the neuron is reset to −V j and held there for a refractory time 2/V j .The neurons produce a spike when V j reaches infinity, i.e. a time interval 1/V j after crossing V p .The exact time of the spike can not be evaluated exactly in numerical simulations, since dt is finite.However, simulations agree very well with the theory provided that 1/V p ≫ dt.
To evaluate the mean membrane potential v, the population average was computed discarding those neurons in the refractory period.The choice V p = −V r is not needed, but in this way the observed mean membrane potential v agrees with the theory, consistent with the definition in (8).The mean synaptic activation (3) was evaluated using the Heaviside step function a τ (t) = Θ(τ − t)/τ with τ = 10 −3 .The instantaneous firing rates were obtained binning time and counting spikes within a sliding time window of size δt = 2 × 10 −2 (Figs. 2, C2, C3, D1) and δt = 4 × 10 −2 (Fig. 3).ξ = J SN r 0 , two equations that permit to find the locus of the saddle-node bifurcations systematically for arbitrary distributions of currents:

Bistability region for uniform and Gaussian distributions
Equations (C5) and (C6) are particularly easy to solve for uniform distributions of currents: for |η − η| < γ 0 otherwise After defining the rescaled parameters η = η/γ and J = J/ √ γ, we find two branches of saddle-node bifurcations emanating from the cusp point at η = −1/3: These functions are plotted in Fig. S1(a).Numerical simulations using the network of QIF neurons confirm the correctness of these boundaries (see the square symbols in Fig. C1).For Gaussian distributions, the solutions of (C5) and (C6) can be numerically found.The results are shown in Fig. C1(b).Again, there is a perfect agreement between Eqs. (C5) and (C6) -derived assuming the Lorentzian Ansatz-and the numerical estimations obtained simulating the network of QIF neurons.

Excitatory networks with external periodic currents
Firing rate equations (12) predict the existence of a stable focus in the shaded region of Fig. 1.Trajectories attracted to this fixed point display oscillations in the firing rate and mean membrane potential due to the transient synchronous firing of the QIF neurons.
When an external periodic current is injected to all neurons in the network, the spiral dynamics around the fixed point is responsible for the bursting behavior observed in Fig. 2, as well as for the emergence of the macroscopic chaos shown in Fig. 3.Here we investigate whether similar phenomena occur when an external periodic current of the same frequency is injected in an excitatory network with either uniform or Gaussian-distributed currents.
We introduced a sinusoidal forcing I(t) = I 0 sin(ωt) observing a behavior qualitatively identical to the one reported in the main text with Lorentzian g(η).Under a low-frequency forcing the system displays periodic bursting, see Fig. C2, provided the parameters are set inside the bistable region (see the black triangles in Fig. C1).Furthermore, note that the range of firing rates and mean membrane potentials in Figure C2 are similar to those of Fig. 2 -for Lorentzian distributions of currents.
As for the simulations shown in Figure (3), we next increase the frequency of the injected current up to ω = π to investigate the emergence of macroscopic chaos in a network with Gaussian-distributed currents.Fig. C3 shows the emergence of an apparently chaotic state.The observed dynamics is similar to that of Fig. 3, which suggests it is chaotic.This type of chaos persists in the thermodynamic limit N → ∞.On top of this, highly-dimensional but weakly chaotic dynamics is probably also present due to 'residual' finite-size effects.

APPENDIX D: IDENTICAL QIF NEURONS DRIVEN BY INDEPENDENT NOISE TERMS
We compare now the results obtained above and in the main text for quenched heterogeneity g(η) with the results for identical neurons g(η) = δ(η − η) driven by noise.Specifically, the inputs currents are now taken as where ξ j are independent white noise terms with expected values ξ j (t) = 0, and ξ j (t)ξ k (t ′ ) = 2D δ jk δ(t − t ′ ).
In the thermodynamic limit, the density ρ(V, t) obeys the Fokker-Planck equation: The Lorentzian ansatz is not a solution of this equation, but we demonstrate here that the phenomena observed with quenched heterogeneity arise also with independent noise sources.This qualitative similarity at the macroscopic level between quenched Lorentzian heterogeneity and Gaussian noise has been noted in previous work [60].A subtle point in Eq. (D2) is that its nondimensionalization entails a different (compared to ∆) rescaling with D: Ṽ = V /D 1/3 , η = η/D 2/3 , J = J/D 1/3 , t = tD 1/3 (implying r = r/D 1/3 ).
Numerical simulations reveal the existence of a region of bistability, see Fig. D1 Additionally, in order to investigate and compare the dynamical behavior of the identical noise-driven neurons with that of the FREs (12), an external current of intensity I 0 = 3 is injected to all neurons at time t = 0, like in Fig. 2. In Fig. D1(b,c) the time series of the firing rate and the mean membrane potential clearly display damped oscillations after the injection of the current, confirming the existence of a stable focus, exactly as observed in the FREs (12).The existence of a stable focus reflects the presence of transient spike synchrony in the network, as seen by the raster plot in Fig. D1(d).Remarkably, the raster plot and the time series closely resemble those of Fig. 2. The resemblance is not just qualitative, but rather there is a near quantitative fit between the network of QIF neurons driven by Gaussian noise and the FREs (12), which were derived assuming quenched Lorentzian heterogeneity (orange curves).

APPENDIX E: MODEL GENERALIZATIONS
To illustrate the potential of the LA for investigating more sophisticated networks, here we provide the low-dimensional FREs corresponding to an excitatory network of QIF neurons with independently distributed currents and synaptic weights, and to two interacting populations of excitatory and inhibitory QIF neurons with distributed currents.

Excitatory population with heterogeneous currents and synaptic weights
As a first example, let us assume that both the currents η and the synaptic weights J are distributed -this type of heterogeneity was also considered in [42].The input currents read then [68] I j = η j + J j r(t) + I(t).

FIG. 2 .
FIG. 2. (color online).The transient dynamics of an ensemble of QIF model neurons (1)-(2) are exactly described by the FREs (12).The instantaneous firing rate (Panels a,b) and mean membrane potential (Panels c,d) of the QIF neurons and the FREs are depicted in black and orange, respectively.Panels (e,f) show the raster plots of 300 randomly selected QIF neurons of the N = 10 4 in the ensemble.Left Panels (a,c,e,g): At time t = 0, a current I = 3 is applied to all neurons, and set to zero again at t = 30; stimulus I(t) shown in panel (g).Right panels (b,d,f,h): At time t = 0 a sinusoidal current is applied to all neurons I(t) = I0 sin(ωt), with I0 = 3, ω = π/20; stimulus I(t) shown in panel (h).Parameters: J = 15, η = −5, ∆ = 1.

FIG. 3 .
FIG. 3. (color online).Firing rate model predicts the existence of macroscopic chaos in a ensemble QIF neurons (1) with an injected periodic current I(t) = I0 sin(ωt).(a) Phase diagram with the regions where chaos is found, obtained using the rate model(12).In the black-shaded region, there is only one chaotic attractor, whereas in the cyan region the chaotic attractor coexists with a periodic attractor.Dotted lines correspond to the bistability boundary of Fig.1(a), depicted to facilitate comparison.(b) Chaotic trajectory obtained simulating the FREs (12); the Lyapunov exponent is λ = 0.183 . ... (c) Chaotic trajectory obtained from the QIF neurons -parameters corresponding to the (yellow) triangle symbol in panel (a).(d) Time series for the rate model (orange) and QIF model (black).(e) Raster plot corresponding to 300 randomly selected neurons.(f) Same raster plot as in (e), ordering neurons according to their intrinsic current η k .Parameters: I0 = 3, ω = π; and η = −2.5,J = 10.5 (triangle symbol in (a)).
FIG. C1. (color online).Phase diagrams for (a) uniform and (b) Gaussian distributions of currents g(η).The cyan-shaded region represents the bistable region, with the solid lines corresponding to saddle-node bifurcations analytically obtained from Eqs. (C5) and (C6).Square symbols are estimations of the bifurcations loci obtained from by direct numerical simulations of N = 10 4 QIF neurons.Triangle symbols indicate parameter values used in numerical simulations of Figs.C2 and C3.
(a), analogous to the one observed for quenched heterogeneity, cf.Figs.1(a) and C1.Obviously, true bistability only holds in the thermodynamic limit, while what we observe are rather exceedingly long residence times close to each fixed point (see figure caption for details).