Dynamics of globally delay-coupled neurons displaying subthreshold oscillations

Cristina Masoller, M. C. Torrent, Jordi García-Ojalvo

Abstract

We study an ensemble of neurons that are coupled through their time-delayed collective mean field. The individual neuron is modelled using a Hodgkin–Huxley-type conductance model with parameters chosen such that the uncoupled neuron displays autonomous subthreshold oscillations of the membrane potential. We find that the ensemble generates a rich variety of oscillatory activities that are mainly controlled by two time scales: the natural period of oscillation at the single neuron level and the delay time of the global coupling. When the neuronal oscillations are synchronized, they can be either in-phase or out-of-phase. The phase-shifted activity is interpreted as the result of a phase-flip bifurcation, also occurring in a set of globally delay-coupled limit cycle oscillators. At the bifurcation point, there is a transition from in-phase to out-of-phase (or vice versa) synchronized oscillations, which is accompanied by an abrupt change in the common oscillation frequency. This phase-flip bifurcation was recently investigated in two mutually delay-coupled oscillators and can play a role in the mechanisms by which the neurons switch among different firing patterns.

1. Introduction

The subthreshold activity of a neuron is receiving increasing attention because it can carry valuable information about the neuron’s synaptic properties and the coupling architecture of the network (Devor Yarom 2002; Chorev et al. 2007). It has been suggested that electrical coupling via gap junctions can serve to synchronize the rhythmic activity of the neurons (Long et al. 2002). From an information theoretic point of view, subthreshold oscillations can play an important role in the information transfer and signal transduction mechanisms in sensory neurons (Braun et al. 1994).

In neurophysiology, time delays naturally arise from the fact that electrical signals propagate realistically along neurons with a finite speed. At the single neuron level and in ensembles of coupled neurons, delays have been shown to generate instabilities, multi-stability, chaotic oscillations and oscillation death (Foss et al. 1996; Pakdaman et al. 1996; Vibert et al. 1998; Foss & Milton 2000; Sainz-Trapaga et al. 2004; Masoller et al. 2008). When the neurons are coupled through both the instantaneous and the delayed collective mean field, by varying the delay one can control (i.e. enhance or suppress) the synchrony of the neuronal rhythmic activity (Rosenblum & Pikovsky 2004a,b; Rosenblum et al. 2006). Purely noise-induced neuronal oscillations can also be controlled by the application of a delayed feedback loop (Hauschildt et al. 2006).

Here, we aim at addressing the interplay between the natural oscillatory subthreshold activity of a neuron and the delay time of a global mean field coupling. We consider an ensemble of neurons that are globally coupled via their time-delayed collective mean field. The individual neuron is modelled using a Hodgkin–Huxley-type conductance model with parameters chosen such that the uncoupled neuron displays subthreshold oscillations of the membrane potential (the unperturbed natural frequency being approx. 10 Hz). We show that the activity of the neuronal ensemble is controlled by the relation between two time scales: the natural period of the subthreshold oscillations of the uncoupled neuron and the delay of the mean field coupling. Depending on these time scales, various oscillatory regimes are possible: either all neurons fire spikes or display subthreshold activity, or there is clustering, with some neurons displaying subthreshold activity, while the others fire spikes. When there is synchronized activity, the neurons oscillate at the same frequency and can be in-phase or out-of-phase. In the latter case, the neuronal oscillations have a phase difference such that the collective mean field is constant or has small-amplitude oscillations. Also, depending on the relationship between the delay and the natural oscillation period, these regimes coexist (different initial conditions result in different patterns of neuronal activity), or only one of these regimes is possible (the oscillatory activity of the ensemble is the same, regardless of the initial condition).

We also show that the phase-shifted activity of the neuronal ensemble can be interpreted in terms of a phase-flip bifurcation occurring in two delay-coupled oscillators. At the bifurcation point, the synchronized dynamics changes from in-phase to out-of-phase oscillations (or vice versa), and this change is accompanied by an abrupt variation in the common oscillation frequency. The phase-flip bifurcation was originally reported by Schuster & Wagner (1989) in a simple model of two delay-coupled phase oscillators and was recently investigated in detail by Prasad et al. (2006, 2008) in various types of delay-coupled oscillators, in which not only the phase but also the amplitude of the oscillators are affected by the coupling. Prasad et al. considered models of excitable lasers, Fitzhugh–Nagumo neurons, predator–prey models and chaotic Chua electronic circuits, thus demonstrating the broad relevance of this bifurcation. We show here that the phase-flip bifurcation can play an important role in the mechanisms that mediate the switching among different neuronal patterns.

This paper is organized as follows. The model is presented in §2. Section 3 presents the results of simulations. In §4, we show that the behaviour of the delay-coupled neurons can be understood in terms of a simpler model, composed of a set of limit cycle oscillators with delayed mean field coupling. A phase-flip bifurcation is shown to occur when the delay is varied. A summary of the results is presented in §5.

2. Neuron model

The equation for the membrane potential of the ith neuron, Vi, with i=1,…,N is that of the model by Braun et al.(1998, 2000), extended to incorporate delayed mean field coupling: Embedded Image 2.1 where CM is the membrane capacitance, T the temperature, Il,i=g l(ViVl) a passive leak current, INa,i and IK,i the fast sodium and potassium currents and Isd,i and Isr,i additional slow currents. These four currents are given by Embedded Image 2.2 where k stands for Na, K, sd and sr, gk are the maximum conductances of the corresponding channels, and the activation variables ak are given by Embedded Image 2.3 for K and sd, Embedded Image 2.4 and Embedded Image. The steady-state activations are defined by Embedded Image 2.5 The temperature-dependent factors ρ(T) and ϕ(T) are Embedded Image 2.6

The last term on the right-hand side of equation (2.1) accounts for the delayed mean field coupling, with η representing the coupling strength and τ the delay time. The delay takes into account both the conduction delay associated with axonal propagation of electrical signals and that neurons have a finite reaction time to the inputs that they receive. We consider for simplicity that the coupling term is linear, an approach that has also been followed by other authors (e.g. Rosenblum & Pikovsky 2004a). We remark, however, that the behaviour of the ensemble is robust to the specific type of coupling considered, in the sense that more realistic coupling terms (e.g. Rosenblum & Pikovsky 2004b), which may be of chemical or electrical origin (and thus may occur via gap junctions or via synapses), result in similar effects (i.e. enhancement or suppression of synchrony depending on η and τ), as in the case of the simple linear coupling term. Moreover, we point out that equation (2.1) can be re-written, via a redefinition of the parameters of the passive leak current, gl and Vl, such as to write the coupling term in the form of Vi(t)−Vj(tτ), which would more realistically represent gap-junction-mediated coupling. We note that self-feedback circuits, representing recurrent delayed connections, are included in the coupling term in equation (2.1).

For simplicity, in this study, we consider that the coupling strength and the transmission delay time are the same in all neurons. Moreover, we consider that the neurons are identical. We focus on the deterministic dynamics and thus we do not take into account noise, the effect of which is left for future work.

3. Results

The parameters used in the simulations are those used by Braun et al. (1998) and Masoller et al. (2008): CM=1 μF cm−2, g l=0.1 mS cm−2, gNa=1.5 mS cm−2, gK=2 mS cm−2, gsd=0.25 mS cm−2, gsr=0.4 mS cm−2, Vl=−60 mV, VNa=50 mV, VK=−90 mV, Vsd=50 mV, Vsr=−90 mV, sNa=0.25 mV−1, sK=0.25 mV−1, ssd=0.09 mV−1, V0Na=−25 mV, V0K=−25 mV, V0sd=−40 mV, τK=2 ms, τsd=10 ms, τsr=20 ms, θ=0.17, μ=0.012, A1=1.3, A2=3.0, T=35oC and Tc=25oC. For these parameters, the uncoupled neuron displays subthreshold oscillations of period T0∼130 ms. To integrate equation (2.1), it is necessary to specify the initial value of the membrane potential, Vi, on the time interval [−τ,0]. We chose arbitrary initial conditions in such a way that the neurons are oscillating in their natural cycle before the coupling begins to act (i.e. the coupling starts when the neurons are at random phases of the oscillation cycle).

Figure 1 displays the amplitude of the oscillations of the mean field, Embedded Image 3.1 versus the delay of the global coupling for a fixed coupling strength (left: η=0.001; right: η=−0.001). The oscillation amplitude is defined as Embedded Image and the values of Embedded Image and Embedded Image are calculated in a time window of duration ΔT (typically, ΔT=200T0) after letting transients die away. The figure shows a substantial degree of multi-stability among different oscillation amplitudes for the mean field, mainly in the case η<0.

Figure 1.

Amplitude of the mean field oscillations. For each value of τ, results of simulations with 10 random initial conditions are presented. The coupling strength is η=0.001 (left column) and η=−0.001 (right column), and the number of neurons is N=2,3,5,10 from top to bottom for both the columns. For certain delay times, the neuronal ensemble exhibits multi-stability, i.e. there are coexisting stable solutions, with different basins of attraction, that result in more than one mean field oscillation amplitude (i.e. two or more dots for a given value of τ, e.g. τ/T0∼0.2 for N=3 and positive coupling). For other delays, there is monostability, and the oscillation amplitude of the mean field is the same, independent of the initial conditions. Multi-stability is seen for most delay values in the case of negative coupling (right column).

In order to visualize the dynamics at the single neuron level, the amplitude and the frequency of the oscillations of the individual neurons are plotted in figures 2 and 3, respectively. The amplitude of the oscillations of a single neuron is defined in the same way as the amplitude of the oscillations of the mean field: Embedded Image, with Embedded Image and Embedded Image calculated in a time window after transients die away. The frequency of the oscillations of the ith neuron is defined as that of the main peak in the Fourier spectrum of the membrane potential Vi(t). The size of the network increases from N=2 (top) to N=10 (bottom). We present results of simulations with 10 random initial conditions. The delay is limited to τT0 because of the range of biologically plausible values. Typically, the speed of signal propagation through axonal fibres is of the order of 1 m s−1, and can result in time delays up to 100 ms, including the delay owing to both dendritic and synaptic signal processing (Kandel et al. 1991; Dhamala et al. 2004). Typical periods of subthreshold oscillations are also of the order of 100 ms or smaller (Xing et al. 2001).

Figure 2.

Amplitude of the oscillations of the individual neurons. The coupling strength is η=0.001 (left column) and η=−0.001 (right column), and the number of neurons is N=2,3,5,10 from top to bottom for both the columns. The dashed line shows the natural oscillation amplitude of the uncoupled neuron. For positive coupling, all neurons display subthreshold oscillations, the amplitude of the individual oscillations being equal to or slightly lower than the natural oscillation amplitude, approximately 40 mV. For negative coupling, three different behaviours can be clearly distinguished: (i) all neurons fire spikes: this occurs in a delay window near τ/T0∼0.6, where it can be seen that the oscillation amplitude of all neurons is >60 mV; (ii) all neurons display subthreshold oscillations: this occurs in a delay window near τ/T0∼0.8, where it can be seen that the oscillation amplitude of all neurons is slightly >40 mV; and (iii) some neurons fire spikes while others display subthreshold oscillations: this occurs for most delay values, and it can be observed that the oscillation amplitude of some neurons is >60 mV, and of others is slightly >40 mV.

Figure 3.

Frequency of the oscillations of the individual neurons. The coupling strength is η=0.001 (left column) and η=−0.001 (right column), and the number of neurons is the same as in figures 1 and 2. An abrupt change in the frequency of the oscillations is seen at τ/T0∼0.5 (this change is more pronounced for negative coupling). The change in the frequency is accompanied by a variation in the phase relation of the neuronal oscillations, which results in an abrupt variation in the amplitude of the oscillations of the mean field, as can be seen in figure 1. It can be noticed in figure 1 the variation of the mean field amplitude near τ/T0=0.5 for both positive and negative coupling: for τ/T0<0.5, the neurons oscillate out-of-phase, which results in small mean field oscillation amplitude; for τ/T0>0.5, the neurons oscillate in-phase, which results in larger mean field oscillation amplitude.

For η>0, the individual neurons display only subthreshold oscillations, which can be either in-phase, partially out-of-phase or in perfect out-of-phase depending on τ. Typical time traces are displayed in figure 4. An almost constant mean field happens only for large N (N≥5, figure 4c); for smaller N, figure 4a,b, in spite of the fact that the neurons oscillate in a perfect out-of-phase relationship, because of the shape of the oscillation cycle, the mean field is not constant and displays small oscillations, the frequency of which is harmonic of the frequency of the oscillations of the individual neurons.

For η<0, the neurons display a more complex behaviour: depending on τ, either all the neurons fire spikes or display subthreshold oscillations, or some neurons display subthreshold oscillations while the others fire spikes (i.e. the ensemble divides into clusters). The neuronal oscillations can be either in-phase or out-of-phase, depending on the delay, as shown in figure 5. If the oscillations are in-phase, the amplitudes of the oscillations of the mean field and a single neuron are equal; if the oscillations are out-of-phase, the amplitude of the oscillations of the mean field is smaller than the amplitude of the oscillations of one neuron, as can be seen in figure 5.

Figure 4.

Neuronal oscillations for positive coupling: out-of-phase ((a) N=2,τ/T0= 0.45, (b) N=3,τ/T0=0.45, (c) N=5,τ/T0=0.45) and in-phase ((d) N=5, τ/T0= 0.55). Thick lines: mean field; thin lines: individual neurons. η=0.001.

Figure 5.

Neuronal oscillations for negative coupling: partial out-of-phase ((a) τ/T0=0.45), in-phase ((b) τ/T0=0.55 and (c) τ/T0=0.65) and perfect out-of-phase ((d) τ/T0=0.75). Thick lines: mean field; thin lines: individual neurons. N=5, η=−0.001.

For delay values such that the neuronal oscillations are in-phase, Vi(t)=Vj(t)=V (t) ∀i, j, and the dynamics is as that of a single neuron with a recurrent delayed connection, Embedded Image 3.2 which was studied in detail by Masoller et al. (2008).

For negative coupling, there are also delay values for which the neurons display perfect out-of-phase subthreshold oscillations that result in a nearly constant mean field if N is large enough. There is also multi-stability of solutions with the coexistence, for certain delay values, of in-phase and out-of-phase behaviour.

The variation of the oscillation frequency with the delay, displayed in figure 3, is such that when the oscillations are synchronized in-phase, the common frequency tends to decrease with τ, while when the oscillations are out-of-phase, the common frequency is locked, i.e. independent of τ.

4. Limit cycle model

Let us now show that the above-presented results can be qualitatively understood in terms of a simple model, of N limit cycle oscillators coupled through their delayed mean field: Embedded Image 4.1 where Zk is a complex variable and the coupling strength η is real.

Figure 6 displays the amplitude and the frequency of the oscillators (black dots) and of the mean field: Embedded Image (grey dots). It can be seen that there are delay values where the oscillators are in-phase, and there is a bifurcation where the phase relation among the oscillators changes abruptly, resulting in zero mean field. When the oscillators are synchronized in phase, the common frequency decreases with τ, while when there is the bifurcation to phase-shifted oscillations, the common frequency changes abruptly and remains locked, i.e. independent of τ. When the oscillations are phase shifted, the mean field displays small oscillations whose frequency is harmonic of the common frequency (in this case, in figure 6, the mean field frequency is not shown because of the vertical scale).

Figure 6.

Amplitude ((a) η=0.01, (b) η=−0.01) and frequency ((c) η=0.01, (d) η=−0.01) of the limit cycle oscillators (black dots) and of the mean field (grey dots) for positive and for negative coupling strength. N=5, ω0=1.

A similar scenario, but with only two oscillators that are mutually coupled as Embedded Image 4.2 Embedded Image 4.3 was recently studied by Prasad et al.(2006, 2008). The simpler model of two delay-coupled phase oscillators, Embedded Image 4.4 Embedded Image 4.5 was investigated by Schuster & Wagner (1989), who reported that the phase shift between the oscillators, α, alternates between α=0 and α=π, i.e. the synchronized solutions are either in-phase or out-of-phase. For increasing coupling and fixed delay, abrupt jumps in the frequency of the synchronized oscillation occur, which are accompanied by abrupt jumps in the value of α (i.e. phase-flip bifurcations occur as the coupling increases).

We note that the coupling scheme considered here differs from that of Schuster and Wagner and Prasad et al., because our scheme also includes self-feedback time-delayed terms, i.e. for N=2 the equations we considered are of the form Embedded Image 4.6 Embedded Image 4.7

Prasad et al.(2006, 2008) reported a phase-flip bifurcation induced by variation of the time delay: at a critical value of the delay time, the oscillators switch, from being in-phase to out-of-phase. The phase-flip bifurcation, which was shown to occur in different dynamical systems, was always accompanied by an abrupt jump in the frequency of the synchronized oscillators.

We stress that the phase-shifted dynamics resulting from the phase-flip bifurcation is different from the well-known amplitude death and oscillation death phenomena, where the oscillations in each unit of the coupled system cease completely (Reddy et al. 1998; Ozden et al. 2004; Ullner et al. 2007). Here, the amplitude of the individual limit cycle oscillators, displayed with dots in figure 6, remains unperturbed, i.e. nearly constant as the delay is varied. In that sense, one can expect that an even simpler model, composed of globally delay-coupled phase oscillators, can be sufficient to qualitatively explain the behaviour of the population of coupled neurons. A detailed comparison with the behaviour predicted by the time-delayed Kuramoto model (Niebur et al. 1991; Zanette 2000; Jeong et al. 2002), as well as with the behaviour of weakly optically coupled lasers (globally coupled through an external mirror; García-Ojalvo et al. 1999; Kozireff et al. 2000; Vladimirov et al. 2003), is in progress and will be reported elsewhere.

5. Summary and conclusions

Summarizing, we have studied the dynamics of an ensemble of neurons that are coupled through their time-delayed collective mean field. The individual neuron was modelled using a Hodgkin–Huxley-type conductance model with parameters chosen such that the uncoupled neuron displays autonomous subthreshold oscillations of the membrane potential. We stress that important simplifications of the model are that the coupling term is linear, the neurons are identical, noise is neglected and all time delays and coupling strengths between interacting neurons are identical. We found that the network generates a rich variety of oscillatory and firing activities that are mainly controlled by two time scales: the natural period of oscillation at the single neuron level and the delay time of the global coupling. When the neuronal oscillations are synchronized, they can be either in-phase or out-of-phase. The phase-shifted activity was interpreted as the result of a phase-flip bifurcation, also occurring in a set of globally delay-coupled limit cycle oscillators. At the bifurcation point, there is a transition from in-phase to out-of-phase (or vice versa) synchronized oscillations, which is accompanied by an abrupt change in the common oscillation frequency. This phase-flip bifurcation was recently investigated in two mutually delay-coupled oscillators (Prasad et al.2006, 2008) and can play a role in the mechanisms by which the neurons switch among different activity patterns. The study of the influence of network heterogeneities and noise is the subject of future work.

Acknowledgements

This research was supported by ‘Ramon y Cajal’ programme (Spain), the Ministerio de Ciencia e Innovación (Spain, project FIS2006-11452) and the GABA project (European Commission, FP6-NEST 043309).

Footnotes

References

View Abstract