Rings of delay-coupled neurons possess a striking capability to produce various stable spiking patterns. In order to reveal the mechanisms of their appearance, we present a bifurcation analysis of the Hodgkin–Huxley (HH) system with delayed feedback as well as a closed loop of HH neurons. We consider mainly the effects of external currents and communication delays. It is shown that typically periodic patterns of different spatial form (wavenumber) appear via Hopf bifurcations as the external current or time delay changes. The Hopf bifurcations are shown to occur in relatively narrow regions of the external current values, which are independent of the delays. Additional patterns, which have the same wavenumbers as the existing ones, appear via saddle–node bifurcations of limit cycles. The obtained bifurcation diagrams are evidence for the important role of communication delays for the emergence of multiple coexistent spiking patterns. The effects of a short-cut, which destroys the rotational symmetry of the ring, are also briefly discussed.
Networks of delay-coupled systems play an important role in modelling many real-world processes and experimental systems, including interacting lasers [1–6], real neuronal networks [7–9], networks of electronic systems [10,11] and traffic models [12,13]. They can be used for information processing [14,15], secure communication [16–18], random number generators  and other applications.
It has been shown that closed loops of unidirectionally connected neurons are capable of producing a large variety of stable spiking patterns [20–22] with appropriately adjusted propagation time delays or synaptic weights between the nodes. The well-known repressilator, a model of a synthetic gene network , represents a simple (topologically) closed loop with three nodes. The closed loops are also generic components in the nervous system  and naturally appear in the context of excitable cardiac tissue and orientation tuning in the visual cortex [25,26]. The above-mentioned examples show how important cycles (closed directed loops) are for determining dynamics in networks. Among other things, cycles provide feedback as well as round-trip time delays, which naturally lead to the appearance of multiple attractors, an important ingredient for the processing of information in neuronal systems.
It has been reported [22,27] that the dynamics of a single system with delayed feedback and a ring of coupled systems with delay have similar features. Therefore, in this paper, we perform a bifurcation analysis of the following two set-ups: a single Hodgkin–Huxley (HH) neuron with delayed feedback, and a ring of unidirectionally coupled HH neurons. In the ring, the state of a single neuron is described by the vector and evolves accordingly to the system 1.1where the nonlinear function Fj represents the HH equations [28,29] for a single neuron synaptically coupled to its nearest neighbour: Here, the variable V j, j=1,…,N, models the membrane potential of neuron j. The neurons are excitatory coupled (V R=0) via the postsynaptic potentials sj . The gating variables are mj, hj and nj. The last neuron in the ring is coupled to the first one, s1≡sN+1. The parameters K and τ determine the coupling strength and connection delays, respectively. Note that all delays can be assumed to be identical without loss of generality . In this paper, we also assume that all coupling weights K are identical, for the sake of simplicity. The rate functions of the gate variables are determined as follows: The remaining parameters are V Na=50, V K=−77, V L=−54.4 mV; gNa=120, gK=36, gL=0.3 mS cm−2; C=1 μF cm−2 and K=2 mS cm−2 if not otherwise stated. Time and delays are assumed to be in milliseconds.
The structure of the paper is as follows. In §2, we perform the analysis of a single neuron with delayed feedback. We start with the stability analysis of the equilibrium corresponding to the membrane resting potential. Depending on the value of the external current density I, the equilibrium may or may not undergo Hopf bifurcations, leading to the appearance of periodic solutions. We determine the corresponding admissible values of I. Although the Hopf bifurcations are subcritical, and the periodic solutions appear unstable, they stabilize further through the secondary bifurcations by varying I. We also investigate in detail the emergence of periodic solutions with changing time delays τ. Section 3 is devoted to similar questions in a ring of coupled neurons. Contrary to the single HH neuron, the ring may exhibit coexisting periodic patterns with different spatial organization (wavenumbers). We study the appearance of such solutions versus τ and I. Finally, §4 briefly discusses the effect of symmetry breaking, i.e. when one link is added, making a short-cut between two non-neighbouring nodes. In this case, the analysis becomes more complicated, since the rotational symmetry is no longer present.
2. Single neuron with delayed feedback
In this section, we investigate the properties of a single neuron with delayed feedback.
(a) Steady state and its stability
The membrane resting potential can be found as the solution of the scalar equation 2.1where , , and are steady-state values of the gating variables and the postsynaptic potential given by For the considered parameter values, the scalar equation (2.1) can be easily studied numerically. It exhibits a unique root corresponding to the steady state . The bifurcation diagram in §2c shows this steady state depending on the values of current density I as a thick black line.
In the following, we study the stability of the steady state and its bifurcations. The eigenvalues are given as solutions of the characteristic equation where is the 5×5 identity matrix, , and . Specifically, we have where
Although the time delay that we consider is not very large, we will use in the following an asymptotic technique, which is developed for systems with large delays [31–35]. It appears that this technique will give quite useful and easily representable information about the eigenvalues of the steady state in our case for moderate delays.
We start with the ansatz 2.2where we introduce two new real variables γ and ω instead of one complex λ. Then, the characteristic equation reads In the case when τ is large, the term can be neglected (see rigorous details in ), and we obtain The latter equation can be solved explicitly with respect to e−γ−iωτ: 2.3Finally, by taking the absolute values of both sides, we eliminate τ and find γ as a function of ω: 2.4The function γ(ω) plays an important role in determining the eigenvalues. In particular, the curve in the complex plane λ(ω)=γ(ω)/τ+iω contains leading eigenvalues for large τ. This curve is known as the asymptotic continuous spectrum (ACS) [31,35,36]. The ACS γ(ω) also provides useful information in the case of moderate delays, because the bifurcation points are given exactly by the zeros of ACS γ(ω)=τ Re λ=0. Strictly speaking, there may exist another part of the spectrum, the strong unstable spectrum, which has different asymptotics from (2.2), but this spectrum typically does not lead to bifurcations.1
Figure 1 (left panel) illustrates the curves of ACS (lines) as well as numerically computed eigenvalues (circles) using the DDE-Biftool software  for different values of the current density I. The following main features can be observed:
1. Destabilization of the stable steady state at I=8.4 μA cm−2 (figure 1a–c, left panel). One observes that the steady state is stable for I<8.4 μA cm−2, where all eigenvalues have negative real parts (figure 1a). With the increasing of I over 8.4 μA cm−2, the steady state destabilizes (figure 1b,c). Moreover, this destabilization is ruled by ACS. In spite of the fact that the exact eigenvalues are computed for τ=30 ms, this destabilization point is practically independent of the delay, because it is determined by the parameter values at which the ACS is tangent to the imaginary axis. Indeed, for larger time delays, the eigenvalues are just more densely distributed on the curve, thus still having negative real parts for I<8.4 μA cm−2. The destabilization of steady state leads to the Hopf bifurcation and the appearance of periodic spiking solutions.
2. Region where the steady state is unstable and ACS is unstable for I∈IH1=[8.4,10.5] (figure 1c–e, left panel). With further increasing of I, ACS develops a cusp (figure 1d) and an eigenvalue, which is unstable and does not belong to ACS, emerges. This is the strong unstable eigenvalue, which is present in figure 1e–i. Since the ACS is unstable, changing delays τ may lead to Hopf bifurcations, because the eigenvalues move along the curve and are able to cross the imaginary axis.
3. Steady state is unstable, and ACS is stable for I∈[10.5,143.7] (figure 1g, left panel). Instability of the steady state is ensured by the strong unstable eigenvalues. In this interval, no Hopf bifurcations are possible with varying time delay. Indeed, the eigenvalues cannot cross the imaginary axis, because ACS is not crossing it.
4. Similar to the case 2, the steady state and ACS are unstable for I∈IH2=[143.7,156.4]. Here, Hopf bifurcations are possible with varying delays. Within this region, the strong unstable eigenvalue disappears, becoming a part of the ACS.
5. Stabilization of the steady state at I=156.4 μA cm−2. For I>156.4 μA cm−2, the steady state is again stable.
As follows from the analysis of eigenvalues and ACS, the Hopf bifurcations are possible in relatively narrow regions of the current density I∈IH1∪IH2=[8.4,10.5]∪[143.7,156.4]. It will be shown later that the families of periodic spiking patterns emerge in a Hopf bifurcation within the first interval and terminate at a Hopf bifurcation within the second interval. This holds independently of the delay τ, even though the number of such families depends strongly on τ. Section 2b illustrates this in more detail.
(b) Hopf bifurcation branches
For small fixed delay τ (including τ=0), there is one Hopf bifurcation within each of the intervals IH1 and IH2. Using the numerical bifurcation software DDE-Biftool , we followed these Hopf bifurcations along the parameter τ and obtained the bifurcation curves in the (τ,I) parameter plane (figure 2). One can clearly observe that the branches are confined to the stripes I∈IH1∪IH2, as predicted by the theory in §2a. The typical ‘periodic-like’ structure of the branches can be explained by the following arguments: Each Hopf bifurcation at some delay value τ0 corresponds to the birth of a trivial periodic solution with zero amplitude and some period , where λ0 is the critical eigenvalue. The same periodic solution reappears in the system with time delay τ0+lT(τ0) with arbitrary integer l (it can be checked by substitution, see also ). Thus, the Hopf bifurcation reappears for all delay values τ0+lT(τ0) as well. Using this rule, in our case, one can continue the Hopf branch analytically for all values τ>12 ms from a given numerically computed primary branch for 0≤τ<12 ms. It also explains the ‘periodic-like’ structure in figure 2. For more detail concerning this reappearance, we refer to . It follows, in particular, that the number of Hopf bifurcations grows with increasing τ. More specifically, the number of Hopf bifurcations (i.e. the number of cross sections of the Hopf curves with the vertical line at τ) grows linearly with τ.
(c) Bifurcation diagrams
The bifurcation diagrams in figure 3 show the branches of periodic solutions, which emanate from the Hopf bifurcations. One can observe that the typical branch connects two Hopf points. For realistic time delays up to 60 ms, multiple Hopf bifurcations, which necessarily occur with increasing τ, do not occur yet. As a result, we observe only one branch of periodic solutions. Despite this, the coexistence of multiple periodic solutions does develop via another mechanism, which can be clearly observed in figure 3. With the increase in the delay, additional loops are created on the branch. The appearance and growth of the first loop is shown in detail in figure 3 for τ changing from 9 to 15 ms. For even larger delays, more and more loops are created. Additionally, one observes that some part of the loops becomes stable (solid lines in figure 3), whereas other parts are unstable (dashed lines). The number of such loops gives the maximal number of stable coexisting periodic solutions. Note that such a structure of branches leads to multiple saddle–node bifurcations of limit cycles with changing I.
Coexisting periodic solutions differ mainly by their period and basically represent different harmonics of the delay interval. Figure 4 shows profiles of periodic solutions for τ=40 ms and I=−10 μA cm−2. One observes the coexistence of stable solutions with approximately two (figure 4a), four (figure 4b,c) and six (figure 4d,e) spikes per delay interval. Note also that there are multiple solutions with the same number of spikes, as in (b) and (c), or (d) and (e).
3. Ring of N unidirectionally coupled Hodgkin–Huxley neurons
(a) Stability of the homogeneous steady state
Let us consider the ring of N unidirectional coupled, identical HH neurons. In this system, there is a unique homogeneous steady state , where was obtained in §2a. The characteristic equation, which governs the stability of the steady state, reads 3.1In a similar way, as we did in §2a, by making the ansatz λ=γ/τ+iω and neglecting terms that are proportional to γ/τ, we obtain an equation for ACS: 3.2The matrix in (3.2) is circulant. Thus, it can be diagonalized by a discrete Fourier transform . Therefore, the characteristic equation can be factorized  as 3.3where φl=2πl/N, l=1,…,N. The solution to equation (3.3) reads where Y (ω) is explicitly given by (2.3). By taking the absolute values, we obtain N identical solutions 3.4where γ(ω) is the ACS for the single system with delayed feedback.
Hence, the ACS for the ring of unidirectionally coupled systems consists of identical curves coinciding with the ACS for the single system. Note that this result is independent of the specific system, i.e. of the matrices A and B. In our specific case, there exist N identical branches of ACS for the ring (3.4). Note that this result has been shown earlier for the case of Stuart–Landau oscillators .
What are the consequences for the spectrum of eigenvalues of the steady state, especially in comparison with that of the single system? First, the qualitative properties of the spectrum are the same as in the single system. More exactly, the destabilization occurs at I≈8.4 μA cm−2, the ACS is unstable for I∈IH1 leading to possibilities of Hopf bifurcation, etc. In short, the conclusions 1–5 from §2a hold for the ring as well. Second, owing to the overlapping of N=5 curves of ACS, there are five times more eigenvalues located on the curve for the same delay value. Owing to the more dense location of eigenvalues, more Hopf bifurcations occur within the intervals IH1 and IH2, leading to the appearance of multiple branches of periodic solutions. The numerically computed typical spectrum, along with ACS, is shown in figure 1 (right panel). A comparison of the right and left panels in figure 1 shows the similarity of the spectrum of the single system with feedback and the ring.
(b) Bifurcation diagrams: coexisting rotating waves
We have seen that a single HH neuron with feedback exhibits two Hopf bifurcations as the current density I is varied and τ is small or intermediate. This has led to the appearance of one branch of periodic solutions in the diagrams in figure 3. In a ring consisting of N neurons, there exist N times more eigenvalues on the ACS, as shown in §3a. Consequently, there occur at least 2N Hopf bifurcations as I is varied, leading to at least N branches of spiking periodic solutions. Although for very large values of τ additional Hopf bifurcations appear, the corresponding values of τ are unrealistically large and we concentrate on the intermediate values of τ, where just N branches exist.
Numerical continuation of the system with N=5, τ=20 ms and K=3 mS cm−2 shows (figure 5a) how five branches of spiking periodic solutions originate within the theoretically predicted interval IH1 and terminate within the interval IH2. Furthermore, one observes the appearance of loops on every branch, similar to the case of the single HH neuron (figure 3). As a rule, every loop on each branch leads to a pair of stable and unstable solutions.
In order to demonstrate the multistability and identify the type of the spiking solutions, we illustrate in figure 5b all found periodic solutions for the fixed value I=10 μA cm−2. We plot the temporal frequency f of the solutions versus k, the number of spiking fronts, which is defined as k=N(φi+1−φi)/2π=N(ti+1−ti)/T, where φi+1−φi is the phase difference and ti+1−ti is the time difference between the nearest spiking times of the neighbouring neurons. Similar coexistence diagrams have been obtained by a direct integration [20,21]. In our case, for these fixed parameter values, we observe 32 coexisting spiking solutions, 14 of which are stable (diamonds in figure 5b). Some of them have the same number of firing fronts k, i.e. the same spatial structure along the ring, which is representing the space, but differ by their temporal frequency. By analysing the origin of the coexisting solutions within the diagram in figure 5a, we discover that the spiking solutions of the same spatial structure (k) originate from the same branch, and their coexistence comes from the looping of the branches, i.e. the saddle–node bifurcations along the parameter I. The solutions with different spatial structure originate from different branches.
4. Symmetry breaking by a short-cut
In this section, we briefly discuss the effect of a short-cut that is introduced in the ring. As an example, we consider an additional link between the first and fourth neurons, as shown in figure 6. The strength of the additional connection Ksc is considered as a control parameter. In particular, for Ksc=0, we start from the symmetric situation without a short-cut. Our goal is to investigate how non-zero Ksc changes the dynamics, in particular, how it influences the periodic spiking solutions. The time delay along the short-cut is τsc.
We observed that increasing Ksc may effect crucially the existing periodic solutions. Depending on the short-cut position, the branches of periodic solutions may destabilize, or even terminate, as the short-cut strength increases. Yet, some of the solutions remain stable. Figure 7 show examples of such stable branches and their properties. The phase differences of the neighbouring neurons Δφj,j−1=2π(tj−tj−1)/T are plotted versus Ksc, where tj are the spiking times (time moment when the corresponding voltage variables reach their maximum). For Ksc=0, the rotational symmetry implies that the phase differences are the same for all j: Δφj,j−1=Δφ=2πk/N. Such solutions are also called rotating waves. As a reference solution, we selected k=−1, hence Δφ=−2π/5. With the increase of the weight of the short-cut connection Ksc, we observe that the selected solution remains stable and its properties are modified. More specifically, the phase differences form clusters as follows: the difference Δφ54 (dashed line in figure 7) splits from the group and all other phase differences remain in the cluster (solid lines in figure 7). The deviation of Δφ54 can be explained by the fact that neuron 5 is additionally influenced by neuron 1. Interestingly, all other neurons still possess an almost constant phase shift with respect to their neighbours.
We have performed a bifurcation analysis of a closed unidirectional ring of coupled HH neurons with delayed connections. Stability analysis of the membrane resting potential state shows that its destabilization is possible through the Hopf bifurcations for values of the external current density I within two relatively narrow intervals IH1 and IH2. We have shown that the branches of periodic solutions emanate from the steady state in the interval IH1 and terminate in the interval IH2 as the external current density is varied. For a single HH system with delayed feedback, there is only one such branch for physically plausible values of delays, whereas for a ring of N coupled systems, there are N branches. Different branches are shown to contain spiking periodic solutions with different wavenumbers, stable as well as unstable. Moreover, each of the branches changes its geometrical form for intermediate delays, by developing loops. These loops lead to the coexistence of multiple stable spiking solutions with the same wavenumber. In this way, a large number of coexisting spiking solutions emerge in the closed loops of HH neurons, a phenomenon that was recently reported .
This work was supported by DFG in the framework of the Collaborative Research Center SFB910.
We acknowledge fruitful discussions with Leonhard Lücken, Jan Philipp Pade and Michael Zaks.
One contribution of 15 to a Theme Issue ‘Dynamics, control and information in delay-coupled systems’.
↵1 Note that γ(ω) describes a continuous set instead of the discrete spectrum, which delay differential equations typically possess . This follows from the fact that we considered only the absolute value of equation (2.3), instead of solving the complex equation (2.3). Further, the discrete values of ω can be obtained by solving another part of equation (2.3), for example, by considering its argument.
- © 2013 The Author(s) Published by the Royal Society. All rights reserved.