Subsections

• • 8.1.0.1 Example: Phase diagram of instabilities
  • 8.1.0.2 Example: Oscillations in random networks

8.1 Instability of the Asynchronous State

In Section 6.4 and throughout Chapter 7, we have assumed that the network is in a state of asynchronous firing. In this section, we study whether asynchronous firing can indeed be a stable state of a fully connected population of spiking neurons - or whether the connectivity drives the network towards oscillations. For the sake of simplicity, we restrict the analysis to SRM₀ neurons; the same methods can, however, be applied to integrate-and-fire neurons or general SRM neurons.

For SRM₀ neurons, the membrane potential is given by

$$\eta \hat{{t}}_{i}^{}$$ (8.1)

where $\eta$(t - $\hat{{t}}_{i}^{}$) is the effect of the last firing of neuron i (i.e., the spike itself and its afterpotential) and h(t) is the total postsynaptic potential caused by presynaptic firing. If all presynaptic spikes are generated within the homogeneous population under consideration, we have

$$\sum_{j}^{} \sum_{f}^{} \epsilon_{0}^{} \int_{0}^{\infty} \epsilon_{0}^{}$$ (8.2)

Here $\epsilon_{0}^{}$(t - t_j^(f)) is the time course of the postsynaptic potential generated by a spike of neuron j at time t_j^(f) and w_ij = J₀/N is the strength of lateral coupling within the population. The second equality sign follows from the definition of the population activity, i.e., A(t) = N^-1$\sum_{j}^{}$$\sum_{f}^{}$$\delta$(t - t_j^(f)); cf. Chapter 6.1. For the sake of simplicity, we have assumed in Eq. (8.2) that there is no external input.

The state of asynchronous firing corresponds to a fixed point A(t) = A₀ of the population activity. We have already seen in Chapter 6.4 how the fixed point A₀ can be determined either numerically or graphically. To analyze its stability we assume that for t > 0 the activity is subject to a small perturbation,

$$\omega \lambda$$ (8.3)

with A₁ $\ll$ A₀. The perturbation in the activity induces a perturbation in the input potential,

$$\omega \lambda$$ (8.4)

with h₀ = J₀ $\hat{{\epsilon}}$(0) A₀ and h₁ = J₀ $\hat{{\epsilon}}$($\omega$ - i$\lambda$) A₁, where

$$\hat{{\epsilon}} \omega \lambda \hat{{\epsilon}} \omega \lambda \psi \omega \lambda \int_{0}^{\infty} \epsilon_{0}^{} \omega \lambda$$ (8.5)

is the Fourier transform of $\epsilon_{0}^{}$ and $\psi$(.) denotes the phase shift between h and A.

The perturbation of the potential causes some neurons to fire earlier (when the change in h is positive) others to fire later (whenever the change is negative). The perturbation may therefore build up ($\lambda$ > 0, the asynchronous state is unstable) or decay back to zero ($\lambda$ < 0, the asynchronous state is stable). At the transition between the region of stability and instability the amplitude of the perturbation remains constant ($\lambda$ = 0, marginal stability of the asynchronous state). These transition points, defined by $\lambda$ = 0, are determined now.

We start from the population integral equation A(t) = $\int_{{-\infty}}^{t}$P_I(t|$\hat{{t}}$) A($\hat{{t}}$) d$\hat{{t}}$ that has been introduced in Chapter 6.3. Here P_I(t|$\hat{{t}}$) is the input-dependent interval distribution, i.e., the probability density of emitting a spike at time t given that the last spike occured at time $\hat{{t}}$. We have seen in Chapter 7 that the linearized population activity equation can be written in the form

$$\Delta \int_{{-\infty}}^{t} \hat{{t}} \Delta \hat{{t}} \hat{{t}} {{\text{d}}\over {\text{d}}t} \int_{0}^{\infty} \mathcal {L} \Delta$$ (8.6)

cf. Eq. (7.3). Here P₀(t - $\hat{{t}}$) is the interval distribution during asynchronous firing and $\mathcal {L}$ is the kernel from Table 7.1. We use $\Delta$A(t) = A₁ e^{i$\omega$t+$\lambda$t} and $\Delta$h(t) = h₁ e^{i$\omega$t+$\lambda$t} in Eq. (8.6). After cancellation of a common factor A₁exp(i$\omega$t) the result can be written in the form

$$\omega {J_0\,A_0\, \hat{\epsilon}(\omega)\, \hat{{\mathcal{L}}}(\omega) \over 1 - \hat{P}(\omega)} \omega \Phi \omega$$ (8.7)

$\hat{{P}}$($\omega$) and $\hat{{{\mathcal{L}}}}$($\omega$) is the Fourier transform of the interval distribution P₀(t - $\hat{{t}}$) and the kernel $\mathcal {L}$, respectively. The second equality sign defines the real-valued functions S_f($\omega$) and $\Phi$($\omega$). Equation (8.7) is thus equivalent to

$$\omega \Phi \omega \pi$$ (8.8)

Solutions of Eq. (8.8) yield bifurcation points where the asynchronous firing state looses its stability towards an oscillation with frequency $\omega$.

We have written Eq. (8.8) as a combination of two requirements, i.e., an amplitude condition S_f($\omega$) = 1 and a phase condition $\Phi$($\omega$)  mod  2$\pi$ = 0. Let us discuss the general structure of the two conditions. First, if S_f($\omega$)$\le$1 for all frequencies $\omega$, an oscillatory perturbation cannot build up. All oscillations decay and the state of asynchronous firing is stable. We conclude from Eq. (8.7) that by increasing the absolute value | J₀| of the coupling constant, it is always possible to increase S_f($\omega$). The amplitude condition can thus be met if the excitatory or inhibitory feedback from other neurons in the population is sufficiently strong. Second, for a bifurcation to occur we need in addition that the phase condition is met. Loosely speaking, the phase condition implies that the feedback from other neurons in the network must arrive just in time to keep the oscillation going. Thus the axonal signal transmission time and the rise time of the postsynaptic potential play a critical role during oscillatory activity (Mattia and Del Giudice, 2001; Treves, 1993; Abbott and van Vreeswijk, 1993; Brunel, 2000; Neltner et al., 2000; Gerstner, 1995; Vreeswijk, 2000; Gerstner, 2000b; Brunel and Hakim, 1999; Ernst et al., 1995; Gerstner and van Hemmen, 1993; Tsodyks et al., 1993).

8.1.0.1 Example: Phase diagram of instabilities

Let us apply the above results to SRM₀ neurons with noise in the reset. We assume that neurons are in a state of asynchronous firing with activity A₀. As we have seen in Chapter 5, the interval distribution for noisy reset is a Gaussian centered at T₀ = 1/A₀. The filter function $\mathcal {L}$ is a $\delta$-function, $\mathcal {L}$(x) = $\delta$(x)/$\eta{^\prime}$; cf. Table 7.1. Hence Eq. (8.7) is of the form

$${i\omega\over \eta'} {J_0 \, A_0 \,\hat{\epsilon}(\omega) \over 1 - \hat{{\mathcal{G}}}_\sigma (\omega) \,e^{-i\omega T_0 }} \omega \Phi \omega$$ (8.9)

where $\hat{{{\mathcal{G}}}}_{\sigma}^{}$($\omega$) = exp{ - $\sigma^{2}_{}$$\omega^{2}_{}$/2} is the Fourier transform of a Gaussian with width $\sigma$.

In order to analyze Eq. (8.9) numerically we have to specify the response kernel. For the sake of simplicity we choose a delayed alpha function,

$$\epsilon_{0}^{} {s-\Delta^{\rm ax}\over \tau^2} \left(\vphantom{ - {s-\Delta^{\rm ax}\over \tau} }\right. {s-\Delta^{\rm ax}\over \tau} \left.\vphantom{ - {s-\Delta^{\rm ax}\over \tau} }\right) \mathcal {H} \Delta^{{\rm ax}}_{}$$ (8.10)

The Fourier transform of $\epsilon$ as defined in Eq. (8.5) has an amplitude |$\hat{{\epsilon}}$($\omega$)| = (1 + $\omega^{2}_{}$ $\tau^{2}_{}$)^-1 and a phase $\psi$($\omega$) = $\omega$ $\Delta^{{\rm ax}}_{}$ +2 arctan($\omega$ $\tau$). Note that a change in the delay $\Delta^{{\rm ax}}_{}$ affects only the phase of the Fourier transform and not the amplitude.

Figure 8.1 shows S_f as a function of $\omega$ T₀. Since S_f = 1 is a necessary condition for a bifurcation, it is apparent that bifurcations can occur only for frequencies $\omega$ $\approx$ $\omega_{n}^{}$ = n 2$\pi$/T₀ with integer n where T₀ = 1/A₀ is the typical inter-spike interval. We also see that higher harmonics are only relevant for low levels of noise. For $\sigma$$\to$ 0 the absolute value of the denominator of (8.9) is 2| sin($\omega$T₀/2)| and bifurcations can occur for all higher harmonics. At a high noise level, however, the asynchronous state is stable even with respect to perturbations at $\omega$ $\approx$ $\omega_{1}^{}$.

Figure 8.1: Amplitude condition for instabilities in the asynchronous state. The amplitude S_f is plotted as a function of the normalized frequency $\omega$ T₀ for two different values of the noise: $\sigma$ = 1ms (solid line) and $\sigma$ = 0.1ms (dashed line). Instabilities of the asynchronous firing state are possible at frequencies where S_f > 1. For low noise S_f crosses unity (dotted horizontal line) at frequencies $\omega$ $\approx$ $\omega_{n}^{}$ = n 2$\pi$/T₀. For $\sigma$ = 1ms there is a single instability region for $\omega$ T₀ $\approx$ 1. For the plot we have set T₀ = 2$\tau$.

A bifurcation at $\omega$ $\approx$ $\omega_{1}^{}$ implies that the period of the perturbation is identical to the firing period of individual neurons. Higher harmonics correspond to instabilities of the asynchronous state towards cluster states (Kistler and van Hemmen, 1999; Ernst et al., 1995; Gerstner and van Hemmen, 1993; Golomb et al., 1992; Golomb and Rinzel, 1994): each neuron fires with a mean period of T₀, but the population of neurons splits up in several groups that fire alternatingly so that the overall activity oscillates several times faster; cf. Section 8.2.3.

Figure 8.1 illustrates the amplitude condition for the solution of Eq. (8.9). The numerical solutions of the full equation (8.9) for different values of the delay $\Delta^{{\rm ax}}_{}$ and different levels of the noise $\sigma$ are shown in the bifurcation diagram of Fig. 8.2. The insets show simulations that illustrate the behavior of the network at certain combinations of transmission delay and noise level.

Let us consider for example a network with transmission delay $\Delta^{{\rm ax}}_{}$ = 2 ms, corresponding to a x-value of $\Delta^{{\rm ax}}_{}$/T₀ = 0.25 in Fig. 8.2. The phase diagram predicts that, at a noise level of $\sigma$ = 0.5 ms, the network is in a state of asynchronous firing. The simulation shown in the inset in the upper right-hand corner confirms that the activity fluctuates around a constant value of A₀ = 1/T₀ = 0.125 kHz.

If the noise level of the network is significantly reduced, the system crosses the short-dashed line. This line is the boundary at which the constant activity state becomes unstable with respect to an oscillation with $\omega$ $\approx$ 3 (2$\pi$/T₀). Accordingly, a network simulation with a noise level of $\sigma$ = 0.1 exhibits an oscillation of the population activity with period T^osc $\approx$ T₀/3 $\approx$ 2.6 ms.

Keeping the noise level constant but reducing the transmission delay corresponds to a horizontal move across the phase diagram in Fig. 8.2. At some point, the system crosses the solid line that marks the transition to an instability with frequency $\omega_{1}^{}$ = 2$\pi$/T₀. Again, this is confirmed by a simulation shown in the inset in the upper left corner. If we now decrease the noise level, the oscillation becomes even more pronounced (bottom right).

Figure 8.2: Stability diagram (center) for the state of asynchronous firing in a SRM₀ network as a function of noise $\sigma$ (y-axis) and delay $\Delta^{{\rm ax}}_{}$ (x-axis). Parameters are J₀ = 1 and $\tau$ = 4 ms. The threshold $\vartheta$ was adjusted so that the mean inter-spike interval is T₀ = 2$\tau$. The diagram shows the borders of the stability region with respect to $\omega_{1}^{}$,...,$\omega_{4}^{}$. For high values of the noise, the asynchronous firing state is always stable. If the noise is reduced, the asynchronous state becomes unstable with respect to an oscillation either with frequency $\omega_{1}^{}$ (solid border lines), or $\omega_{2}^{}$ (long-dashed border lines), $\omega_{3}^{}$ (short-dashed border lines), or $\omega_{4}^{}$ (long-short dashed border lines). Four insets show typical patterns of the activity as a function of time taken from a simulation with N = 1000 neurons. Parameters are $\sigma$ = 0.5 ms and $\Delta^{{\rm ax}}_{}$ = 0.2 ms (top left); $\sigma$ = 0.5 ms and $\Delta^{{\rm ax}}_{}$ = 2.0 ms (top right); $\sigma$ = 0.1 ms and $\Delta^{{\rm ax}}_{}$ = 0.2 ms (bottom left); $\sigma$ = 0.1 ms and $\Delta^{{\rm ax}}_{}$ = 2.0 ms (bottom right). Taken from (Gerstner, 2000b).

In the limit of low noise, the asynchronous network state is unstable for virtually all values of the delay. The region of the phase diagram in Fig. 8.2 around $\Delta^{{\rm ax}}_{}$/T₀ $\approx$ 0.1 which looks stable hides instabilities with respect to the higher harmonics $\omega_{6}^{}$ and $\omega_{5}^{}$ which are not shown. We emphasize that the specific location of the stability borders depends on the form of the postsynaptic response function $\epsilon$. The qualitative features of the phase diagram in Fig. 8.2 are generic and hold for all kinds of response kernels.

The numerical results apply to the response kernel $\epsilon_{0}^{}$(s) defined in (8.10) which corresponds to a synaptic current $\alpha$(s) with zero rise time; cf. (4.2) and (4.34). What happens if $\alpha$ is a double exponential with rise time $\tau_{{\rm rise}}^{}$ and decay time $\tau_{{\rm syn}}^{}$? In this case, the right-hand side of (8.9) has an additional factor [1 + i $\omega$ $\tau_{{\rm rise}}^{}$]^-1 which leads to two changes. First, due to the reduced amplitude of the feedback, instabilities with frequencies $\omega$ > $\tau_{{\rm rise}}^{{-1}}$ are suppressed. The tongues for the higher harmonics are therefore smaller. Second, the phase of the feedback changes. Thus all tongues of frequency $\omega_{n}^{}$ are moved horizontally along the x-axis by an amount $\Delta$/T₀ = - arctan($\omega_{n}^{}$ $\tau_{{\rm rise}}^{}$)/(n 2$\pi$).

What happens if the excitatory interaction is replaced by inhibitory coupling? A change in the sign of the interaction corresponds to a phase shift of $\pi$. For each harmonic, the region along the delay axis where the asynchronous state is unstable for excitatory coupling (cf. Fig. 8.2) becomes stable for inhibition and vice versa. In other words, we simply have to shift the instability tongues for each frequency $\omega_{n}^{}$ horizontally by an amount $\Delta$/T₀ = 1/(2n). Apart from that the pattern remains the same.

8.1.0.2 Example: Oscillations in random networks

Figure 8.3: Oscillations with irregular spike trains. The activity A(t) (bottom left) of a population of 5000 integrate-and-fire neurons exhibits oscillations which are not evident in the spike trains of 50 individual neurons (top left), but which are confirmed by a significant oscillatory component in the autocorrelation function (bottom right). The spike trains have a broad distribution of interspike intervals (top right); taken from Brunel and Hakim (1999).

Our discussion of random-connectivity networks in Chapter 6.4.3 has been focused on the stationary state of asynchronous firing. The stability analysis of the asynchronous state in such randomly connected networks is completely analogous to the approach sketched in Eqs. (8.3) - (8.6) except that the linearization is performed on the level of the density equations (Brunel, 2000; Brunel and Hakim, 1999). Close to the asynchronous state, the activity can be written as A(t) = A₀ + A₁(t) and the membrane potential distribution as p(u, t) = p₀(u) + p₁(u, t). Here p₀(u) is the stationary distribution of membrane potential in the state of asynchronous firing [cf. Eqs. (6.27) and (6.28)] and p₁(u, t) is a small time-dependent perturbation. The stability analysis requires a linearization of the Fokker-Planck equation (6.21) with respect to p₁ and A₁.

For short transmission delays, the asynchronous state A(t) $\equiv$ A₀ can loose its stability towards an oscillation with a frequency that is much faster than the single-neuron firing rate. Brunel (2000) distinguishes two different variants of such fast oscillations. First, as in the previous example there are cluster states where the neuronal population splits into a few subgroups. Each neuron fires nearly regularly and within a cluster neurons are almost fully synchronized; cf. Section 8.2.3. Second, there are synchronous irregular states where the global activity oscillates while individual neurons have a broad distribution of interspike intervals; cf. Fig. 8.3. We will come back to synchronous irregular states in Section 8.3.