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.
Single cell recordings [1] and electro-encephalography [2,3] suggest the existence of chaotic dynamics in the brain. Consistent with chaotic dynamics, in-vivo experiments have demonstrated that cortical circuits are sensitive to weak perturbations [4,5]. Remarkably, the misplacement of even a single spike in a cortical network has a marked effect on the timing of subsequent spikes in the network [6]. Chaotic states in extended dynamical systems can be classified as synchronous or asynchronous, depending on the spatial patterns of the dynamics. In synchronous chaos the temporal fluctuations exhibit spatial correlations. If the temporal fluctuations are spatially incoherent, the chaotic state is classified as asynchronous. EEG measures the activity of a large population of neurons. Therefore, it is probable that chaoticity observed in EEGs reflects synchronous chaos in brain regions of rather large size. Models of local cortical circuits exhibiting synchronous chaos have been studied in [7–12]. A computational advantage of synchronous chaos in the brain is that it enables neuronal populations to respond quickly to changes in their external inputs [7] and facilitates the access of the network to states (e. g. limit cycles or fixed points) that encode different stimuli [3]. A large body of experimental data, however, has reported that cortical neurons exhibit very weak correlations [13,14] and thus are more compatible with asynchronous than with synchronous chaos. Moreover, recent studies have demonstrated that the richness, the complexity and the high dimension of the dynamics in systems operating in asynchronous chaos endows them with remarkable computational capabilities [15–18]. The present paper focuses on the mechanisms underlying the emergence of asynchronous chaos in local neuronal circuits. Asynchronous chaos was studied in a seminal work by Sompolinsky, Crisanti and Sommers (SCS) [19], who investigated a large network of N neuronal-like units fully connected with random weights drawn from a zero mean Gaussian distribution (called hereafter as the SCS model). The dynamics of the network are those of a “rate” model [20], in which the activity of a unit, S (t), is characterized by a continuous variable which is a non-linear function, S = ϕ (h), of the total input to the unit. In the SCS model the activity variables take values between [-1,1] and the function ϕ (h) is sigmoidal and odd. Using Dynamical Mean-Field Theory (DMFT) SCS showed that if the standard deviation of the weight distribution is sufficiently large, the dynamics bifurcate from fixed point to asynchronous chaos. The SCS model in its original form or in its discrete time version has been used in numerous studies in theoretical and computational neuroscience [15–17,21–25]. However, the connectivity of the SCS model violates Dale’s Law, whereby in biological networks a given neuron is either excitatory or inhibitory [26]. Also, the equation of the SCS model dynamics are invariant under the transformation h → −h, a symmetry not fulfilled in more realistic neuronal network models. More importantly, as this is the case frequently for rate models, the physiological meanings of the dynamical “neuronal” variables and of the parameters are not clear in the SCS network. Should these variables and the time constant of their dynamics—which sets the time scale of the chaotic fluctuations—be interpreted as characterizing neurons, or synapses? In this paper we address the following general and fundamental issues: To what extent are asynchronous chaotic states generic in networks of spiking neurons? How does this depend on 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].
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 Jij = JCij (i, j = 1,2…N), with J ≤ 0, and Cij = 1 with probability K/N and 0 otherwise. The outgoing synapses of neuron j obey τ s y n d S j (t) d t = - S j (t) + J ∑ t j s < t δ (t - t j s) (1) where Sj (t) is the synaptic current at time t and τsyn the synaptic time constant. When neuron j fires a spike (time t j s), Sj increments by J. Thus, the total input to neuron i, hi (t) = I+∑j Jij Sj (t), satisfies: τ s y n d h i (t) d t = - h i (t) + I + ∑ j ∑ t j s < t J i j δ (t - t j s) (2) We assume K ≫ 1, hence the number of recurrent inputs per neuron is K ± O (K). Scaling J and I as: J = − J 0 / K, I = K I 0, the time-averaged synaptic inputs are O (K) 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–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: τ s y n d h i (t) d t = - h i (t) + I + ∑ j J i j r j (t) (3) r i (t) = g (h i (t) ) (4) where ri (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. 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. We now consider EI spiking networks with recurrent feedback interactions between the two populations. The synaptic strengths and time constants are J 0 α β / K and ταβ (α, β ∈ {E, I}). Assuming slow synapses, the dynamics can be reduced to four sets of equations for the four types of synaptic inputs, h α β i (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).
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–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–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 / N with the largest real part is real, and it is equal to SD [50,51] (more generally, the eigenvalues of A / N 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–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–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. GABAA inhibition may be too fast for this [57] (see however [58]), but GABAB[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–62]). Hence, the combination of GABAA 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].
The two population network of leaky integrate-and-fire (LIF) neurons considered in this work consists of NE excitatory (E) and NI inhibitory neurons. The subthreshold dynamics of the membrane potential, V i α, of neuron i in population α (i = 1, …, Nα; α, β ∈ {E, I}) obeys: τ m d V i α (t) d t = - V i α (t) + I α + J α E ∑ j C i j α E S j α E (t) - ∑ j J α I C i j α I S j α I (t) (11) where τm is the membrane time constant (we take τm = 10 msec for both populations), C i j α β 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 NE = NI = 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 j α β, which describe the synapses connecting neuron j in population β to population α, follow the dynamics: τ α β d S j α β d t = - S j α β + ∑ t j β δ (t - t j β) (12) where ταβ is the synaptic time constant and the sum is over all the spikes emitted at times t j β < t. Eqs (11), (12) are supplemented by a reset condition. If at time tsp, V i α (t s p) = 1, the neuron emits a spike and V i α (t s p +) = 0. For simplicity we do not include the neuronal refractory period. We assume that the connectivity is random with all the C i j α β uncorrelated and such that C i j α β = 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: J α β = J 0 α β / K and I α = I 0 α K[29]. The dynamics of the network of the one-population spiking LIF neurons considered in the first part of the paper are: τ m d V i (t) d t = - V i (t) + I + J ∑ j C i j S j (t) (13) supplemented with the reset condition at threshold. The elements of the connectivity matrix, Cij, are uncorrelated and such that Cij = 1 with probability K/N and 0 otherwise. All neurons are inhibitory, thus J < 0. The synaptic dynamics are: τ s y n d S j d t = - S j + ∑ t j δ (t - t j) (14) where τsyn is the synaptic time constant of the inhibition and the sum is over all the spikes emitted at times tj < t. The interaction strength and the feedforward inputs scale with K as: J = − J 0 / K and I = I 0 K with J0 > 0. 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. 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: τ α β d S i α β d t = - S i α β + g (J β E ∑ j C i j β E S j β E - J β I ∑ j C i j β I S j β I + I β) (15) where g (x) is the transfer function of the neuron (the f − I curve) [20]. In particular, for the LIF networks, g (x) = - 1 τ m log (1 - 1 / x) H (x - 1) (16) 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 i α β ≜ J α β ∑ j C i j α β S j α β, the dynamics of h i α β are given by τ α β d h i α β d t = - h i α β + ∑ j J α β C i j α β g (h j β E (t) - h j β I (t) + I β) (17) We will denote by h i β the total input into neuron i in population β: h i β = h i β E − h i β I + I β. For networks comprising only one population of inhibitory spiking neurons we will drop the superscript β = I and denote this input by hi. The dynamics then yield: τ s y n d h i d t = - h i + I - J ∑ j = 1 N C i j g (h j) (18) where τsyn is the inhibitory synaptic time constant.
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 K as the external input and J = J 0 / K as the coupling strength. In this case, a functional integral derivation shows that these dynamics can be written as: τ s y n d h i (t) d t = - h i (t) + η i (t), i = 1, . . . , N (19) where ηi (t) is a Gaussian noise: η i (t) = μ + J 0 q z i + ξ i (t) (20) with zi, i. i. d Gaussian quenched variables with zero mean and unit standard deviation (SD), ξi (t) are Gaussian noises with ⟨ξi (t) ⟩t = 0, and ⟨ξi (t) ξj (t+τ) ⟩t = Cξ (τ) δi, j where ⟨ ⋅ ⟩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: μ = K (I 0 - J 0 [ ⟨ g (h i (t) ) ⟩ ]) (21) q = [ ⟨ g (h) ⟩ 2 ] (22) C ξ (τ) = J 0 2 ([ ⟨ g (h (t) ) g (h (t + τ) ⟩ ] - q) (23) where ⟨ ⋅ ⟩ and [⋅] stand for averaging over noise and quenched disorder, respectively. Thus the quantities q and μ obey: q = ∫ - ∞ ∞ [ ∫ - ∞ ∞ g (μ + J 0 q z + σ 0 - J 0 2 q ξ) D ξ ] 2 D z (24) and: 1 J 0 (I 0 - μ K) = ∫ - ∞ ∞ g (μ + σ 0 z) D z (25) where σ (τ) = [⟨h (t) h (t+τ) ⟩] − μ2 is the population-averaged autocovariance (PAC) of the input to the neurons and we define: σ0 = σ (0) and D x = e − x 2 2 2 π d x. In the limit K → ∞, μ must remain finite. This implies that the population averaged firing rate, [⟨g (h) ⟩] = I0/J0 does not depend on the specifics of the transfer function of the neurons and varies linearly with I0. 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 Σ x y = (σ 0 σ σ σ 0) and E[⋅] stands for averaging over temporal noise and quenched disorder. Defining [ x y ] = [ σ 0 - | σ | 0 | σ | 0 σ 0 - | σ | sign (σ) · | σ | ] · [ ξ θ z ] where ξ, θ and z are independent Gaussian variables with zero mean and unit variance yields E [ g (μ + x) g (μ + y) ] = = E [ E [ g (μ + σ 0 - | σ | ξ + | σ | z) | z ] E [ g (μ + σ 0 - | σ | θ + sign (σ) · | σ | z) | z ] ] = = ∫ - ∞ ∞ [ ∫ - ∞ ∞ g (μ + σ 0 - | σ | ξ + | σ | z) D ξ · ∫ - ∞ ∞ g (μ + σ 0 - | σ | θ + sign (σ) · | σ | z) D θ ] D z (26) A straightforward derivation shows that σ (τ) obeys: τ s y n 2 d 2 σ d τ 2 = = σ - J 0 2 ∫ - ∞ ∞ ∫ - ∞ ∞ ∫ - ∞ ∞ g (μ + σ 0 - | σ | ξ + | σ | z) g (μ + σ 0 - | σ | θ + sign (σ) · | σ | z) D ξ D θ D z (27) with initial conditions: σ (0) = σ 0; d σ d τ (0) = 0 (28) where the last condition results from σ (τ) = σ (−τ). Eq (27) can be rewritten as: τ s y n 2 d 2 σ d τ 2 = - ∂ V (σ; σ 0) ∂ σ (29) where the “potential” V (σ; σ0) which depends parametrically on σ0 is: V (σ; σ 0) = = - σ 2 2 + J 0 2 ∫ - ∞ ∞ ∫ - ∞ ∞ ∫ - ∞ ∞ G (μ + σ 0 - | σ | ξ + | σ | z) G (μ + σ 0 - | σ | θ + sign (σ) · | σ | z) D ξ D θ D z (30) with G (x) = ∫g (x) dx. Note that for positive σ this equation yields V (σ; σ 0) = - σ 2 2 + J 0 2 ∫ - ∞ ∞ [ ∫ - ∞ ∞ G (μ + σ 0 - σ ξ + σ z) D ξ ] 2 D z (31) Therefore the quantity E = 1 2 (τ s y n d σ d τ) 2 + V (σ; σ 0) (32) is conserved under the dynamics, Eq (29). Hence: 1 2 (τ s y n d σ d τ) 2 + V (σ; σ 0) = V (σ 0; σ 0) (33) 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′ (σ), V′′ (σ) and V′′′ (σ). For illustrative purpose, we consider a sigmoid transfer function, g (x) = ϕ (x) ≜ 1 2 [ 1 + erf (x 2) ]. In this case we have G (x) = Φ (x) ≜ x 2 [ 1 + erf (x 2) ] + e - x 2 2 2 π Using the identities: ∫ - ∞ ∞ ϕ (a + b z) D z = ϕ (a 1 + b 2) and ∫ - ∞ ∞ ϕ (a + b z) z D z = b 1 + b 2 e - a 2 2 (1 + b 2) 2 π the potential V (σ) can be written as: V (σ) = - σ 2 2 + J 0 2 ∫ - ∞ ∞ (1 + σ 0 - | σ |) Φ (μ + | σ | z 1 + σ 0 - | σ |) Φ (μ + sign (σ) · | σ | z 1 + σ 0 - | σ |) D z Fig 17A1–3 plots V for σ ∈ (−σ0, σ0) for J0 = 4, fixed I0 = 1 and different values of σ0. When V′ (σ0) > 0 (Fig 17A1), the solution to Eq (29), σ (τ), decreases monotonically from σ0 to −σ0 that it reaches in finite time with a strictly negative velocity; this solution does not correspond to an autocovariance function. For σ0 such that V′ (σ0) = 0 (Fig 17A2) the solution is σ (τ) = σ0. It corresponds to a fixed point of the dynamics, Eq (18) in which all the inputs to the neurons are constant in time, h i (t) = h i 0, and h i 0 has a Gaussian distribution. Finally, for σ0 such that V′ (σ0) < 0 (Fig 17A3), there is no solution to Eq (33) with σ (0) = σ0. Fig 17B1–3 plots V for J0 = 15. For small σ0, the solution Eq (33) does not correspond to an autocovariance function. As σ0 increases, V (σ) becomes non-monotonic in the vicinity of σ = σ0 with local maxima and minima at σ = σmax and σ = σmin, respectively (Fig 17B2). However, here also the solution for σ (τ) does not correspond to an autocovariance because σ0 is the global maximum in the range σ ∈ [−σ0, σ0]. For σ 0 = σ 0 *, such that V (σ m a x; σ 0 *) = V (σ 0 *; σ 0 *) (Fig 17B3) an acceptable solution appears, in which σ decays monotonically from σ 0 * and converges to σmax as τ → ∞, i. e. σmax = σ∞. This solution corresponds to a chaotic state of the network. If σ0 is further increased beyond σ 0 *, V (σmax, σ0) > V (σ0) (Fig 17B4), and the solution exhibits oscillations around σmin. For σ0 ≈ 11. 77, V′ (σ0) = 0, and the solution corresponds to a fixed point (Fig 17B5). Finally, for σ0 larger, V′ (σ0) is negative (Fig 17B6) and there is no solution to Eq (18) with σ (0) = σ0. A bifurcation between these behaviors occurs at some critical value, Jc, such that for J0 < Jc the self-consistent solutions of Eq (29) are either oscillatory or constant as a function of τ, whereas for J0 > Jc 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 J0. For sufficiently large J0 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 J0 in between these two regimes, the network displays oscillatory patterns of activity. However, for increasing network sizes, N, the range of J0 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 Jc. A similar phenomenology occurs for other non-linear positive monotonically increasing transfer functions. In summary, for a fixed feedforward input, I0, there are two regimes in the large N limit: for J0 < Jc: 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: μ = K (I 0 - J 0 ∫ - ∞ ∞ g (μ + σ z) D z) (34) σ = J 0 2 ∫ - ∞ ∞ [ g (μ + σ z) ] 2 D z (35) For a transfer function, g (x), which is zero when x is smaller than some threshold T (functions without threshold correspond to T = −∞), the distribution of the neuronal firing rates, ri, in this state is given by: p m (x) = d d x [ Pr (r i ≤ x) ] = = 1 2 π σ e - μ 2 2 σ · δ (x - T) + 1 2 π σ e - (g - 1 (x) - μ) 2 2 σ · 1 g ′ (g - 1 (x) ) · H (x - T) (36) for J0 > Jc: the stable state is chaotic. The distribution of time average inputs is Gaussian with mean μ and variance σ ∞ = J 0 2 q and the autocovariance of the inputs is determined by Eq (29) which depends on σ0. The quantities μ, σ0 and σ∞ are determined by the self-consistent equations: σ ∞ = J 0 2 ∫ - ∞ ∞ [ ∫ - ∞ ∞ g (μ + σ 0 - σ ∞ ξ + σ ∞ z) D ξ ] 2 D z (37) and σ 0 2 - σ ∞ 2 2 = = J 0 2 ∫ - ∞ ∞ G (μ + σ 0 z) 2 D z - - J 0 2 ∫ - ∞ ∞ [ ∫ - ∞ ∞ G (μ + σ 0 - σ ∞ ξ + σ ∞ z) D ξ ] 2 D z (38) together with Eq (25).
A DMFT approach can also be developed to investigate the dynamics of the two population network model, Eq (17). To that end, the last term in Eq (17) is written as a Gaussian random process with mean μαβ and autocorrelation function Cαβ (τ) and derives the self-consistent equations that these quantities satisfy. The quantity μαβ is therefore μ α β = J 0 α β K E [ g (h β) ] where: h i β = h i β E - h i β I + I β (39) is the net input to neuron i in population β. The synaptic inputs h i α β is also a Gaussian random process. We denote its mean over time and over all the neurons in population α by μαβ = E[hαβ (t) ] and its PAC by σαβ (τ) = E[hαβ (t) hαβ (t+τ) ]− (μαβ) 2. Taking I β = I 0 β ⋅ K we can write the mean of h j β as μ β = μ β E - μ β I + I β = = K (J 0 α E E [ g (h E) ] - J 0 α I E [ g (h I) ] + I 0 β) (40) The PAC of h j β then reads: σ β (τ) ≜ E [ h β (t) h β (t + τ) ] - (μ β) 2 = = σ β E (τ) + σ β I (τ) We can now write the balance condition in the large K limit: I 0 α + J 0 α E r E - J 0 α I r I = μ α K (41) where r β = E [ g (h α) ] = ∫ - ∞ ∞ g (μ α + σ 0 β z) e - z 2 2 2 π d z (42) is the neuronal firing rate averaged over cells in population α. Here, σ 0 β = σ β (0). We can also express Cαβ (τ) in terms of σα (τ) as: C α β (τ) = E [ ∑ j J i j α β g (h j β (t) ) ∑ j J i j α β g (h j β (t + τ) ) ] = (J 0 α β) 2 C ˜ β (σ β (τ) ) + (μ α β) 2 (43) where: C ˜ β (σ β) = = ∫ - ∞ ∞ ∫ - ∞ ∞ ∫ - ∞ ∞ g (μ β + σ 0 β - σ β ξ + sign (σ β) | σ β | z) g (μ β + σ 0 β - σ β θ + | σ β | z) D θ D ξ D z (44) Let us denote by Δαβ (τ) the autocorrelation of hαβ (t). We can express the relation between Cαβ (τ) and Δαβ (τ) by their Fourier transforms as Δαβ (ω) = H (ω) H* (ω) Cαβ (ω), where H (ω) = 1/ (1+iταβ ω). Transforming back to the time domain yields: (τ α β) 2 d 2 Δ α β d τ 2 = Δ α β - C α β (45) Since Δαβ = σαβ+ (μαβ) 2 we get: (τ α β) 2 d 2 σ α β d τ 2 = σ α β - (J 0 α β) 2 C ˜ β (46) Thus we get a set of self-consistent equations for the four PACs σαβ. The relevant soutions have to satisfy the four boundary conditions: lim τ → ∞ d σ α β (τ) d τ = 0 (47) 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. 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: σ α β = (J 0 α β) 2 C ˜ β (48) where C ˜ β 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 d 2 σ I E d τ 2 = - d V d σ I E where d V d σ I E = - σ I E + (J 0 I E) 2 C ˜ E (σ I E) (49) The instability of the fixed point occurs when, V′ (σIE) and V′′ (σIE), the first and the second derivatives of V with respect to σIE, vanishes. Using Eq (49) one has: V ′ ′ (σ I E) = - 1 + (J 0 I E) 2 d C ˜ E d σ E · d σ E d σ I E (50) Since σα = σαE+σαI: d σ E d σ I E = (J 0 E E) 2 d C ˜ E d σ E · d σ E d σ I E + (J 0 E I) 2 d C ˜ I d σ I · d σ I d σ I E (51) and d σ I d σ I E = 1 + (J 0 I I) 2 d C ˜ I d σ I · d σ I d σ I E (52) where d C ˜ β d σ β = ∫ - ∞ ∞ [ ∫ - ∞ ∞ g ′ (μ β + σ 0 β - σ β ξ + σ β z) e - ξ 2 2 2 π d ξ ] 2 e - z 2 2 2 π d z From Eqs (51) – (52) one gets: d σ E d σ I E = (J 0 E I) 2 d C ˜ I d σ I (1 - (J 0 E E) 2 d C ˜ E d σ E) (1 - (J 0 I I) 2 d C ˜ I d σ I) and: V ′ ′ (σ I E) = - 1 + (J 0 I E) 2 d C ˜ E d σ E (J 0 E I) 2 d C ˜ I d l σ I (1 - (J 0 E E) 2 d C ˜ E d σ E) (1 - (J 0 I I) 2 d C ˜ I d σ I) (53) Thus at chaos onset, together with Eq (41), J 0 α β, σα and μα obey: σ α = (J 0 α E) 2 C ^ (μ E, σ E) + (J 0 α I) 2 C ^ (μ I, σ I) (54) 1 = (J 0 E E) 2 C ^ ′ (μ E, σ E) + (J 0 I I) 2 C ^ ′ (μ I, σ I) + + [ (J 0 E I J 0 I E) 2 - (J 0 E E J 0 I I) 2 ] C ^ ′ (μ E, σ E) C ^ ′ (μ I, σ I) (55) where: C ^ (μ, σ) = ∫ - ∞ ∞ [ g (μ + | σ | z) ] 2 e - z 2 2 2 π d z C ^ ′ (μ, σ) = ∫ - ∞ ∞ [ g ′ (μ + | σ | z) ] 2 e - z 2 2 2 π d z For instance for the threshold-linear transfer function we have C ^ (μ, σ) = F 2 (μ, σ) C ^ ′ (μ, σ) = ϕ (μ σ) and r α = F 1 (μ α, σ α) where Fi (a, b) are defined in Eq (28). It should be noted that if the transition to chaos occurs for the same parameters for which the fixed point loses stability and that this is controlled by a real eigenvalue crossing zero, the location of the transition will not depend on the synaptic time constant. If this is the case, Eq (54) will characterize the location of the transition to chaos in the parameter space of the network in general and not only under the assumption of the separation of time scales under which we have established this condition. Let us denote the fixed point solution of the dynamics, Eq (17), by: h α β (t) = h ‾ α β. Writing h α β (t) = h ‾ α β + δ h α β with δ h α β ≪ h ‾ α β, linearizing the dynamics and looking for solution of the form δ h ∝ eλt) one gets: λ τ E E δ h E E = - δ h E E + J 0 E E C ˜ E E (δ h E E - δ h E I) λ τ I E δ h I E = - δ h I E + J 0 I E C ˜ I E (δ h E E - δ h E I) λ τ E I δ h E I = - δ h E I - J 0 E I C ˜ E I (δ h I E - δ h I I) λ τ I I δ h I I = - δ h I I - J 0 I I C ˜ I I (δ h I E - δ h I I) (56) where the C ˜ α β (α = E, I, β = E, I) are N×N sparse matrices with elements C ˜ i j E E = g ′ (h ¯ j E E - h ¯ j E I) C i j E E C ˜ i j I E = g ′ (h ¯ j E E - h ¯ j E I) C i j I E C ˜ i j E I = g ′ (h ¯ j I E - h ¯ j I I) C i j E I C ˜ i j I I = g ′ (h ¯ j I E - h ¯ j I I) C i j I I (57) (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: det [ I - J 0 I E J 0 E I C ˜ I E (I - J 0 E E C ˜ E E) - 1 C ˜ E I (I - J 0 I I C ˜ I I) - 1 ] = 0 (58) If J 0 E E = 0, one can rewrite Eq (58) as: det [ I - M ] = 0 (59) with: M = J 0 I I C ˜ I I + J 0 I E J 0 E I C ˜ I E C ˜ E I (60) Let us assume that J 0 I I is fixed and such that for small enough J 0 I E J 0 E I the fixed point is stable. When increasing, J 0 I E J 0 E I the fixed point loses stability when the value of J 0 I E J 0 E I 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: 1 = (J 0 I I) 2 C ^ ′ (μ I, σ I) + (J 0 E I J 0 I E) 2 C ^ ′ (μ E, σ E) C ^ ′ (μ I, σ I) (61) Interestingly, this condition means that the variance of the elements of the matrix N 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: λ m a x = (J 0 I I) 2 C ^ ′ (μ I, σ I) + (J 0 E I J 0 I E) 2 C ^ ′ (μ E, σ E) C ^ ′ (μ I, σ I) (62) 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 non-square systems. Numerical calculations of integrals was done with MATLAB function trapz. The population average autocovariance (PAC) functions of neuronal quantities fi (t) (i = 1…N) were computed as σ (τ) = σ (k Δ τ) = = 1 N ∑ i = 1 N 1 N t - | k | ∑ n = 0 N t - 1 f i (n Δ τ) f ( (n + k) Δ τ) - [ 1 N 1 N t ∑ i = 1 N ∑ n = 0 N t - 1 f i (n Δ τ) ] 2 where Nt is the number of time samples for the calculation of the PAC. In all figures fi (t) = hi (t) except in Fig 16 where f i (t) = I α + h i α E (t) − h i α 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. 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 by h ⃗ * the network state at that time, we then ran two copies of the dynamics, one with initial conditions h ⃗ 1 (t = 0) = h ⃗ * and the other with slightly perturbed initial conditions, h ⃗ 2 (t = 0) = h ⃗ * + ϵ / N (∣ ∣ h ⃗ 1 (0) − h ⃗ 2 ( (0) ∣ ∣ = ϵ, where ∣∣⋅∣∣ is the l2 norm). Monitoring the difference, d ⃗ (t) = h ⃗ 1 (t) − h ⃗ 2 (t) we computed T r e s e t (1) = min (arg (∣ ∣ d ⃗ (t) ∣ ∣ = D m a x), T m a x) and D r e s e t (1) = ∣ ∣ d ⃗ (T r e s e t (1) ) ∣ ∣. We then reinitialized the dynamics of the second network copy to h ⃗ 2 (T r e s e t (1) ) + d ⃗ (T r e s e t (1) ) ∣ ∣ d ⃗ (T r e s e t (1) ) ∣ ∣ ⋅ ϵ. We iterated the process n times and estimate the Lyapunov exponent according to: Λ = ∑ i = 1 n ln (D r e s e t (i) ϵ) ∑ i = 1 n T r e s e t (i) A similar method was used for two population networks, Eq (15), the only difference being that the vector h ⃗ now had dimension 4N. Throughout the article we take n = 100, Tmax = 5τsyn, Dmax = 10−3 and ϵ = 10−6. The Lyapunov exponent values reported in this article are averages over 5 realizations of the networks. Fig 10D in the main text plots the lines in the J0 − I0 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 J0,100 realizations of the network. For each realization, we computed the population average of the temporal variance the synaptic inputs, ρ: ρ = 1 N ∑ i = 1 N [ 1 N t o t ∑ k = 0 N t o t - 1 h i (k Δ t) 2 - (1 N t o t ∑ k = 0 N t o t - 1 h i (k Δ t) ) 2 ] where Ntot 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, Fu, was fitted with a logistic function: Fu (J0) = 100[1+exp (− (J0 − Jm) /ΔJ) ]−1. The thick red line and red dots plot the values of Jm vs. γ, and the dashed lines are the values of J0 for which Fu = 95 and Fu = 5.
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 pave the way in the study of the physiological mechanisms and computational significance of chaotic states in neuronal networks.
Asynchronous Rate Chaos in Spiking Neuronal Circuits