Subsections

• 9.3.1 Traveling fronts and waves (*)
• 9.3.2 Stability (*)

9.3 Patterns of spike activity

We have seen that the intricate interplay of excitation and inhibition in locally coupled neuronal nets can result in the formation of complex patterns of activity. Neurons have been described by a graded-response type formalism where the `firing rate' is given as a function of the `average membrane potential'. This approach is clearly justified for a qualitative treatment of slowly varying neuronal activity. In the context of spatio-temporal patterns of neuronal activity, however, a slightly closer look is in order.

In the following we will dismiss the firing rate paradigm and use the Spike Response Model instead in order to describe neuronal activity in terms of individual action potentials. We start with a large number of SRM neurons arranged on a two-dimensional grid. The synaptic coupling strength w of neurons located at $\vec{{r}}_{i}^{}$ and $\vec{{r}}_{j}^{}$ is, as hitherto, a function of their distance, i.e., w = w(|$\vec{{r}}_{i}^{}$ - $\vec{{r}}_{j}^{}$|). The response of a neuron to the firing of one of its presynaptic neurons is described by a response function $\epsilon$ and, finally, the afterpotential is given by a kernel named $\eta$, as customary. The membrane potential of a neuron located at $\vec{{r}}_{i}^{}$ is thus

$$\vec{{r}}_{i}^{} \int_{0}^{\infty} \eta \vec{{r}}_{i}^{} \sum_{j}^{} \vec{{r}}_{i}^{} \vec{{r}}_{j}^{} \int_{0}^{\infty} \epsilon \vec{{r}}_{j}^{}$$ (9.44)

with S($\vec{{r}}_{i}^{}$, t) = $\sum_{f}^{}$$\delta$(t - t_i^f) being the spike train of the neuron at $\vec{{r}}_{i}^{}$. Note that we have neglected distance-dependent propagation delays; constant (synaptic) delays, however, can be absorbed into the response function $\epsilon$. Note also that for the sake of notational simplicity we have included the afterpotential of all spikes and not only the afterpotential of the last spike. This, however, will not affect the present discussion because we restrict ourselves to situations where each neurons is firing only once, or with long inter-spike intervals.

Spikes are triggered whenever the membrane potential reaches the firing threshold $\vartheta$. This can be expressed in compact form as

$$\vec{{r}}_{i}^{} \delta \vec{{r}}_{i}^{} \vartheta \left[\vphantom{ \frac{\partial u(\vec{r}_i,t)}{\partial t} }\right. {\frac{{\partial u(\vec{r}_i,t)}}{{\partial t}}} \left.\vphantom{ \frac{\partial u(\vec{r}_i,t)}{\partial t} }\right]_{+}^{}$$ (9.45)

Here, [...]₊ denotes the positive part of its argument. This factor is required in order to ensure that spikes are triggered only if the threshold is crossed with a positive slope and to normalize the Dirac $\delta$ functions to unity.

Figure 9.11 shows the result of a computer simulation of a network consisting of 1000×1000 SRM neurons. The coupling function is mexican-hat shaped so that excitatory connections dominate on small and inhibitory connections on large distances. In a certain parameter regime the network exhibits an excitable behavior; cf. Section 9.2.2. Starting from a random initial configuration, a cloud of short stripes of neuronal activity evolves. These stripes propagate through the net and soon start to form rotating spirals with two, three or four arms. The spirals have slightly different rotation frequencies and in the end only a few large spirals with three arms will survive.

Figure 9.11: Snapshots from a simulation of a network of 1000×1000 neurons. [Taken from Kistler et al. (1998)].

Let us try to gain an analytic understanding of some of the phenomenon observed in the computer simulations. To this end we suppose that the coupling function w is slowly varying, i.e., that the distance between two neighboring neurons is small as compared to the characteristic length scale of w. In this case we can replace in Eq. (9.45) the sum over all presynaptic neurons by an integral over space. At the same time we drop the indices that label the neurons on the grid and replace both h and S by continuous functions of $\vec{{r}}\,$ and t that interpolate in a suitable way between the grid points for which they have been defined originally. This leads to field equations that describe the membrane potential u($\vec{{r}}\,$, t) of neurons located at $\vec{{r}}\,$,

$$\vec{{r}}\, \int_{0}^{\infty} \eta \vec{{r}}\, \int \vec{{r}}{^\prime} \vec{{r}}\, \vec{{r}}{^\prime} \int_{0}^{\infty} \epsilon \vec{{r}}{^\prime}$$ (9.46)

together with their spike activity S($\vec{{r}}\,$, t),

$$\vec{{r}}\, \delta \vec{{r}}\, \vartheta \left[\vphantom{ \frac{\partial u(\vec{r},t)}{\partial t} }\right. {\frac{{\partial u(\vec{r},t)}}{{\partial t}}} \left.\vphantom{ \frac{\partial u(\vec{r},t)}{\partial t} }\right]_{+}^{}$$ (9.47)

as a function of time t. In the following sections we try to find particular solutions for this integral equation. In particular, we investigate solutions in the form of traveling fronts and waves. It turns out, that the propagation velocity for fronts and the dispersion relation or waves can be calculated analytically; cf. Fig. 9.12. The stability of these solutions is determined in Section 9.3.2.

The approach sketched in Sections 9.3.1 and 9.3.2 (Kistler et al., 1998; Bressloff, 1999; Kistler, 2000) is presented for a network of a single population of neurons, but it can also be extended to coupled networks of excitatory and inhibitory neurons (Golomb and Ermentrout, 2001). In addition to the usual fast traveling waves that are found in purely excitatory networks, additional slow and non-continuous `lurching' waves appear in an appropriate parameter regime (Golomb and Ermentrout, 2001; Rinzel et al., 1998).

9.3.1 Traveling fronts and waves (*)

We start our analysis of the field equations (9.47) and (9.48) by looking for a particular solution in form of a plane front of excitation in a two-dimensional network. To this end we make an ansatz for the spike activity

$$\delta$$ (9.48)

This is a plane front that extends from y = - $\infty$ to y = + $\infty$ and propagates in positive x-direction with velocity v. We substitute this ansatz into the expression for the membrane potential Eq. () and find

$$\eta \int \int \left[\vphantom{\sqrt{(x-x')^2+(y-y')^2} }\right. \sqrt{{(x-x')^2+(y-y')^2}} \left.\vphantom{\sqrt{(x-x')^2+(y-y')^2} }\right] \epsilon$$ (9.49)

Up to now the propagation velocity is a free parameter. This parameter can be fixed by exploiting a self-consistency condition that states that the membrane potential along the wave front equals the firing threshold; cf. Eq. (9.48). This condition gives a relation between the firing threshold $\vartheta$ and the propagation velocity v, i.e.,

$$\int \int \left[\vphantom{\sqrt{(v\,t-x')^2+(y-y')^2} }\right. \sqrt{{(v\,t-x')^2+(y-y')^2}} \left.\vphantom{\sqrt{(v\,t-x')^2+(y-y')^2} }\right] \epsilon \overset{!}{=} \vartheta$$ (9.50)

Note that the afterpotential $\eta$ drops out because each neuron is firing only once. The propagation velocity as a function of the firing threshold is plotted in Fig. 9.12A. Interestingly, there are two branches that correspond to two different velocities at the same threshold. We will see later on that not all velocities correspond to stable solutions.

Figure 9.12: A. Propagation velocity v of a single plane front versus the firing threshold $\vartheta$. B. Dispersion relation of a traveling wave for different values of the firing threshold ( $\vartheta$ = 0, 0.01,..., 0.06; $\vartheta$ = 0 is the outermost curve, $\vartheta$ = 0.06 corresponds to the curve in the center).

The simulations show that the dynamics is dominated in large parts of the net by a regular pattern of stripes. These stripes are, apart from the centers of the spirals, formed by an arrangement of approximatively plane fronts. We can use the same ideas as above to look for such a type of solution. We make an ansatz,

$$\sum_{{n=-\infty}}^{{\infty}} \delta \left(\vphantom{t-\frac{x-n\,\lambda}{v} }\right. {\frac{{x-n\,\lambda}}{{v}}} \left.\vphantom{t-\frac{x-n\,\lambda}{v} }\right)$$ (9.51)

that describes a traveling wave, i.e., a periodic arrangement of plane fronts, with wave length $\lambda$ traveling in positive x direction with (phase) velocity v. Both the phase velocity and the wave length are free parameters that have to be fixed by a self-consistency condition. We substitute this ansatz in Eq. (9.47) and find

$$\sum_{{n=-\infty}}^{{\infty}} \eta \left(\vphantom{t-\frac{x-n\,\lambda}{v} }\right. {\frac{{x-n\,\lambda}}{{v}}} \left.\vphantom{t-\frac{x-n\,\lambda}{v} }\right) \left(\vphantom{x,y,t-\frac{x-n\,\lambda}{v} }\right. {\frac{{x-n\,\lambda}}{{v}}} \left.\vphantom{x,y,t-\frac{x-n\,\lambda}{v} }\right)$$ (9.52)

with

$$\int \int \left[\vphantom{\sqrt{(x-x')^2+(y-y')^2} }\right. \sqrt{{(x-x')^2+(y-y')^2}} \left.\vphantom{\sqrt{(x-x')^2+(y-y')^2} }\right] \epsilon$$ (9.53)

Using the fact that the membrane potential on each of the wave fronts equals the firing threshold we find a relation between the phase velocity and the wave length. This relation can be reformulated as a dispersion relation for the wave number k = 2$\pi$/$\lambda$ and the frequency $\omega$ = 2$\pi$ v/$\lambda$. The dispersion relation, which is shown in Fig. 9.12B for various values of the firing threshold, fully characterizes the behavior of the wave.

9.3.2 Stability (*)

A single front of excitation that travels through the net triggers a single action potential in each neuron. In order to investigate the stability of a traveling front of excitation we introduce the firing time t($\vec{{r}}\,$) of a neuron located at $\vec{{r}}\,$. The threshold condition for the triggering of spikes can be read as an implicit equation for the firing time as a function of space,

$$\vartheta \int \vec{{r'}}\, \vec{{r}}\, \vec{{r}}{^\prime} \epsilon \vec{{r}}\, \vec{{r}}{^\prime}$$ (9.54)

In the previous section we have found that t(x, y) = x/v can satisfy the above equation for certain combinations of the threshold parameter $\vartheta$ and the propagation velocity v.

We are aiming at a linear stability analysis in terms of the firing times (Bressloff, 1999). To this end we consider a small perturbation $\delta$t(x, y) which will be added to the solution of a plane front of excitation traveling with velocity v in positive x-direction, i.e.,

$$\delta$$ (9.55)

This ansatz will be substituted in Eq. (9.55) and after linearization we end up with a linear integral equation for $\delta$t, viz.

$$\int \vec{{r}}{^\prime} \vec{{r}}\, \vec{{r}}{^\prime} \epsilon{^\prime} \delta \vec{{r}}\, \delta \vec{{r}}{^\prime}$$ (9.56)

Due to the superposition principle we can concentrate on a particular form of the perturbation, e.g., on a single Fourier component such as $\delta$t(x, y) = e^$\lambda$ x cos$\kappa$ y. The above integral equation provides an implicit relation between the wave number $\kappa$ of the perturbation with the exponent $\lambda$ that describes its growth or decay as the front propagates through the net. If there is for a certain $\kappa$ a particular solution of

$$\begin{multline} 0 = \int \!\!\!\! \int \!\! {\text{d}}x' \,{\text{d}}y' \, w... ...\cos (\kappa\,y) - {\text{e}}^{\lambda\,x'}\,\cos (\kappa\,y')] \end{multline}$$

with Re($\lambda$) > 0 then the front is unstable with respect to that perturbation.

Figure 9.13 shows the result of a numerical analysis of the stability equation (9.58). It turns out that the lower branch of the v-$\vartheta$ curve corresponds to unstable solutions that are susceptible to two types of perturbation, viz., a perturbation with Im($\lambda$) = 0 and a oscillatory perturbation with Im($\lambda$)$\ne$ 0. In addition, fronts with a velocity larger than a certain critical velocity are unstable because of a form instability with Im($\lambda$) = 0 and $\kappa$ > 0. Depending on the actual coupling function w, however, there may be a narrow interval for the propagation velocity where plane fronts are stable; cf. Fig. 9.13B.

Figure 9.13: Similar plots as in Fig. 9.12A showing the propagation velocity v as a function of the firing threshold $\vartheta$ (thick line). Dashes indicate solutions that are unstable with respect to instabilities with $\kappa$ = 0. These instabilities could also be observed in a one-dimensional setup. Thin dotted lines mark the domain of instabilities that show up only in two dimensions, i.e., form instabilities with $\kappa$ > 0. Wide dots correspond to perturbations with Im($\lambda$) = 0 and narrow dots to Hopf-instabilities with Im($\lambda$) > 0. A. Coupling function w(r) = $\sum_{{i=1}}^{2}$a_iexp(- r²/$\lambda_{i}^{2}$) with $\lambda_{1}^{2}$ = 15, $\lambda_{2}^{2}$ = 100, a₁ = 1.2 and a₂ = - 0.2. B. Similar plot as in A, but with a₁ = 1.1 and a₂ = - 0.1. Note that there is a small interval between v = 0.475 and v = 0.840 that corresponds to stable fronts. [Taken from Kistler (2000)].

The stability of plane waves can be treated in a similar way as that of a plane front, we only have to account for the fact that each neuron is not firing only once but repetitively. We thus use the following ansatz for the firing times {t_n(x, y)| n = 0,±1,±2,...} of a neuron located at (x,y),

$${\frac{{x-n\,\lambda}}{{v}}} \delta$$ (9.57)

with $\delta$t_n(x, y) being a `small' perturbation. If we substitute this ansatz into Eq. (9.53) we obtain

$$\begin{multline} 0 = \sum_n \eta'\left [\frac{(n-m)\,\lambda}{v} \right] \, [... ...bda}{v} \right ] \, [\delta t_m(\vec{r})-\delta t_n(\vec{r}')] \end{multline}$$

in leading order of $\delta$t_n and with $\vec{{r}}\,$ = (x, y), $\vec{{r}}{^\prime}$ = (x', y'). Primes at $\eta$ and $\epsilon$ denote derivation with respect to the argument. This equation has to be fulfilled for all $\vec{{r}}\,$ $\in$ $\mathbb {R}$² and all m $\in$ $\mathbb {Z}$.

For the sake of simplicity we neglect the contribution of the after potential in Eq. (9.60), i.e., we assume that $\eta{^\prime}$[n $\lambda$/v] = 0 for n > 0. This assumption is justified for short-lasting afterpotentials and a low firing frequency.

As before, we concentrate on a particular form of the perturbations $\delta$t_n(x, y), namely $\delta$t_n(x, y) = exp[c (x - n $\lambda$)] cos($\kappa_{n}^{}$ n) cos($\kappa_{y}^{}$ y). This corresponds to a sinusoidal deformation of the fronts in y-direction described by $\kappa_{y}^{}$ together with a modulation of their distance given by $\kappa_{n}^{}$. If we substitute this ansatz for the perturbation in Eq. (9.60) we obtain a set of equations that can be reduced to two linearly independent equations for c, $\kappa_{n}^{}$, and $\kappa_{y}^{}$. The complex roots of this system of equations determines the stability of traveling waves, as it is summarized in Fig. 9.12.

Figure 9.14: A. Dispersion relation (solid lines) of periodic wave trains for various values of the threshold parameter $\vartheta$ = 0 (uppermost trace), 0.1, ..., 0.6 (center). Network parameters as in Fig. 9.13A. The shading indicates solutions that are unstable with respect to perturbations with $\kappa_{n}^{}$ $\in$ {0,$\pi$}. The unshaded (white) region is bounded to the left by a form instability with $\kappa_{n}^{}$ = 0 and $\kappa_{y}^{}$ > 0 (dots). Its right border (dashes) is formed by a Hopf bifurcation with $\kappa_{n}^{}$ = $\pi$ and $\kappa_{y}^{}$ > 0. The right branch of this curve corresponds to the same type of bifurcation but with $\kappa_{n}^{}$ = 0. The remaining dotted lines to the left and to the right indicate form instabilities with $\kappa_{n}^{}$ = $\pi$. The long-dashed lines reflect bifurcations with $\kappa_{y}^{}$ = 0, i.e., bifurcations that would show up in the corresponding one-dimensional setup. B. Similar plot as in A but with a coupling function as in Fig. 9.12B. The firing threshold varies from $\vartheta$ = 0 (leftmost trace) to $\vartheta$ = 1.0 (center) in steps of 0.2. [Taken from Kistler (2000)].