Asynchronous Rate Chaos in Spiking Neuronal Circuits

The brain exhibits temporally complex patterns of activity with features similar to those of chaotic systems. Theoretical studies over the last twenty years have described various computational advantages for such regimes in neuronal systems. Nevertheless, it still remains unclear whether chaos requires specific cellular properties or network architectures, or whether it is a generic property of neuronal circuits. We investigate the dynamics of networks of excitatory-inhibitory (EI) spiking neurons with random sparse connectivity operating in the regime of balance of excitation and inhibition. Combining Dynamical Mean-Field Theory with numerical simulations, we show that chaotic, asynchronous firing rate fluctuations emerge generically for sufficiently strong synapses. Two different mechanisms can lead to these chaotic fluctuations. One mechanism relies on slow I-I inhibition which gives rise to slow subthreshold voltage and rate fluctuations. The decorrelation time of these fluctuations is proportional to the time constant of the inhibition. The second mechanism relies on the recurrent E-I-E feedback loop. It requires slow excitation but the inhibition can be fast. In the corresponding dynamical regime all neurons exhibit rate fluctuations on the time scale of the excitation. Another feature of this regime is that the population-averaged firing rate is substantially smaller in the excitatory population than in the inhibitory population. This is not necessarily the case in the I-I mechanism. Finally, we discuss the neurophysiological and computational significance of our results.


Author Summary
Cortical circuits exhibit complex temporal patterns of spiking and are exquisitely sensitive to small perturbations in their ongoing activity. These features are all suggestive of an underlying chaotic dynamics. Theoretical works have indicated that a rich dynamical reservoir can endow neuronal circuits with remarkable computational capabilities. Nevertheless, the mechanisms underlying chaos in circuits of spiking neurons remain unknown. We combine analytical calculations and numerical simulations to investigate this fundamental issue. Our key result is that chaotic firing rate fluctuations on the time scales of the synaptic dynamics emerge generically from the network collective dynamics. Our results Introduction single neuron properties? How do excitation and inhibition contribute to the emergence of these states? To what extent these chaotic dynamics share similarities with those exhibited by the SCS model? We first study these questions in one population of inhibitory neurons receiving feedforward excitation. We then address them in networks of two populations, one inhibitory and the other excitatory, connected by a recurrent feedback loop. A major portion of the results presented here constitutes the core of the Ph.D thesis of one of the authors (O.H) [27].

One population of inhibitory neurons: General theory
We consider N randomly connected inhibitory spiking neurons receiving an homogeneous and constant input, I. The voltage of each neuron has nonlinear dynamics, as e.g. in the leaky integrate-and fire (LIF model, see Materials and Methods) or in conductance-based models [20].
The connection between two neurons is J ij = JC ij (i, j = 1, 2. . .N), with J 0, and C ij = 1 with probability K/N and 0 otherwise. The outgoing synapses of neuron j obey t syn dS j ðtÞ dt where S j (t) is the synaptic current at time t and τ syn the synaptic time constant. When neuron j fires a spike (time t s j ), S j increments by J. Thus, the total input to neuron i, h i (t) = I+∑ j J ij S j (t), satisfies: We assume K ) 1, hence the number of recurrent inputs per neuron is K AE Oð ffiffiffi ffi K p Þ. Scaling J and I as: J ¼ ÀJ 0 = ffiffiffi ffi K p , I ¼ ffiffiffi ffi K p I 0 , the time-averaged synaptic inputs are Oð ffiffiffi ffi K p Þ and their spatial (quenched) and temporal fluctuations are O(1) [28,29]. Finite neuronal activity requires that excitation and inhibition cancel to the leading order in K. In this balanced state, the mean and the fluctuations of the net inputs are O(1) [28,29]. The properties of the balanced state are well understood if the synapses are much faster than all the typical time constants of the intrinsic neuronal dynamics [30]. Temporally irregular asynchronous firing of spikes is a hallmark of this regime [13,28,29,31,32]. However, this stochasticity does not always correspond to a true chaotic state [28,29,[33][34][35][36]. In fact, this depends on the spike initiation dynamics of the neurons [37]. The opposite situation, in which some of the synapses are slower than the single neuron dynamics, remains poorly understood. This paper mostly focuses on that situation.
When the synaptic dynamics is sufficiently slow compared to the single neuron dynamics, the network dynamics can be reduced to the set of non-linear first order differential equations: t syn dh i ðtÞ dt ¼ Àh i ðtÞ þ I þ X j J ij r j ðtÞ ð3Þ where r i (t) is the instantaneous firing rate of neuron i and g(h) is the neuronal input-output transfer function [20]. These are the equations of a rate model [20,38] in which the activity variables correspond to the net synaptic inputs in the neurons. Eqs (3)-(4) differ from those of the SCS model in that they have a well defined interpretation in terms of spiking dynamics, the time constant has a well defined physiological meaning, namely, the synaptic time constant, the transfer function quantifies the spiking response of the neurons and is thus positive, the interactions satisfy Dale's law and the neuronal connectivity is partial. Dynamical mean-field theory (DMFT). We build on a DMFT [19] to investigate the dynamics, Eqs (3)-(4), in the limit 1 ( K ( N. Applying this approach, we rewrite the last two terms in the right hand side of Eq (3) as a Gaussian noise whose statistics need to be self-consistent with the dynamics. This yields a set of self-consistency conditions which determine the statistics of the fluctuations, from which the synaptic net inputs and the firing rates of the neurons can be calculated. This approach is described in detail in the Materials and Methods section.
The DMFT shows that, for a given transfer function, depending on the parameters J 0 and I 0 , the dynamics either converge to a fixed point state or remain in an asynchronous, time-dependent state. In the fixed point state, the net inputs to the neurons, h 0 i , (i = 1. . .N) are constant. Their distribution across the population is Gaussian with mean μ and variance J 0 2 q. The DMFT yields equations for μ, q, as well as for the distribution of firing rates r 0 i (i = 1. . .N) (Eqs (24)- (25) and (36)). In the time-dependent state, h i (t) exhibit Gaussian temporal fluctuations, which are characterized by a mean, μ = [hh(t)i], and a population-averaged autocovariance (PAC) function, σ(τ) = [hh(t)h(t+τ)i] − μ 2 ([Á] and h Á i denote means over the population and over time, respectively). Solving the set of self-consistent equations which determine σ(τ) and μ (Eqs (25), (27) and (37)- (38)) indicates that σ(τ) decreases monotonically along the flow of the deterministic dynamics, thus suggesting that the latter are chaotic. To confirm that this is indeed the case one has to calculate the maximum Lyapunov exponent of the dynamics (which characterizes the sensitivity of the dynamics to initial conditions [39]) and verify that it is positive. This can be performed analytically in the framework of DMFT [19]. However, this is beyond the scope of the present paper. Therefore, in the specific examples analyzed below we rely on numerical simulations to verify the chaoticity of the dynamics.
For sufficiently small J 0 , the fixed point state is the only solution of the dynamics. When J 0 increases beyond some critical value, J c , the chaotic solution appears. We show in the Materials and Methods section that J c is given by: where q and μ are computed at the fixed point state and Dz ¼ e À z 2 2 ffiffiffi ffi 2p p dz. On the stability of the fixed point state. The NxN matrix characterizing the stability of the fixed point is D ¼ M ffiffi ffi N p À I with I the NxN identity matrix and: where h 0 j is the total input in neuron j at the fixed point. This is a sparse random matrix with, on average, K non zero elements per line or column. In the limit N ! 1, these elements are uncorrelated, have a mean ÀJ 0 ffiffi ffi (for large N, the second moment of the matrix elements is equal to their variance). Interestingly, Eq (5) means that the SD of the elements of M crosses 1 (from below) at J c . As J 0 increases, the fixed point becomes unstable when the real part of one of the eigenvalues crosses 1. Note that that for large K, D always has a negative eigenvalue, which is Oð ffiffiffi ffi K p Þ.
In the specific examples we investigate below, simulations show that when the chaotic state appears the fixed point becomes unstable. This implies that for J < J c given by Eq (5) the real parts of all the eigenvalues of M ffiffi ffi N p are smaller than 1 and that for J = J c , the real part of one of the eigenvalues, the eigenvalue with maximum real part, crosses 1. This suggests the more general conjecture that in the limit 1 ( K ( N the eigenvalue with the largest real part of M= ffiffiffiffi N p is: Below we compare this formula to results from numerical diagonalization of M= ffiffiffiffi N p .

One population of inhibitory neurons: Examples
The above considerations show that when synapses are slow, the dynamics of inhibitory networks is completely determined by the transfer function of the neurons. Therefore, to gain insights into the way dynamics become chaotic in such systems we proceed by investigating various spiking models that differ in the shape of their transfer functions.
Sigmoidal transfer functions. Neurons in a strong noise environment can be active even if their net inputs are on average far below their noiseless threshold, whereas when these inputs are large the activity saturates. The transfer functions of the neurons can therefore be well approximated by a sigmoid. Specifically here we consider the dynamics described in Eqs (3)-(4) with a sigmoidal transfer function: This form of the sigmoid function makes analytical calculations more tractable. Fig 1A shows that for J 0 = 4, I 0 = 1, the simulated network dynamics converge to a fixed point. This is not the case for J 0 = 6 and J 0 = 15 (Fig 1B, 1C). In these cases the activities of the neurons keep fluctuating at large time. Note also that the mean level of activity is different for the three neurons. This is a consequence of the heterogeneities in the number of inputs the neurons receive.
These differences in the network dynamics for these three values of J 0 are consistent with the full phase diagram of the DMFT in the parameter space I 0 − J 0 . Fig 2A depicts the results obtained by solving numerically the self-consistent equations that define chaos onset with g(x) = ϕ(x) (Eqs (17)- (18) in S2 Text). In the region above the line a chaotic solution exists whereas it does not exist below it. Simulations indicate that in the region above the line, very small perturbations from the fixed point state drive the network toward the time dependent state. In other words, the fixed point solution is unstable above the line: the bifurcation to the time dependent state is thus supercritical.
The instability of the fixed point on this line is also confirmed by direct diagonalization of the matrix M= ffiffiffiffi N p (see Eq (6)). To this end, we solved numerically the mean field equations for different values of J 0 to obtain μ and q, randomly sampled h 0 i values from the distribution defined by μ and q to generate the random matrix matrix M= ffiffiffiffi N p , and then computed numerically the spectrum of the matrix (for N = 10000). Examples of the results are plotted in Fig 3A for two values of J 0 , one below and one above the critical value J c . In both cases, the bulk of the spectrum is homogeneously distributed in the disk of radius λ max centered at the origin. Fig 3B plots λ max computed numerically (dots) and compare the results to our conjecture, Eq (7) (solid line). The agreement is excellent. The instability of the fixed point corresponds to λ max crossing 1.
To verify the chaoticity of the time dependent state predicted by the DMFT in the region above the bifurcation line we simulated the dynamics and computed numerically the largest Lyapunov exponent, Λ, for different values of I 0 and J 0 (see Materials and Methods for details). The results plotted in Fig 2A (red dots and inset) show that Λ crosses zero near the DMFT bifurcation line and is positive above it. Therefore the dynamics observed in simulations are chaotic in the parameter region above this line as predicted by the DMFT. We solved numerically the parametric self-consistent differential equation which determined the PAC, σ(τ), (Eqs (25), (29) and (37)- (38)) for different values of J 0 and I 0 . An example of the results is plotted in Fig 2B. It shows that numerical simulations and DMFT predictions are in very good agreement. Moreover, simulations with increasing values of N and K indicate that the small deviations from the DMFT predictions are due to finite N and K effects; a detailed study of these effects is reported in S1 Text. Fig 4A shows the bifurcation diagram of the PAC amplitude, σ 0 − σ 1 . For J 0 below the bifurcation point (BP) the PAC amplitude is zero, which corresponds to the fixed point state (solid blue line). At the bifurcation the fixed point loses stability (dashed blue line) and a chaotic state with a strictly positive PAC amplitude emerges (black line).
We studied analytically the critical behavior of the dynamics at the onset of chaos. We solved perturbatively the DMFT equations for 0 < δ = J 0 − J c ( 1, as outlined in the Materials and Methods section and in S2 Text. This yields (σ(τ) − σ 1 ) / δ α /cosh 2 (τ/τ dec ), with α = 1 and a decorrelation time scaling like τ dec / δ β with β = −1/2. Therefore at the onset of chaos, the PAC amplitude vanishes and the decorrelation time diverges. We show in the Materials and Methods section that this critical behavior with exponents α = 1, β = −1/2, is in fact a general property of the model, Eqs (3)-(4), whenever g(h) is twice differentiable. It should be noted that in the SCS model the PAC also vanishes linearly at chaos onset. However, the critical exponent of the decorrelation time is different (β = −1) [19].
The inset in Fig 4A compares the PAC amplitude obtained by numerically solving Eq (27) (black line) with the corresponding perturbative result (red line) for small δ. The agreement is excellent. In fact, the perturbative calculation provides a good estimate of the PAC even if δ is as large as 0.2J c (Fig 4A, main panel and Fig 4B). More generally, the PAC can be well fitted with the function (σ 0 − σ 1 ) Á cosh − 2 (τ/τ dec ) (Fig 4C, inset) providing an estimate of the decorrelation time, τ dec , for all values of J 0 . Fig 4C plots τ dec vs. σ 0 − σ 1 for I 0 = 1. It shows that the formula t dec / 1= ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi s 0 À s 1 p we derived perturbatively for small δ provides a good approximation of the relationship between the PAC amplitude and the decorrelation time even far above the bifurcation. for inhibitory population rate model with g(x) = ϕ(x). The matrix was diagonalized numerically for N = 10000, K = 400, I 0 = 1 and different values of J 0 . A: The bulk of the spectrum for J 0 = 6 (blue) and for J 0 = 1.12 (red). Left: The imaginary parts of the eigenvalues are plotted vs. their real parts for one realization of M. This indicates that the support of the spectrum is a disk of radius λ max . Right: Histograms of N eig /R (one realization of M) where N eig is the number of eigenvalues with a modulus between R and R+ΔR (ΔR = 0.0428 (top), 0.0093 (bottom)) for J 0 = 6 (top) and J 0 = 1.12 (bottom). The distribution of eigenvalues is uniform throughout the spectrum support. B: The largest real part of the eigenvalues (black dots), λ max , is compared with the conjecture, Eq (7) (solid line). The fixed point loses stability when λ max crosses 1. Threshold power-law transfer function. We next consider the dynamics of the network (Eqs (3)-(4)) with a transfer function where γ > 0 and H(x) = 1 for x > 0 and 0 otherwise. Non-leaky integrate-and-fire neurons [40] (see also S3 Text) and θ-neurons [41][42][43][44] correspond to γ = 1 and γ = 1/2, respectively. The transfer functions of cortical neurons in-vivo can be well fitted by a power-law transfer function with an exponent γ % 2 [45,46]. 4 as we show analytically in S4 Text. For γ < 1/2, the integral in the right hand side of Eq (5) diverges. Equivalently, the elements of the stability matrix have infinite variance. Therefore, the DMFT predicts a chaotic dynamics as soon as J 0 > 0.
To compare these predictions with numerical simulations, we simulated different realizations of the network (N = 32000, K = 400, I 0 = 1) for various values of J 0 . For each value of J 0 and γ we determined whether the dynamics converge to a fixed point or to a time dependent state as explained in the Materials and Methods section. This allowed us to compute the fraction of networks for which the dynamics converge to a fixed point. The solid red line plotted in Fig 5B corresponds to a fraction of 50% whereas the dotted red lines correspond to fractions of 5% (upper line) and 95% (lower line). We also estimated the Lyapunov exponent, Λ, for each values of J 0 and γ. The blue line in Fig 5B corresponds to the location where Λ changes sign according to our estimates (see Materials and Methods for details).  (11) in S2 Text) in the limit δ ! 0. C: Blue dots: Decorrelation time, τ dec vs. PAC amplitude. The PAC, σ(τ) − σ 1 , was obtained by solving numerically the DMFT equations and τ dec was estimated by fitting the result to the function A/cosh 2 (τ/τ dec ). Red: In the whole range, J 0 2 [5,7] considered, τ dec can be well approximated by t dec ¼ 4:97= ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi s 0 À s 1 p . This relation becomes exact in the limit σ 0 − σ 1 ! 0. Inset: Numerical solution of the DMFT equations for J 0 = 6.65 (blue dots) and the fit to A/cosh 2 (τ/τ dec ) (red). The fit is very good although this is far from bifurcation. For γ > % 0.6, the fraction of networks with an unstable fixed point varies sharply from 0 to 100% in the vicinity of the bifurcation line predicted by the DMFT. Moreover, for these values of γ, the spectrum of the matrix M= ffiffiffiffi N p is homogeneously distributed in the disk of radius λ max centered at the origin and the values of λ max agrees with Eq (7). This is shown in Fig  However, as γ ! (1/2) + , the discrepancies between DMFT and simulations become more pronounced. Very close to γ = (1/2) + there is a whole range of values of J 0 for which the DMFT predicts chaos whereas in numerical simulations the dynamics always converge to a fixed point. This discrepancy can be understood by observing that the integral over the Gaussian measure in Eq (5) corresponds to a population average over neurons. When γ ! (1/2) + , the region where z is just above À m J 0 ffi ffi q p dominates the integral; in other words, the neurons with positive close-to-threshold net inputs are those that make the largest contribution to the destabilization of the fixed point. On the other hand, the DMFT shows that these neurons become extremely rare as γ ! (1/2) + : in that limit μ c increases sharply, thus shifting the center of the Gaussian distribution to very large positive values. Therefore, we would need to simulate outrageously large networks to obtain a quantitative agreement with the DMFT predictions for the locations of the bifurcation to chaos. Similar arguments explain why when γ < 1/2 we find a transition from fixed point to chaos in numerical simulations for J 0 < % 0.9 although according to the DMFT the fixed point is always unstable since the integral in Eq (5) diverges.
Numerical diagonalization of M ffiffi ffi N p shows that when γ < % 0.6 (i) the eigevalues in the bulk of the spectrum are distributed in a disk centered at the origin and that this distribution is less and less homogeneous as γ ! (1/2) + (ii) the eigenvalue λ max governing the instability exhibits substantial deviations from Eq (7) especially for large J 0 ( Fig 6C) (iii) λ max exhibits large sample to sample fluctuations (results not shown). We conjecture that these features are due to large finite N and K effects and stem from the fact that the SD of the elements of M ffiffi ffi N p diverges when γ ! (1/2) + . We studied the dynamics in detail for γ = 1. The DMFT predicts that J c ¼ ffiffi ffi 2 p for all I 0 and K (K large). As already mentioned, the simulations agree well with this result (Fig 5B). We studied analytically the dynamics for J 0 close to this transition (Fig 7A-7C). To this end, we solved the self-consistent DMFT equations in the limit δ = J 0 − J c ! 0 + . The perturbative calculation, explained in S4 Text, is less straightforward than in the case of a sigmoid transfer function. This stems from the fact that at the threshold, the threshold-linear transfer function is only differentiable once. It yields that σ − σ 1 * δ α σ s (τ/δ β ) with α = 2, β = −1/2 and the function σ s (x)) has to be determined numerically. The function σ s is plotted in Fig 7B. It can be well fitted to the function A[cosh(x/x dec )] −1 with A = 12.11 and x dec = 2.84 (see Fig 7B, inset). In particular, for small δ, the amplitude and the decorrelation time of the PAC are related by τ dec / 1/(σ 0 − σ 1 ) 1/4 . Note that the amplitude of the PAC vanishes more rapidly (α = 2) than for sigmoidal transfer functions (α = 1) whereas the decorrelation time diverges with the same critical exponent (β = −1/2) in the two cases.  (27). Unlike what we found for the sigmoid transfer function, δ must be very small (δ < % 0.03J c ) to achieve a good quantitative agreement. It should be noted, however, that the quality of the fit of σ − σ 1 to A[cosh(τ/τ dec )] −1 does not deteriorate by much even far from the bifurcation ( The agreement is reasonably good but not perfect. We show in S1 Text that the discrepancy between the two improves as the network size increases but that finite size effects are stronger here than in the rate model with sigmoid transfer function.
Leaky integrate-and-fire (LIF) inhibitory networks. Our objective here is to obtain further insights into the relevance of the chaotic behavior exhibited by rate dynamics, Eqs (3)-(4), Inset: The function (σ(τ) − σ 1 )/δ 2 (black) can be well fitted to A/cosh(x/x dec ) (red dots, A = 12.11, x dec = 2.84). C: Decorrelation time, τ dec vs. PAC amplitude (blue). The function σ(τ) − σ 1 was obtained by integrating numerically Eq (29) and τ dec was estimated by fitting this function to A/cosh(τ/τ dec ). Red: In the whole range of J 0 considered (J 0 2 [1.4, 1.9] the relation between τ dec and σ 0 − σ 1 can be well approximated by y ¼ 5:29= ffiffi ffi to understand spiking network dynamics. The dynamics of one population of LIF spiking neurons reduces to Eqs (3)-(4) with the transfer function in the limit where the synapses are much slower than the cell membrane time constant, τ m . Our goal is twofold: 1) to study the emergence of chaos in this rate LIF rate model and 2) to compare it to full spiking dynamics and characterize the range of values of the synaptic time constant for which the two dynamics are quantitatively or at least qualitatively similar. Figs 8, 9 depict typical patterns of neuronal activity in simulations of the inhibitory spiking LIF model. For strong and fast synapses (τ syn = 3 ms, Fig 8A), neurons fire spikes irregularly and asynchronously (Fig 8A). Fig 8B shows that when τ syn = 100 ms the population average firing rate remains essentially the same (* 14.1 Hz) and the network state stays asynchronous. The spiking patterns, however, change dramatically: with large τ syn neurons fire irregular bursts driven by slowly decorrelating input fluctuations (Fig 9A, blue). Fig 9B shows that reducing J 0 increases the firing rate, reduces the amplitude of the fluctuations (Fig 9B, inset) and slows down their temporal decorrelation. Eventually, for small enough J 0 , σ(τ) becomes flat and the fluctuations are negligible. Fig 10 compares the dynamics of the rate to those of the spiking LIF networks. Panels A,B show that for J 0 = 2, I 0 = 0.3 and τ syn = 100 ms, σ(τ), the distributions of the time averages of neuronal firing rates and net inputs, hr i i and hh i i, are essentially the same in the simulations of the two networks. When reducing τ syn down to τ syn > % 15 ms, the function σ(τ/τ syn ) measured in the spiking network simulations, changes only slightly. In fact, this function is remarkably similar to what is found for the corresponding function in the DMFT and in simulations of the LIF rate model (Fig 11A). Fitting σ(τ) with the function B+A[cosh(τ/τ dec )] −1 yields τ dec % 2.45Áτ syn .
How small can τ syn be for the two models to still behave in a quantitatively similar manner? Simulations show that this value increases with the mean activity of the network (see examples in Fig 11) but that for reasonable firing rates, fewer than several several tens of Hz, the fluctuations have similar properties in the two models even for τ syn % 20 ms.  The results for the rate model are also plotted (black). The firing rates are * 15 Hz in A and C, * 10 Hz in B and * 30 Hz in D, in good agreement with the prediction from the balance condition ([hri] = 100I 0 /J 0 Hz). As the population firing rate increases, a larger τ syn is needed for good agreement between the spiking and the rate model.  We conducted extensive numerical simulations of the inhibitory LIF rate and spiking models (N = 40000, K = 800) to compute their phase diagrams in the I 0 − J 0 parameter space. The results for the rate model are plotted in Fig 12. For sufficiently small J 0 the dynamics always converge to a fixed point whereas for sufficiently large J 0 the network always settles in a state in which the activity of the neurons keeps fluctuating at large time. We show in S5 Text that in this regime the maximum Lyapunov exponent is strictly positive, therefore the dynamics are chaotic. Between these two regimes, whether the dynamics converge to a fixed point or to a chaotic state depends on the specific realization of the connectivity matrix. The fraction of networks for which the convergence is to a fixed point depends on J 0 . The range of J 0 where this fraction varies from 95% to 5% is notably large as shown in Fig 12. Additional simulation results on this issue are depicted in S5 Text. The counterpart of this behavior in the spiking network is that when J 0 is small, neurons fire regular spikes tonically whereas for sufficiently large J 0 they fire highly irregular bursts. The transition between the two regimes occurs for similar values of J 0 in the rate and in the spiking networks. In both networks this transition is driven by the neurons with low firing rates; i.e., with larger numbers of recurrent inputs. These neurons are the first to become bursty as J 0 increases (see S6 Text).
In Fig 13A we plot the bifurcation diagram of the model as obtained in the numerical solution of the DMFT equations (black line) and as calculated in simulations of the rate model (blue dots) and of the spiking network with τ syn = 25 ms (red ×'s) and τ syn = 7.5 ms (green ×'s). The rate model simulations are in good agreement with DMFT for 0.8 < % J 0 < % 2. For larger J 0 the discrepancy becomes significant and increases with J 0 . This is because of finite K effects that grow stronger as J 0 increases as shown in the right inset in Fig 13A, for J 0 = 3 (blue) and J 0 = 4 (red). Fig 13A also shows that, as discussed above, the amplitude of the PACs obtained in simulations of the LIF rate and spiking networks are barely different provided the synaptic time constant is sufficiently large. Finally, according to the DMFT the fixed point should be always unstable since for the LIF transfer function the elements of the stability matrix always have an infinite variance or, equivalently, the integral in Eq (5) always diverges. This can be seen in the close-up in the left inset of Fig 13A, indicating that the PAC amplitude is non-zero for small J 0 and that it approaches 0 very slowly as J 0 decreases. By contrast, in numerical simulations in the same range of J 0 , the dynamics are not chaotic for most of the realizations of the network: they converge to a fixed point, as shown in Fig 12. The explanation for this difference is as for the rate model with threshold power-law transfer function with γ < 1/2 (see above).

Two asynchronous chaos mechanisms in excitatory-inhibitory recurrent networks
We now consider EI spiking networks with recurrent feedback interactions between the two populations. The synaptic strengths and time constants are J ab 0 = ffiffiffi ffi K p and τ αβ (α, β 2 {E, I}). Assuming slow synapses, the dynamics can be reduced to four sets of equations for the four types of synaptic inputs, h i ab ðtÞ (Materials and Methods, Eq (17)). The DMFT yields self-consistent equations for the statistics of these inputs. These equations can be analyzed straightforwardly for the fixed point state. In contrast to purely inhibitory networks where the fixed point loses stability only via a bifurcation to chaos, it can now also lose stability via a Hopf bifurcation. This depends on the synaptic time constants. When this happens the network develops synchronous oscillations which break the balance of excitation and inhibition (the oscillation amplitude diverges for large K).
We focus here on instabilities which lead to chaos. Their locations in the 6 dimensional parameter space (4 synaptic strengths, 2 external inputs) of the model can be derived for a general transfer function (Eqs (54)-(55)). Differential equations for the PAC functions, σ αβ (τ), can also be derived in the chaotic regime. However, general analytical characterization of their solutions is particularly difficult. Leaving such study for future work, we mostly focus below on numerical simulations. Our key result is that in EI networks asynchronous chaos emerges in two ways, one driven by I-I interactions (II mechanism) and the other by the EIE loop (EIE mechanism).
EI network with threshold-linear transfer function. We first study a EI network in which all the neuronal transfer functions are threshold-linear. Fig 14 plots for different K the phase diagram of the DMFT of this model in the J IE 0 À J II 0 parameter space, when J EE 0 ¼ 0 and I E = I I = 1, J EI 0 ¼ 0:8.(The phase-diagram for a non-zero of value J EE 0 , J EE 0 ¼ 1:5, is plotted and briefly discussed in S7 Text). On the lines, the dynamics bifurcate from fixed point (below the lines) to chaos (above). As J II 0 decreases the lines go to infinity. Numerical simulations indicate the maximum Lyapunov exponent changes sign very close to these lines (compare red line and red dots) in good agreement with DMFT. For any finite K, the instability line exhibits a reentrance, crossing the J II 0 -axis at J II 0 ¼ ffiffi ffi 2 p , where the instability occurs in a purely inhibitory network; in this respect, the limit K ! 1 is singular. Solving the self-consistent equations for the average firing rates, r E and r I , one finds that the two populations can have a comparable firing rate for large J II 0 when J IE 0 is not too large. As J II 0 becomes small, the activity in the E population becomes much lower than in the I population. In fact, for K ! 1, r E vanishes on the line J II 0 ¼ I I I E J EI 0 ¼ 0:8 and is zero for J II 0 < I I I E J EI 0 (white region in Fig 14). In other words, in the latter case, inhibition is not balanced by excitation in the E population.
As shown above, in the single inhibitory population case with threshold-linear transfer functions the transition to chaos occurs at J 0 ¼ ffiffi ffi is small, but when J II 0 is large this happens only for h IE . By contrast, dividing τ II by 10 (purple) has very little effect when J II 0 is small but the fluctuations of all the inputs are substantially faster when J II 0 is large. Fig 15 also demonstrates the effect of changing τ αβ on the PAC of the net inputs to the E neurons, h i E ðtÞ ¼ I E þ h i EE ðtÞ À h i EI ðtÞ (corresponding results for the I population are shown in S8 Text). The PAC in the reference case is plotted in gray. For large J II 0 , a ten-fold increase in τ II causes the PAC width to become ten times larger and the PAC amplitude increases (Fig 15A, blue; see also inset). For a ten-fold decrease in τ II (purple) compared to reference, the width of the PAC is smaller but by a smaller factor whereas its amplitude is greatly reduced. By contrast, a ten-fold increase in τ IE has no noticeable effect, either on the width or on the amplitude of the PAC (black). Fig 15B plots the PAC of the total input to the E population for small J II 0 . Here, decreasing τ II by a factor of 10 (purple line) weakly affects the width as well as the amplitude of the PAC. In contrast, a ten-fold increase of τ IE (black) widens the PAC by a comparable factor (see also inset). A similar widening occurs if τ EI is increased ten-fold (see S8 Text).
This phenomenology can be understood as follows. In the large J II 0 regime, the II interactions play the key role in the generation of chaos. Therefore, the time scale of the fluctuations in the activity of the I neurons is essentially determined by τ II . Thus if the latter is 10 times larger than reference, the I inputs to the E neurons are slowed down by the same factor. At the same time, the filtering effect of the EI synapses becomes weaker and thus the amplitude of the PAC of the net input in the E neurons increases. The effect of decreasing τ II stems from the filtering effect of the EI synapses which is now stronger than in the reference case. Finally, changing τ IE has no noticeable effect since the fluctuations are generated by the II interactions. By contrast, when J II 0 is small, II interactions are not sufficient to generate chaotic fluctuations. In this regime, the EIE loop drives these fluctuations if J IE 0 is sufficiently large. That is why the time scale of the activity fluctuations depends primarily on τ IE and to a much smaller extent on τ II .
These results point to the existence of two mechanisms for chaos emergence in two population networks; they differ by the type of the dominant interactions (EIE or II) and therefore on the synaptic time constants which settle the time scale of the activity fluctuations. Another difference is that in the EIE mechanism, the E population is always significantly less active than the I population. This is not the case in the II mechanism.
Two-population spiking LIF network. We ran a similar analysis for LIF networks. Fig 16A, 16C plot the PACs of h i E ðtÞ for the LIF spiking and rate models (PACs of h i I ðtÞ are shown in S9 Text). In all panels J EE 0 ¼ 0, J IE 0 ¼ 3, J EI 0 ¼ 0:8 and τ EI = 3 ms. For J II 0 ¼ 4 ( Fig  16A), increasing τ II slows down the fluctuations. By contrast, changing τ IE only has a very mild effect (S10 Text). This is because the fluctuations are essentially driven by the II interactions. For τ II > % 15 ms, the fluctuation statistics are quantitatively similar in the spiking and the rate models: in both, the decorrelation time, τ dec % 2τ II (Fig 16A, inset). Moreover, simulations indicate that the dynamics of the rate model are chaotic (Λ % 1.7/τ II ). The trace in Fig 16B shows that with large τ II (=100 ms) the spiking pattern is bursty. The membrane potential between bursts exhibit slow fluctuations because they are generated by the slow II connections. Fig 16C plots the PACs of h i E ðtÞ for J II 0 ¼ 1. Here also, the LIF rate model operates in a chaotic regime (Λ % 120s −1 ). In the spiking model the PACs exhibit a slow time scale but also a fast one (the sharp peak around τ = 0). These correspond to the slow and fast fluctuations observable in the voltage traces in Fig 16D. Increasing τ IE while keeping τ EI = τ II = 3 msec has a substantial effect on the slow component but hardly affects the fast component. When plotted vs. τ/τ IE , the slow components of the PACs all collapse onto the same curve (Fig 16C, inset). This indicates that the EIE loop is essential in generating the slow, but not the fast, fluctuations. Fitting this slow component with the function AÁ[cosh(τ/τ dec )] −1 yields τ dec % 2.4τ IE . Furthermore, increasing τ II suppresses the fast fluctuations and amplifies the slow ones. These two effects saturate simultaneously when τ II % 10 ms (S11 Text). Thus, it can be inferred that fast fluctuations are mostly generated by II interactions. Their amplitude is suppressed as τ II is increased because they become more filtered. Concomitantly, the slow fluctuations become amplified. This is because fast fluctuations smooth the effective transfer function of the E neurons in the low firing rate regime. Thus, their suppression increases the gain of this transfer function. This explains the quantitative differences between the PACs in the spiking and the rate LIF network when II synapses are fast and why these differences are lessened as τ II increases (S11 Text).
In the simulations reported in Fig 16 there is no recurrent excitation in the E population (J EE 0 ¼ 0). Moreover, all the excitatory synapses to the I population are slow. Both assumptions were made to reduce the number of parameters in order to simplify the analysis. However, in cortical networks in general, fast (AMPA) and slow (NMDA) excitation coexist (in fact AMPA synapses are required to open the NMDA receptors). Moreover, recurrent excitation is thought to be in general substantial (see however [47]). Results depicted in S12 Text show that the EIE loop can induce slow rate fluctuations in our network when it combines slow and fast excitatory synapses and when substantial recurrent excitation is present in the E population.

Discussion
Networks of neurons operating in the so-called balanced regime exhibit spiking activity with strong temporal variability and spatial heterogeneity. Previous theoretical studies have investigated this regime assuming that excitatory and inhibitory synapses are sufficiently fast compared to the neuronal dynamics. The nature of the balanced state is now fairly well understood in this case. By contrast, here we focused on networks in which some of the synapses are slow. To study the dynamics in these networks, we reduced them to a rate dynamics that we investigated by combining Dynamical Mean-Field Theory and simulations. Our key result is that when synaptic interactions are sufficiently strong and slow, chaotic fluctuations on the time scales of the synaptic dynamics emerge naturally from the network collective behavior. Moreover, the nature of the transition to chaos and the behavior in the chaotic regime are determined only by the neuronal f − I curve and not by the details of the spike-generation mechanism.
We identified two mechanisms for the emergence of asynchronous chaos in EI neuronal networks. One mechanism relies on II interactions whereas in the other the EIE feedback loop plays the key role. These mechanisms hold in rate models (Eq (3)) as well as in LIF spiking networks. By computing the maximum Lyapunov exponent, we provided direct evidence that in rate models these states are indeed chaotic. For LIF spiking networks, we argued that when the synapses are sufficiently slow, the observed activity fluctuations are chaotic since their statistics are quantitatively similar to those observed in the corresponding rate model. This similarity persists for synaptic time constants as small as the membrane time constant. This is in agreement with [33][34][35] which relied on numerical integration of the LIF model to compute the Lyapunov spectra of networks of various sizes and increasing synaptic time constants. They found that the LIF dynamics are chaotic only if the synapses are sufficiently slow.
In these two mechanisms, the dynamics of the synaptic currents play the key role whereas dependence on the intrinsic properties of the neurons only occurs via their nonlinear instantaneous input-output transfer function. Since the synaptic currents are filtered versions of the neuronal spike trains, and that the temporal fluctuations of the activity occur on the time scales of the synaptic currents, it is natural to qualify the dynamical regime as rate chaos. Although the features of the bifurcation to chaos may depend on the shape of the transfer function, as we have shown, the qualitative features of the chaotic state are very general, provided that the synaptic currents are sufficiently slow. Rate chaos is therefore a generic property of networks of spiking neurons operating in the balanced regime. We show in S3 Text that rate chaos occurs also in networks of non-leaky integrate-and-fire spiking neurons. In that case, the statistics of the fluctuations are similar to those of the model in Eq (3) with a threshold-linear transfer function. We also found rate chaos in biophysically more realistic network models in which the dynamics of the neurons and of the synapses are conductance-based (results not shown). In these cases, the dynamics of the synaptic conductances give rise to the chaotic fluctuations.
Quantitative mappings from spiking to rate models have been derived for networks in stationary asynchronous non chaotic states [38] or responding to external fluctuating inputs [48]. Spiking dynamics also share qualitative similarities with rate models for networks operating in synchronous states [9-11, 38, 43]. To our knowledge, the current study is the first to report a quantitative correspondance between spiking and rate model operating in chaotic states.
The SCS model [19] has been widely used to explore the physiological [22,49] and computational significance of chaos in neuronal networks. Recent works have shown that because of the richness of its chaotic dynamics, the SCS model has remarkable learning capabilities [15][16][17][18]. Our work paves the way for an extension of these results to networks of spiking neurons with a connectivity satisfying Dale's law, which are biologically more realistic than the SCS model.
Another interesting implication of our work is in the field of random matrices. Given a dense NxN random matrix, A, with i.i.d elements with zero mean and finite standard deviation (SD), in the large N limit, the eigenvalue of A= ffiffiffiffi N p with the largest real part is real, and it is equal to SD [50,51] (more generally, the eigenvalues of A= ffiffiffiffi N p are uniformly distributed within a disk of radius SD centered at the origin [50,51]). Several results regarding the spectra (bulk and outliers) of dense random matrices with structures reflecting Dale's law have been derived recently [52][53][54]. Less is known when the matrices are sparse. A byproduct of our approach are two conjectures for the maximal eigenvalue of such sparse random matrices, namely Eqs (7) and (62) that we verified numerically.
Neuronal spiking statistics (e.g., firing rate, spike counts, inter-spike intervals) exhibit a very broad range of time scales during spontaneous or sensory evoked activity in-vivo (see e.g [55,56]). Fluctuations on time scales larger than several 100s of millisecond can be accounted for by neuromodulation which changes the global excitability of the cortical network or changes in behavioral state. Very fast fluctuations are naturally explained in the framework of the standard model of balance of excitation and inhibition [28][29][30]. By contrast, it is unclear how to explain modulations in the intermediate temporal range of a few 10s to several 100s of milliseconds. In fact, the standard framework of balanced networks predicts that fluctuations on this time scale are actively suppressed because the network state is very stable. Our work extends this framework and shows two mechanisms by which modulations in this range can occur. In the II mechanism, inhibitory synapses must be strong and slower than 10 − 20 ms. GABA A inhibition may be too fast for this [57] (see however [58]), but GABA B [59] are sufficiently slow. In contrast, the EIE mechanism is achieved when inhibition in fast. It requires slow recurrent excitation to inhibitory neurons, with a time constant of a few to several tens of ms, as is typically the case for NMDA receptors (see e.g [60][61][62]). Hence, the combination of GABA A and NMDA synapses can generate chaotic dynamics in the cortex and fluctuations in activity on a time scale of several tens to a few hundreds of ms.
Note added in production: Following a request from the editors after formal acceptance of our article, we note that a recent paper [63] claims that spiking networks with instantaneous delayed synapses exhibit an asynchronous state similar to the chaotic state of the SCS model. However, this claim is incorrect and has been shown to rely on flawed analysis [64].

Models
Two population leaky integrate-and-fire spiking network The two population network of leaky integrate-and-fire (LIF) neurons considered in this work consists of N E excitatory (E) and N I inhibitory neurons. The subthreshold dynamics of the membrane potential, V a i , of neuron i in population α (i = 1, . . ., N α ; α, β 2 {E, I}) obeys: where τ m is the membrane time constant (we take τ m = 10 msec for both populations), C ab ij and J αβ are respectively the connectivity matrix and the strength of the connections between the (presynaptic) population β and (postsynaptic) population α and I α the external feedforward input to population α. For simplicity we take N E = N I = N. However, all the results described in the paper are also valid when the number of neurons is different in the populations (provided both numbers are large)., The variables S ab j , which describe the synapses connecting neuron j in population β to population α, follow the dynamics: where τ αβ is the synaptic time constant and the sum is over all the spikes emitted at times t b j < t. Eqs (11), (12) are supplemented by a reset condition. If at time t sp , V a i ðt sp Þ ¼ 1, the neuron emits a spike and V a i ðt þ sp Þ ¼ 0. For simplicity we do not include the neuronal refractory period. We assume that the connectivity is random with all the C ab ij uncorrelated and such that C ab ij ¼ 1 with probability K/N and 0 otherwise. Hence each neuron is connected, on average, to K neurons from its population as well as to K neurons from the other population. When varying the connectivity K we scale the interaction strength and the feedforward inputs according to:

Network of inhibitory leaky integrate-and-fire neurons
The dynamics of the network of the one-population spiking LIF neurons considered in the first part of the paper are: supplemented with the reset condition at threshold. The elements of the connectivity matrix, C ij , are uncorrelated and such that C ij = 1 with probability K/N and 0 otherwise. All neurons are inhibitory, thus J < 0. The synaptic dynamics are: where τ syn is the synaptic time constant of the inhibition and the sum is over all the spikes emitted at times t j < t. The interaction strength and the feedforward inputs scale with K as:

Network of non-leaky integrate-and-fire neurons
We consider briefly this model in S3 Text. The network architecture as well as the synaptic dynamics are as above. The single neuron dynamics of non-leaky integrate-and-fire (NLIF) neurons are similar to those of LIF neurons except for the first terms on the right-hand side of Eqs (11), (13) which are now omitted.

Rate dynamics for spiking networks with slow synapses
If the synapses are much slower than the membrane time constant, the full dynamics of a spiking network can be approximated by the dynamics of the synapses driven by the instantaneous firing rates of the neurons, namely: where g(x) is the transfer function of the neuron (the f − I curve) [20]. In particular, for the LIF networks, with H(x) = 1 for x > 0 and H(x) = 0 otherwise. For the NLIF networks, the transfer function is threshold-linear: g(x) = xH(x).
Defining h ab i ≜J ab P j C ab ij S ab j , the dynamics of h ab i are given by We will denote by h b i the total input into neuron i in population β: For networks comprising only one population of inhibitory spiking neurons we will drop the superscript β = I and denote this input by h i . The dynamics then yield: where τ syn is the inhibitory synaptic time constant.

Dynamical Mean-Field Theory of the Single Inhibitory Population
A Dynamical Mean-Field Theory (DMFT) can be developed to investigate the rate model, Eq (17), for a general transfer function under the assumption, 1 ( K ( N. Here we provide a full analysis of a one-population network of inhibitory neurons whose dynamics are given in Eq (18). We take I ¼ I 0 ffiffiffi ffi K p as the external input and J ¼ J 0 = ffiffiffi ffi K p as the coupling strength. In this case, a functional integral derivation shows that these dynamics can be written as: where η i (t) is a Gaussian noise: with z i , i.i.d Gaussian quenched variables with zero mean and unit standard deviation (SD), ξ i (t) are Gaussian noises with hξ i (t)i t = 0, and hξ i (t)ξ j (t+τ)i t = C ξ (τ)δ i, j where h Á i t stands for averaging over time. Therefore, in general, the inputs to the neurons display temporal as well as quenched fluctuations.
The self-consistent equations that determine the mean, temporal correlations and quenched fluctuations yield: where h Á i and [Á] stand for averaging over noise and quenched disorder, respectively. Thus the quantities q and μ obey: and: where σ(τ) = [hh(t)h(t+τ)i] − μ 2 is the population-averaged autocovariance (PAC) of the input to the neurons and we define: σ 0 = σ(0) and Dx ¼ e À x 2 2 ffiffiffi ffi 2p p dx. In the limit K ! 1, μ must remain finite. This implies that the population averaged firing rate, [hg(h)i] = I 0 /J 0 does not depend on the specifics of the transfer function of the neurons and varies linearly with I 0 . This is a key outcome of the balance between the feedforward excitatory and the recurrent inhibitory inputs to the neurons.
To express C ξ (τ) in terms of σ, we note that the vector (h(t), h(t+τ)) T is a bivariate Gaussian, so in fact we need to calculate E[g(μ+x)g(μ+y)] where (x, y) T has zero mean and a covariance matrix with G(x) = R g(x)dx. Note that for positive σ this equation yields Therefore the quantity is conserved under the dynamics, Eq (29). Hence: To simplify notations, we drop the parameter σ 0 and denote the potential by V(σ). The first, second and third order derivatives of the potential with respect to σ are denoted V 0 (σ), V 00 (σ) and V 000 (σ).
A bifurcation between these behaviors occurs at some critical value, J c , such that for J 0 < J c the self-consistent solutions of Eq (29) are either oscillatory or constant as a function of τ, whereas for J 0 > J c they are either oscillatory or decay monotonically. A stability analysis of these different solutions is beyond the scope of this paper; instead, we rely on numerical simulations of the full dynamics. They indicate that the network dynamics always reach a fixed point for sufficiently small J 0 . For sufficiently large J 0 the fixed point is unstable and the network settles in a state in which σ(τ) decays monotonically with τ. Simulations also show that the maximum Lyapunov exponent in these cases is positive (see below); i.e. the network is in a chaotic state. For values of J 0 in between these two regimes, the network displays oscillatory patterns of activity. However, for increasing network sizes, N, the range of J 0 in which oscillations are observed vanishes (not shown). Therefore for large N the bifurcation between a fixed point and chaos occurs abruptly at some critical value J c . A similar phenomenology occurs for other non-linear positive monotonically increasing transfer functions.
In summary, for a fixed feedforward input, I 0 , there are two regimes in the large N limit: 1. for J 0 < J c : the stable state is a fixed point. The distribution of the inputs to the neurons is a Gaussian whose mean, μ, and variance, σ are determined by the self-consistent mean-field equations: For a transfer function, g(x), which is zero when x is smaller than some threshold T (functions without threshold correspond to T = −1), the distribution of the neuronal firing rates, r i , in this state is given by: 2. for J 0 > J c : the stable state is chaotic. The distribution of time average inputs is Gaussian with mean μ and variance s 1 ¼ J 0 2 q and the autocovariance of the inputs is determined by Eq (29) which depends on σ 0 . The quantities μ, σ 0 and σ 1 are determined by the self-consistent equations: and s 0 2 À s 1 together with Eq (25).
(1+iτ αβ ω). Transforming back to the time domain yields: Since Δ αβ = σ αβ +(μ αβ ) 2 we get: Thus we get a set of self-consistent equations for the four PACs σ αβ . The relevant soutions have to satisfy the four boundary conditions: In general, these dynamical equations cannot be written like those of a particle in some potential. This makes the study of their solutions substantially more difficult than in the one population case.

Separation of time scales
A potential function can be written for the DMFT if the time scale of one type of synapses is substantially larger than the others, which makes it possible to consider the latter as instantaneous. We carry out this analysis below assuming τ IE ) τ EI , τ EE , τ II .
Setting all the synapses except those from E neurons to I neurons to be instantaneous implies that except for σ IE one has: whereC b is defined in Eq (44). Since τ IE is now the only time scale we can take τ IE = 1. Also, σ EE , σ EI , σ II and the potential V are now functions of a single variable, σ IE . Therefore, the differential equation for σ IE can be written as The instability of the fixed point occurs when, V 0 (σ IE ) and V 00 (σ IE ), the first and the second derivatives of V with respect to σ IE , vanishes. Using Eq (49) one has: Since σ α = σ αE +σ αI : where theC ab (α = E, I, β = E, I) are N×N sparse matrices with elements (C αβ is the matrix of connectivity between populations β (presynaptic) and α). We are interested in instability onsets at which a real eigenvalue crosses 0. Using Eq (56), it is straightforward to show that such an instability happens if the synaptic strength are such that: If J EE 0 ¼ 0, one can rewrite Eq (58) as: with: Let us assume that J II 0 is fixed and such that for small enough J IE 0 J EI 0 the fixed point is stable. When increasing, J IE 0 J EI 0 the fixed point loses stability when the value of J IE 0 J EI 0 is the smallest for which Eq (59) is satisfied, that is for which the largest real eigenvalue, λ max of the matrix M crosses 1. If this instability also corresponds to chaos onset, Eq (54), this would imply that the condition λ max = 1 is equivalent to: Interestingly, this condition means that the variance of the elements of the matrix ffiffiffiffi N p M is equal to one leading us to conjecture that more generally the eigenvalue of the latter which has the largest real part and is given by:

Numerical simulations Integration of network dynamics and mean-field equation solutions
The integration of differential equations, Eq (15) and Eq (18) (Eq (3) in main text), was performed with a C code using the Euler method with fixed Δt = τ syn /20 (the validity of the results was verified using smaller values of Δt).
Simulations of the LIF spiking networks were done using a second-order Runge-Kutta integration scheme supplemented by interpolation of spike times as detailed in [65]. In all the spiking network simulations the time step was Δt = 0.1 ms.
Self-consistent mean-field equations were solved with MATLAB function fsolve, which implements a 'trust-region-dogleg' algorithm or the Levenberg-Marquardt algorithm for nonsquare systems. Numerical calculations of integrals was done with MATLAB function trapz.

Population-averaged autocovariance
The population average autocovariance (PAC) functions of neuronal quantities f i (t) (i = 1. . .N) were computed as where N t is the number of time samples for the calculation of the PAC. In all figures f i (t) = h i (t) except in Fig 16 where f i ðtÞ ¼ I a þ h aE i ðtÞ À h aI i ðtÞ. All PACs of spiking networks were calculated over 163.84 sec, and averaged over 10 realizations of the connectivity. For models Eq (15) and Eq (18), PACs were calculated over 2048τ syn after discarding 200τ syn of transient dynamics and averaged over 8 realizations.

Largest Lyapunov exponents
To calculate the maximal Lyapunov exponent, Λ, of the inhibitory network, Eq (3), we simulated the system for a sufficiently long duration (200τ syn ) so that it settled on the attractor of the dynamics. Denoting byh Ã the network state at that time, we then ran two copies of the dynamics, one with initial conditionsh 1 ðt ¼ 0Þ ¼h Ã and the other with slightly perturbed initial conditions,h 2 ðt ¼ 0Þ ¼h Ã þ = ffiffiffiffi N p (jjh 1 ð0Þ Àh 2 ðð0Þ jj¼ , where jjÁjj is the l 2 norm).
Monitoring the difference,dðtÞ ¼h 1 ðtÞ Àh 2 ðtÞ we computed T ð1Þ reset ¼ minðargðjjdðtÞ jj¼ D max Þ; T max Þ and D ð1Þ reset ¼jjdðT ð1Þ reset Þ jj. We then reinitialized the dynamics of the second network copy toh 2 ðT ð1Þ reset Þ þd Fraction of networks with a stable fixed point in rate dynamics Fig 10D in the main text plots the lines in the J 0 − I 0 phase diagrams of the threshold-power law rate model, for which 5%,50%,95% of randomly chosen networks have dynamics which converge to a fixed point. To compute these lines we simulated, for each value of γ and J 0 , 100 realizations of the network. For each realization, we computed the population average of the temporal variance the synaptic inputs, ρ: where N tot is the total number of time steps of the simulations after discarding a transient with a duration of 256τ syn . The fixed point was considered to be unstable if ρ > 10 −9 . The fraction of unstable networks, F u , was fitted with a logistic function: S9 Text. The two mechanisms underlying asynchronous chaos in two-population LIF networks: Results for inhibitory neurons. (PDF) S10 Text. Two-population LIF rate and spiking models: In the II mechanism the PAC depends very mildly on τ IE . (PDF) S11 Text. Two-population LIF rate and spiking models: In the EIE mechanism the slow component of the PAC depends very mildly on τ II . (PDF) S12 Text. Two-population integrate-and-fire network with recurrent EE excitation, AMPA and NMDA synapses and fast inhibition.