Transition to Chaos in Random Neuronal Networks

Firing patterns in the central nervous system often exhibit strong temporal irregularity and considerable heterogeneity in time-averaged response properties. Previous studies suggested that these properties are the outcome of the intrinsic chaotic dynamics of the neural circuits. Indeed, simplified rate-based neuronal networks with synaptic connections drawn from Gaussian distribution and sigmoidal nonlinearity are known to exhibit chaotic dynamics when the synaptic gain (i.e., connection variance) is sufficiently large. In the limit of an infinitely large network, there is a sharp transition from a fixed point to chaos, as the synaptic gain reaches a critical value. Near the onset, chaotic fluctuations are slow, analogous to the ubiquitous, slow irregular fluctuations observed in the firing rates of many cortical circuits. However, the existence of a transition from a fixed point to chaos in neuronal circuit models with more realistic architectures and firing dynamics has not been established. In this work, we investigate rate-based dynamics of neuronal circuits composed of several subpopulations with randomly diluted connections. Nonzero connections are either positive for excitatory neurons or negative for inhibitory ones, while single neuron output is strictly positive with output rates rising as a power law above threshold, in line with known constraints in many biological systems. Using dynamic mean field theory, we find the phase diagram depicting the regimes of stable fixed-point, unstable-dynamic, and chaotic-rate fluctuations. We focus on the latter and characterize the properties of systems near this transition. We show that dilute excitatory-inhibitory architectures exhibit the same onset to chaos as the single population with Gaussian connectivity. In these architectures, the large mean excitatory and inhibitory inputs dynamically balance each other, amplifying the effect of the residual fluctuations. Importantly, the existence of a transition to chaos and its critical properties depend on the shape of the single-neuron nonlinear input-output transfer function, near firing threshold. In particular, for nonlinear transfer functions with a sharp rise near threshold, the transition to chaos disappears in the limit of a large network; instead, the system exhibits chaotic fluctuations even for small synaptic gain. Finally, we investigate transition to chaos in network models with spiking dynamics. We show that when synaptic time constants are slow relative to the mean inverse firing rates, the network undergoes a transition from fast spiking fluctuations with constant rates to a state where the firing rates exhibit chaotic fluctuations, similar to the transition predicted by rate-based dynamics. Systems with finite synaptic time constants and firing rates exhibit a smooth transition from a regime dominated by stationary firing rates to a regime of slow rate fluctuations. This smooth crossover obeys scaling properties, similar to crossover phenomena in statistical mechanics. The theoretical results are supported by computer simulations of several neuronal architectures and dynamics. Consequences for cortical circuit dynamics are discussed. These results advance our understanding of the properties of intrinsic dynamics in realistic neuronal networks and their functional consequences.

Popular Summary

The firing patterns of neurons in the central nervous system often exhibit strong temporal irregularity and spatial heterogeneity. Previous studies have suggested that these properties are the outcome of the intrinsic chaotic dynamics of neural circuits, and simplified models of large neuronal networks with random Gaussian connections have been shown to exhibit a transition from a stationary state to slow, chaotic fluctuations. However, the transition to chaos in more realistic neuronal circuits has not yet been established. Here, we prove the universality of the transition to chaos using a mathematical framework for a broad class of networks with realistic connectivity architectures; we reveal the dependence of this transition on the form of the nonlinear input-output transfer function and synaptic time constants.

We theoretically study the onset of chaos using rate-based dynamics of sparsely connected networks composed of multiple neuronal subpopulations distinguished by their type (excitatory or inhibitory) and by the gain of their connectivities. We find that emergent dynamic balance between mean excitatory and inhibitory inputs amplifies the spatiotemporal fluctuations in the system, giving rise to a transition to chaos similar to that in a single population with Gaussian connectivity. Importantly, the sharpness of the neuronal transfer function has a crucial impact on the onset of chaos. Furthermore, in spiking networks, a sharp onset of chaos exists only in the limit of a large synaptic time constant. Scaling analysis, borrowed from critical phenomena, provides a quantitative description of the crossover from stationary rates to fluctuating rates for biologically realistic parameters.

We expect that our results will advance our understanding of the intrinsic dynamics of neuronal circuits in biologically relevant settings.

Figure 1

Network schematics. Arrows (circles) denote excitatory (inhibitory) connections. (a) The general network architecture studied here consists of multiple excitatory ($E$) and inhibitory ($I$) subpopulations driven by external inputs. Individual connections are randomly diluted or are Gaussian distributed. Subpopulations are distinguishable by the different statistics of connectivity. (b) The simplest model consists of a single population with random inhibitory recurrent connections with an external excitatory input. (c) Two randomly connected $E$ and $I$ populations.

Figure 2

Phase diagrams for a single inhibitory population. (a) No stable fixed point for sharp transfer functions. Simulation results for the normalized variance $q(\infty )$ for networks with $g=0.2$ and $g=0.1$ and different values of the power law, $\nu$, characterizing the rise of the transfer function, Eq. (33). Below $\nu =0.5$, the variance is nonzero, implying that there are temporal fluctuations even for low values of $g$. Simulation is performed on an inhibitory network ($N=5000$) with randomly diluted connections and mean activity of $m=0.1$. (b) Phase-space diagram for a threshold-linear ($\nu =1$) network of inhibitory neurons with mean connectivity $\overline{g}<0$ and variance $g$. The dashed line shows an example of a transition in a diluted network (with $K=650$) where, because $\overline{g}=-\sqrt{K}g$ [see Eq. (49)], the network lies on a line in the phase diagram. (c) Phase diagram for a threshold-quadratic ($\nu =2$) network with an excitatory external field $h^0=1$. The dashed line shows uniform instability of the fixed point, and the solid line shows the transition to the chaotic phase (shaded area). The dotted (black) line shows the instability of the chaotic phase found by simulations (Gaussian connectivity, $N=6000$, $h^0=1$). The inset (d) is the same as (c) for larger values of $|\overline{g}|$. The dashed line shows the existence for a diluted network as in (b).

Figure 3

Calculation of the largest Lyapunov exponent for a threshold-linear network. (a) Numerical integration of Eq. (45). The value of $x$ is determined through the requirement that the maximum of the potential at nonzero $q$ equals its value at zero $q$. Compare the form of the potential with the exact value of $x$ (solid line) with those calculated with $x$ deviating by $\pm 1\%$ from the correct value (dashed lines). (b) Numerical solution for the normalized variance $q(\infty )$ as a function of $g$. (c)–(f) The normalized autocorrelation function $q(\tau )$ found by integrating the equation of motions (45) using the correct value of $x$ for two values of $g$ (top), and the corresponding quantum potential $W=-{\partial }^{2}V/\partial q^2$ (bottom). The ground-state energies ${\epsilon }_0$ (horizontal lines) were found numerically. In both cases, they are negative, implying a positive Lyapunov exponent.

Figure 4

Numerical simulations of a single population with a linear-threshold transfer function. (a) Distribution of the local fields $h^i$ in the fixed-point regime. The thick curve shows normal distribution given the theoretical mean and variance obtained from the mean field theory. The dashed (red) vertical line marks the firing threshold, which is taken to be zero. (b,c) Activities of two networks with gains $g=2.1$ and $g=2.8$, respectively. Bold lines show spatially averaged local fields; thin lines are local fields of a sample of four neurons. (d) The normalized ac function, $\mathrm{\Delta }(\tau )/{\mathrm{\Delta }}_0=1-q(\tau )$, for two values of gain parameters as in (b) and (c). Solid black lines are the solutions of the equation of DMFT, Eq. (45), superimposed on simulation results (dashed lines). We show simulations performed on a network ($N=6800$) with Gaussian-distributed connections and $\overline{g}=\sqrt{K}g$, where $K=680$ and the external field $h^0=1$.

Figure 5

Robustness to changes in the synaptic distribution. (a) Autocorrelation function for sparse ($p=0.05$), dense ($p=0.8$), and Gaussian connectivity networks ($N=7000$, $h^0=1$). The gains were chosen such that all networks have the same variance for the connectivity distribution ($g=2.8$ for the Gaussian distribution and $g\sqrt{1-p}=2.8$ for the randomly diluted networks). (b) Normalized autocorrelation function $1-q(\tau )=\mathrm{\Delta }(\tau )/{\mathrm{\Delta }}_0$ compared with the solution of the DMFT.

Figure 6

Transition to chaos in a system of two populations. Parameters used: $J_{EE}=J_{IE}=W_E=\alpha g$, $J_{II}=-g$, $J_{EI}=-1.11g$, and $W_I=0.44g$. Note that $g$ is a global gain multiplying all the synapses; $\alpha$ denotes the excitatory efficacy. (a) The critical value of $g$ as a function of the excitatory efficacy $\alpha$. For $\alpha =0$, the network is identical to the single inhibitory network and the phase transition occurs at $g=\sqrt{2}$. (b) Largest eigenvalue of the stability matrix (53) for a network with $\alpha =0.55$ [marked by “${}^{^}$” in panel (a)]. MF predicts a phase transition when ${\mathrm{\Lambda }}_{+}=1$ (asterisks). (c) Normalized variance $q_k(\infty )$ for inhibitory (blue triangles) and excitatory (red circles) populations as a function of the gain calculated from network simulations with $\alpha =0.55$. It was calculated from simulation of a balanced diluted network with $N_E=N_I=3500$, $K=700$, connectivity parameters as in panel (b), and external activity $m_0=1$. Theory predicts a phase transition at $g=1.21$ (asterisks), as seen in panel (b).

Figure 7

Fluctuations of an E-I balanced network with linear-threshold transfer function. Same parameter set as in Fig. 6. (a) Traces of the local fields of six neurons from excitatory (red) and inhibitory (blue) populations. Dashed lines show mean input into each population. (b) Time-lagged autocorrelation function of two populations computed from the simulation. (c) Trace of the fluctuations in the spatially averaged activities $\delta m_k(t)\equiv m_k(t)-{⟨m_k(t)⟩}_t$ of both populations [same color codes as (a)]. Simulations of network are similar to the one used in Fig. 6, with $g=1.6$.

Figure 8

Critical behavior of a single inhibitory population with linear transfer function. Normalized variance (a) network relaxation time (b), and largest Lyapunov exponent (c) as a function of the distance from the critical point $\epsilon =\frac{1}{2}{g}^2-1$. Circles show the average over 20 simulations (Gaussian connectivity, $N=6000$, $h^0=1$), and dashed lines show theoretical predictions (see Appendix pp5). (d) Rescaled one-dimensional quantum potential $F(\tau /\sqrt{\epsilon )}$, Eq. (63), for the Hamiltonian in Eq. (29).

Figure 9

Critical behavior of a single inhibitory population with quadratic ($\nu =2$) transfer function. Normalized variance (a), network relaxation time (b), and LE (c) for small $\epsilon \equiv g^2/g_c^2-1$ gain above the critical point. Circles show the average over 20 simulations (Gaussian connectivity in the balance regime, $N=6000$, $K=1200$, $h^0=1$), and dashed lines show theoretical predictions (see Appendix pp5).

Figure 10

Critical behavior of a network with two populations and a linear transfer function. Diamond (red) and circles (blue) denote excitatory and inhibitory populations, respectively. (a) The normalized variances of the two populations, $q_E(\infty )$ and $q_I(\infty )$, show the same critical behavior (a linear rise from zero) as the single population [Eq. (62)]. Because of the linearity, the normalized autocorrelations are nearly equal to each other near the critical point (see Appendix pp5). (b) Network relaxation times of both populations are similar at the scaling limit. (c) LE of the entire network calculated over all populations, Eq. (27), scale as a single population, Eq. (64). We show simulations conducted on a network with a diluted connectivity matrix ($N=3600$, $K=600$ for each population) and a parameter set as in Figs. 7 and 6.

Figure 11

Network of spiking neurons with inhomogeneous Poisson statistics. Simulation results for a network that follows the dynamic equation in Eq. (65). (a) Above the critical point of the underlying rate dynamics, $g_c=\sqrt{2}$, the rate term dominates the DMF equation, Eq. (72), and the autocorrelation calculated from simulation ($r{\tau }_s=30$, solid line) is fully explained by autocorrelations of the rates, found by Eq. (45) (dashed line). (b) For values of $g$ well below the rate transition, the rates do not have intrinsic fluctuations, and the fluctuations in the fields are solely due to the Poisson firing. MFT results of Eq. (74) (dashed line) agree with simulations ($r{\tau }_s=30$, solid line). In intermediate values of $g$, rate and spike fluctuations interact nonlinearly. (c) Normalized variance of the local fields $q(\infty )$ from simulations of networks at different gain levels in the transition region. For higher levels of $r{\tau }_s$, the crossover between variance dominated by the Poisson process and variance dominated by the fields dynamics becomes sharper. (d) The coefficient of variation ($CV$), Eq. (76), for different values of gain. For low values of $g$, the $CV$ approaches 1, indicating a Poisson process with (almost) constant rates. For higher values of the gain, the $CV$ is larger, as expected from an inhomogeneous Poisson process. The crossover between the two regimes becomes sharper as $r{\tau }_s$ increases. Simulations are performed using an inhibitory population with Gaussian-distributed connectivity ($N=3500$, $K=700$) and $\overline{g}=h^0=1$.

Figure 12

Scaling behavior near criticality. (a) Rescaled form for the normalized variance, $(r{\tau }_s{)}^{\alpha }q(\infty )$ and (b) inverse effective relaxation time $(r{\tau }_s{)}^{\gamma }{\tau }_{\mathrm{eff}}^{-1}$, using the scaling variable $z={\tau }_{s}r|\epsilon {|}^{\beta }\mathrm{sign}(\epsilon )$. The data reveal the underlying scaling functions $f(z)$ and $F(z)$, respectively, consistent with Eqs. (77) and (80). Dashed lines show the predicted asymptotic behavior (see text). Simulations are performed on a Gaussian network of various sizes ($N=2000$, 3500, and 6000, $h^0=1$) and with various synaptic time scales (see legend). Effective relaxation was calculated from simulations by taking ${\tau }_{\mathrm{eff}}=\int \tau \delta \mathrm{\Delta }(\tau )d\tau /\int \delta \mathrm{\Delta }(\tau )d\tau$ for $\tau \ge 0$.