## 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 2004*a*,*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 *i*th neuron, *V*_{i}, with *i*=1,…,*N* is that of the model by Braun *et al.*(1998, 2000), extended to incorporate delayed mean field coupling:
2.1
where *C*_{M} is the membrane capacitance, *T* the temperature, *I*_{l,i}=*g*_{ l}(*V*_{i}−*V*_{l}) a passive leak current, *I*_{Na,i} and *I*_{K,i} the fast sodium and potassium currents and *I*_{sd,i} and *I*_{sr,i} additional slow currents. These four currents are given by
2.2
where *k* stands for Na, K, sd and sr, *g*_{k} are the maximum conductances of the corresponding channels, and the activation variables *a*_{k} are given by
2.3
for K and sd,
2.4
and . The steady-state activations are defined by
2.5
The temperature-dependent factors *ρ*(*T*) and *ϕ*(*T*) are
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 2004*a*). 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 2004*b*), 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, *g*_{ l} and *V*_{l}, such as to write the coupling term in the form of *V*_{i}(*t*)−*V*_{j}(*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): *C*_{M}=1 *μ*F cm^{−2}, *g*_{ l}=0.1 mS cm^{−2}, *g*_{Na}=1.5 mS cm^{−2}, *g*_{K}=2 mS cm^{−2}, *g*_{sd}=0.25 mS cm^{−2}, *g*_{sr}=0.4 mS cm^{−2}, *V*_{l}=−60 mV, *V*_{Na}=50 mV, *V*_{K}=−90 mV, *V*_{sd}=50 mV, *V*_{sr}=−90 mV, *s*_{Na}=0.25 mV^{−1}, *s*_{K}=0.25 mV^{−1}, *s*_{sd}=0.09 mV^{−1}, *V*_{0Na}=−25 mV, *V*_{0K}=−25 mV, *V*_{0sd}=−40 mV, *τ*_{K}=2 ms, *τ*_{sd}=10 ms, *τ*_{sr}=20 ms, θ=0.17, *μ*=0.012, *A*_{1}=1.3, *A*_{2}=3.0, *T*=35^{o}C and *T*_{c}=25^{o}C. For these parameters, the uncoupled neuron displays subthreshold oscillations of period *T*_{0}∼130 ms. To integrate equation (2.1), it is necessary to specify the initial value of the membrane potential, *V*_{i}, 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,
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 and the values of and are calculated in a time window of duration *Δ**T* (typically, *Δ**T*=200*T*_{0}) 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.

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: , with and calculated in a time window after transients die away. The frequency of the oscillations of the *i*th neuron is defined as that of the main peak in the Fourier spectrum of the membrane potential *V*_{i}(*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 *τ*≤*T*_{0} 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).

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 4*c*); for smaller *N*, figure 4*a*,*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.

For delay values such that the neuronal oscillations are in-phase, *V*_{i}(*t*)=*V*_{j}(*t*)=*V* (*t*) ∀*i*, *j*, and the dynamics is as that of a single neuron with a recurrent delayed connection,
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:
4.1
where *Z*_{k} 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: (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).

A similar scenario, but with only two oscillators that are mutually coupled as
4.2
4.3
was recently studied by Prasad *et al.*(2006, 2008). The simpler model of two delay-coupled phase oscillators,
4.4
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
4.6
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

One contribution of 14 to a Theme Issue ‘Topics on non-equilibrium statistical mechanics and nonlinear physics’.

- © 2009 The Royal Society