Bistability and criticality in the stochastic Wilson-Cowan model

We study a stochastic version of the Wilson-Cowan model of neural dynamics, where the response function of neurons grows faster than linearly above the threshold. The model shows a region of parameters where two attractive fixed points of the dynamics exist simultaneously. One fixed point is characterized by lower activity and scale-free critical behavior, while the second fixed point corresponds to a higher (supercritical) persistent activity, with small fluctuations around a mean value. When the number of neurons is not too large, the system can switch between these two different states with a probability depending on the parameters of the network. Along with alternation of states, the model displays a bimodal distribution of the avalanches of activity, with a power-law behavior corresponding to the critical state, and a bump of very large avalanches due to the high-activity supercritical state. The bistability is due to the presence of a first-order (discontinuous) transition in the phase diagram, and the observed critical behavior is connected with the line where the low-activity state becomes unstable (spinodal line).

Figure 1

(a) Activation functions $f\left(s\right)$ for $\gamma =0,1,2,3,4$ (from bottom to top). (b) Plot of the two terms of Eq. (2), $\alpha \mathrm{\Sigma }$ (dashed line) and $(1-\mathrm{\Sigma })f\left(s\right)$ (solid curves), where $s=w_0\mathrm{\Sigma }+h$, for $\alpha =\beta =0.1{\mathrm{ms}}^{-1}$, $w_0=0.9$, $h=0$, and $\gamma =0,1,2,3,4$ (from bottom to top). Attractive fixed points correspond to the crossing points where the solid curve has a lower derivative than the dashed line. For $\gamma <2.33$ there is only one attractive fixed point at $\mathrm{\Sigma }=0$, while for $\gamma \ge 2.33$ another fixed point appears (discontinuously) at higher values of $\mathrm{\Sigma }$. Inset: Zoom of the crossing of the curves at $\mathrm{\Sigma }=0$. (c) Values of $\mathrm{\Sigma }$ at the attractive fixed points, for $w_0=0.9$, $h=0$, and $0<\gamma <4$. At $\gamma =2.33$ a fixed point with $\mathrm{\Sigma }>0$ appears discontinuously.

Figure 2

Comparison between the dynamics of the network simulated by (a) the continuous-time Markov process and (b) the Langevin equations (1), for $\gamma =2$, $w_0=1$, $N={10}^4$ (dark blue), and $N={10}^7$ (light cyan). In (a) the state of individual neurons is simulated using the event-driven Gillespie algorithm, while in (b) one uses the Gaussian noise approximation to simulate directly the variables $k$ and $l$, representing the number of active excitatory and inhibitory neurons. While for $N={10}^4$ fluctuations are larger and the system switches frequently between different attractive fixed points, for $N={10}^7$ fluctuations are small and the switching is very unlikely (no switching is observed in the time window of the figure).

Figure 3

Mean field phase diagram for $h\to 0$ and $\alpha =\beta$. The red dashed line is the continuous transition between the phase of low activity where ${\mathrm{\Sigma }}_0=0$, and the phase of high activity where ${\mathrm{\Sigma }}_0>0$. The lower blue solid line corresponds to the points where a solution ${\mathrm{\Sigma }}_0>0$ appears discontinuously, and the upper blue solid line to the points where the solution ${\mathrm{\Sigma }}_0=0$ becomes unstable. Between the blue solid lines there are two attractive fixed points of the network dynamics.

Figure 4

Probability distribution of the firing rate $R$ for $N_E=N_I={10}^4$, on the critical line $w_0=1$ and different values of $\gamma$.

Figure 5

Firing rate $R\left(t\right)$ as a function of time for $N_E=N_I={10}^4$, on the critical line $w_0=1$ and (from left to right) $\gamma =1,2,3$. The red dashed line marks the value of the minimum in probability distribution of $R$ (see Fig. 4). Lower row: Firing rate $R\left(t\right)$ with an enlarged scale to put in evidence the short bursts of activity that characterize the critical low-activity state.

Figure 6

Probability distribution of the lifetime of [(a), (c)] low-activity and [(b), (d)] high-activity states, for $N={10}^4$, $2\times {10}^4$, $w_0=1$, and different values of $\gamma$.

Figure 7

(a) Lifetimes of low- and high-activity states for $N={10}^4$, $w_0=1$, as a function of $\gamma$. (b) Lifetimes of low- and high-activity states for $\gamma =3$, $w_0=1$, as a function of the number of neurons. (c) The same for $\gamma =4$, $w_0=0.84$.

Figure 8

Hysteresis loops for a system of $N={10}^5$ neurons, with (a) $\gamma =0$ and (b) $\gamma =4$. The parameter $w_0$ increases from $w_0=0.7$ to $w_0=1.1$ in 50 s, and then decreases back to the value $w_0=0.7$ in another 50 s. The lower red (upper blue) line shows the firing rate per neuron when $w_0$ is increased (decreased).

Figure 9

Avalanche distributions, defined as histograms of the probability that an avalanche has a size or duration within some interval, divided by the width of the interval (logarithmic spacings). The parameters are taken on the second-order transition line and upper spinodal, $w_0=1$, with $N={10}^4$ and $\gamma$ between 0 and 4. The width of the time bins is $\delta =1$ ms. (a) Size distribution $P\left(S\right)$. The exponent was measured in the case $\gamma =3$ for sizes $50<S<5000$. (b) Duration distribution $P\left(T\right)$. The exponent was measured in the case $\gamma =3$ for durations $5<T<300$. (c) Mean size $\langle S\rangle \left(T\right)$ as a function of the duration of the avalanche. (d) Dependence of $P\left(S\right)$ on the number of neurons, for $w_0=1$, $\delta =1$ ms, $\gamma =3$, and $N=100$ (leftmost) to $N=14000$ (rightmost). (e) Dependence of $P\left(S\right)$ on the width of time bins $\delta$, for $N={10}^4$, $w_0=1$, $\gamma =3$. Inset: Power-law exponent $\alpha$ measured for avalanches with size $50<S<5000$. (f) Dependence of $P\left(T\right)$ on the width of time bins $\delta$, for $N={10}^4$, $w_0=1$, $\gamma =3$. Exponents were measured with the powerlaw package [35].

Figure 10

Avalanche distributions on the lower spinodal, $w_0=w\left(\gamma \right)$, defined with a finite threshold in the number of spikes in a time bin, $n_{\text{th}}=N\delta R_0$, where $R_0$ is the firing rate per neuron at the high-activity fixed point. (a) Size distribution $P\left(S\right)$ for a system with $N={10}^4$, $\gamma =2,3,4$, and the width of the time bins is $\delta =1$ ms. The exponent was measured in the case $\gamma =3$ for sizes $50<S<5000$. (b) Duration distribution $P\left(T\right)$. The exponent was measured in the case $\gamma =3$ for durations $5<T<100$.

Figure 11

Average normalized shape $V\left(t\right)/V_0$ of avalanches having durations in a small interval around $T$, for (a) $\gamma =0$ and (b) $\gamma =3$.