Dynamics of a large system of spiking neurons with synaptic delay

We analyze a large system of heterogeneous quadratic integrate-and-fire (QIF) neurons with time delayed, all-to-all synaptic coupling. The model is exactly reduced to a system of firing rate equations that is exploited to investigate the existence, stability, and bifurcations of fully synchronous, partially synchronous, and incoherent states. In conjunction with this analysis we perform extensive numerical simulations of the original network of QIF neurons, and determine the relation between the macroscopic and microscopic states for partially synchronous states. The results are summarized in two phase diagrams, for homogeneous and heterogeneous populations, which are obtained analytically to a large extent. For excitatory coupling, the phase diagram is remarkably similar to that of the Kuramoto model with time delays, although here the stability boundaries extend to regions in parameter space where the neurons are not self-sustained oscillators. In contrast, the structure of the boundaries for inhibitory coupling is different, and already for homogeneous networks unveils the presence of various partially synchronized states not present in the Kuramoto model: Collective chaos, quasiperiodic partial synchronization (QPS), and a novel state which we call modulated-QPS (M-QPS). In the presence of heterogeneity partially synchronized states reminiscent to collective chaos, QPS and M-QPS persist. In addition, the presence of heterogeneity greatly amplifies the differences between the incoherence stability boundaries of excitation and inhibition. Finally, we compare our results with those of a traditional (Wilson Cowan-type) firing rate model with time delays. The oscillatory instabilities of the traditional firing rate model qualitatively agree with our results only for the case of inhibitory coupling with strong heterogeneity.

Figure 1

Phase diagram for identical neurons, $\mathrm{\Delta }=0$. Shaded region: The asynchronous state (${\mathbf{a}}_{+}$) is stable. Slantwise hatched region: full synchrony is unstable. Horizontally hatched region: The fully synchronized state does not exist and the only attractor is the global rest state ${\mathbf{q}}_{-}$. The orbit of fully synchronized self-sustained oscillations is created at the dashed black line (at $\overline{\eta }<0$), Eq. (14). Blue and red lines are the loci of the sub- and super-critical Hopf-like instabilities of incoherence Eqs. (10). Solid green line: saddle-node bifurcation. The vertical dashed green line separates the oscillatory from the excitable regime of the QIF neuron.

Figure 2

Macroscopic (columns 1 and 2) and microscopic (columns 3–5) dynamics of QPS (rows a and b), M-QPS (row c) and collective chaos (row d), see Table 1. Column 1: time series of the mean firing rate. Blue lines correspond to numerical simulations of the FRE Eqs. (9), while red dotted lines are obtained computing the mean firing rate of a population of $N=2000$ QIF neurons. Column 2: $(r,v)$ phase portraits, numerically obtained using Eqs. (9). In panel (b2) two coexisting periodic attractors are shown: QPS-asym(I) (solid) and QPS-asym(II) (dashed)—see also inset of Fig. 3. Panels (b1, b3–b5) correspond to QPS-asym(I). Columns (3–5) show the dynamics of a population of $N=2000$ QIF neurons. Column 3: raster plots. Neurons are ordered according to their firing time at the beginning of the simulation (due to the first order nature of the QIF model, this ordering is preserved in time). Columns 4 and 5 show return ISI plots and ISI distributions of an arbitrary neuron $j$. The return plots of panels (a4, b4) are closed curves, indicating quasiperiodic microscopic dynamics in the QPS-sym and QPS-asym. The corresponding ISI histograms (a5, b5) show two (QPS-sym) or three (QPS-asym) peaks. In the M-QPS, neurons are quasiperiodic with three characteristic frequencies—the return plots of panel (c4) is a closed surface in 3D, and therefore its projection in 2D fills a defined region of the space. Parameters: (row a) $J=-9.2$ (row b) $J=-9.5$, (row c) $J=-10.3$, (row d) $J=-10.6$. We use $\sqrt{\overline{\eta }}=3.6$ in all simulations.

Figure 3

Four largest Lyapunov exponents for two alternative bifurcation sequences in a range of negative $J$ values and fixed $\sqrt{\overline{\eta }}=3.6$. For each solution, the continuation was carried out either increasing or decreasing parameter $J$ adiabatically. In the top panel the vertical dashed lines indicate, from right to left: a supercritical Hopf bifurcation (SC-H), a transcritical bifurcation (TC), a Neimark-Sacker bifurcation (NS), and a subcritical Hopf bifurcation (SB-H). In the bottom panel the vertical dashed lines indicate, from right to left: a saddle-node bifurcation (SN), a Neimark-Sacker bifurcation (NS), and the onset of chaos (C). The inset shows a sketch of the bifurcation diagram connecting the two bifurcation sequences.

Figure 4

Poincaré sections of the FRE (9) for $\sqrt{\overline{\eta }}=3.6$, and for three different values of the inhibitory coupling strength: (a) $J=-10.3$; (b) $J=-10.48$; (c) $J=-10.6$. The Poincaré surface is $v=0,\dot{v}<0$.

Figure 5

Numerical exploration of the partially synchronized states (QPS,M-QPS, collective chaos) near the supercritical Hopf bifurcation in phase diagram Fig. 1. In the light gray region the largest Lyapunov exponent is zero, and QPS is stable. The purple dots correspond to two vanishing Lyapunov exponents, indicating quasiperiodic dynamics. In the cyan region the dynamics is chaotic. The vertical dashed black line at $\sqrt{\overline{\eta }}=3.6$ corresponds to the range of parameters explored in Fig. 3.

Figure 6

Phase diagram for populations of heterogeneous neurons, $\mathrm{\Delta }=0.1$. Dark shaded region: Incoherence (fixed point) is the only stable state. Light shaded region: Incoherence (fixed point) coexist with a partially synchronous state (limit cycle). Brown region: Two forms of asynchrony (high and a low activity fixed points) coexist with a partially synchronous state. Green lines are saddle-node bifurcations, and black lines correspond to Hopf boundaries. Note that here, in contrast with Fig. 1, the Hopf boundaries are not represented in Blue/Red (we do not explicitly specify whether these boundaries are subcritical or supercritical). The boundary between light and dark shaded regions was obtained numerically.

Figure 7

Enlarged view of the region of multistability located at $\overline{\eta }<0$ in Fig. 6. Black line: Hopf bifurcation (subcritical). Green lines: saddle-node bifurcations. In the dark shaded region, only the quiescent, low activity state is stable. In the light shaded region, incoherence coexists with a collective oscillatory state—self-sustained due to recurrent excitation. In the brown region the low activity fixed point coexists with a high activity fixed point and with the oscillatory state. In the small dark purple region only the two high and low activity fixed points are attracting. Right panels: Sketches of the Poincaré section in different regions (assuming that it coincides with the one-dimensional manifold that connects different fixed points). The thick lines indicate two-dimensional manifolds, and periodic orbits are indicated by a point surrounded by a small circle.

Figure 8

Macroscopic (columns 1 and 2) and microscopic (columns 3–5) dynamics of (row a) PS-II states, (row b) M-PS states, and (row c) collective chaos for heterogeneous neurons—see Fig. 2 and Table 1. Column 1: Time series of the firing rate for the FRE Eqs. (9) (blue) and for a population of $N=2000$ QIF neurons Eqs. (1) (red dotted). Column 2 shows the corresponding attractors, obtained using the FRE. In rows (a) and (b), the dynamics is periodic but, in contrast with the identical case, here the limit cycle is asymmetric due to the presence heterogeneity. Column 3 shows the raster plots corresponding to numerical simulations of a population of $N=2000$ QIF neurons Eqs. (1), and columns 4 and 5 show the corresponding return plots and ISI histograms, respectively. In the raster plots, each neuron index $j$ corresponds to a specific ${\eta }_j$ value (see text). For the computation of return plots and ISI histograms we used neuron $j=500$. In panels (a4) and (b4) one can see that the neuron behaves quasiperiodically, with two and three incommensurable frequencies, respectively. Note also the three peaks in panel (a5) due to the asymmetry of the limit cycle. In all panels we use $\mathrm{\Delta }=0.1$, $\sqrt{\overline{\eta }}=3.5$, and (row a) $J=-9.60$; (row b) $J=-10.70$; (row c) $J=-11.30$.

Figure 9

Time-averaged coupling-modified ISIs as a function of the intrinsic ISI for a population of 2000 QIF neurons in three different states: (a) PS-II, (b) M-PS, and (c) collective chaos. The red dots are obtained with direct simulations of the population of QIF neurons, while the blue line is obtained forcing each neuron with the FRE. Note the multiple plateaus in the middle panel. Parameters are as in Fig. 8: $\sqrt{\overline{\eta }}=3.5$ and (a) $J=-9.60$; (b) $J=-10.70$; (c) $J=-11.30$.

Figure 10

The four largest Lyapunov exponents for $\mathrm{\Delta }=0.1$ and $\sqrt{\overline{\eta }}=3.5$. The stability regions of the different attractors are indicated by vertical gray dashed lines.

Figure 11

Increasing the level of heterogeneity $\mathrm{\Delta }$ reveals different synchronization scenarios for excitation and inhibition (see text). Black, dark blue, blue and light blue lines correspond, respectively, to the Hopf boundaries of Eqs. (9) with $\mathrm{\Delta }=0.1,5,10,20$. These boundaries determine the regions of stability of the incoherent/asynchronous states. In the shaded regions incoherence is stable for $\mathrm{\Delta }=0.1$. In the dark shaded region the only attractor is incoherence.

Figure 12

Oscillations emerge only for inhibitory coupling in the traditional firing rate model Eq. (17). In the gray region, limited by the black line ($\mathrm{\Delta }=0.1$), the fixed point determined by Eq. (16) is stable and looses stability via a Hopf bifurcation—compare with Fig. 11. The dark blue and blue curves correspond to $\mathrm{\Delta }=5\text{and}10$, respectively. The green dashed boundary corresponds to the case $\mathrm{\Delta }=0$ and is a straight line.