# When long-range zero-lag synchronization is feasible in cortical networks

^{1}Donders Institute for Brain, Cognition and Behavior, Radboud University Nijmegen (Medical Centre), Nijmegen, Netherlands^{2}School of Psychology, Centre for Computational Neuroscience and Cognitive Robotics, University of Birmingham, Birmingham, UK

Many studies have reported long-range synchronization of neuronal activity between brain areas, in particular in the beta and gamma bands with frequencies in the range of 14–30 and 40–80 Hz, respectively. Several studies have reported synchrony with zero phase lag, which is remarkable considering the synaptic and conduction delays inherent in the connections between distant brain areas. This result has led to many speculations about the possible functional role of zero-lag synchrony, such as for neuronal communication, attention, memory, and feature binding. However, recent studies using recordings of single-unit activity and local field potentials report that neuronal synchronization may occur with non-zero phase lags. This raises the questions whether zero-lag synchrony can occur in the brain and, if so, under which conditions. We used analytical methods and computer simulations to investigate which connectivity between neuronal populations allows or prohibits zero-lag synchrony. We did so for a model where two oscillators interact via a relay oscillator. Analytical results and computer simulations were obtained for both type I Mirollo–Strogatz neurons and type II Hodgkin–Huxley neurons. We have investigated the dynamics of the model for various types of synaptic coupling and importantly considered the potential impact of Spike-Timing Dependent Plasticity (STDP) and its learning window. We confirm previous results that zero-lag synchrony can be achieved in this configuration. This is much easier to achieve with Hodgkin–Huxley neurons, which have a biphasic phase response curve, than for type I neurons. STDP facilitates zero-lag synchrony as it adjusts the synaptic strengths such that zero-lag synchrony is feasible for a much larger range of parameters than without STDP.

## Introduction

Coupling between different oscillators and pacemakers can generate a large range of different behaviors and has been a topic of study in many different conditions, for example in cardiac pacemaking and chemical oscillations (see e.g., Goldbeter, 1996; Koch and Segev, 1998; Roxin et al., 2005). A special case is the interaction between neurons, which gives rise to neuronal oscillations in particular frequency bands. Neuronal gamma band synchronization has been reported in many species and in a large number of brain structures for a variety of sensory and motor tasks (Gray et al., 1989; Fries et al., 2001, 2008; Pesaran et al., 2002; Schoffelen et al., 2005). This large-scale synchronization of multiple cortical areas has been postulated as a potential mechanism for integration and coordination of neuronal activity in cognitive tasks (Engel et al., 1992; Singer and Gray, 1995; Fries, 2005).

The first studies on this topic presented experimental evidence that the relative phase of gamma oscillations in widely separated brain areas is near zero (Frien et al., 1994; Roelfsema et al., 1997; Castelo-Branco et al., 1998; Rodriguez et al., 1999; Gross et al., 2004). This result was remarkable since synchronization requires interactions between distant brain areas, which come with considerable delays due to axonal conduction and synaptic transmission. Many further studies hence investigated how distant oscillatory brain regions can synchronize at zero-lag in spite of non-negligible delays. Several theoretical studies have argued that direct mutual pulse-coupling between two dynamical systems with delays and excitatory synapses cannot easily lead to zero-lag synchrony (Ernst et al., 1995, 1998; Goel and Ermentrout, 2002; Zeitler et al., 2009). Therefore, Fischer et al. (2006) and Vicente et al. (2008) suggested that zero-lag synchronization between brain areas might be mediated by a third (relay) oscillator. A potential candidate for this neuronal relay oscillator is the thalamus (Gollo et al., 2010; Theyel et al., 2010).

In order to obtain a better understanding of the possibilities for zero-lag synchronization of distant brain areas, we have investigated the proposed network of neuronal oscillators coupled indirectly by a relay oscillator (Fischer et al., 2006; Vicente et al., 2008). We have investigated the model with both type I Mirollo–Strogatz neurons, as well as with type II Hodgkin–Huxley neurons. We focused in particular on the robustness of zero phase synchronization as a function of both the delay time between oscillators and the strength of synaptic coupling. Various types of synaptic coupling were investigated, including Spike-Timing Dependent Plasticity (STDP), since STDP has been suggested to contribute to efficient information transmission (Buonomano and Maass, 2009; Lindner et al., 2009; Hennequin et al., 2010). Our results show that zero-lag synchrony can occur, especially for models with Hodgkin–Huxley type II neurons. STDP facilitates zero-lag synchrony as STDP modifies synaptic strengths and thereby allows a larger range of initial synaptic strengths that may lead at zero-lag synchronization.

## Materials and Methods

The model consists of three coupled neuronal oscillators (see Figure 1A). Each oscillator can be considered as a single neuron or as a population of neurons, where the activity of the neurons within each population is assumed to be homogeneous and highly synchronized. As illustrated in Figure 1A, the “relay” or “inner” oscillator (oscillator 2) is coupled bi-directionally with two “outer” oscillators (oscillators 1 and 3). The outer oscillators are only coupled with the relay oscillator but not directly to each other. As starting point we assume that all three oscillators are identical with the same intrinsic firing period *T*_{0}.

**Figure 1. Schematic of the model and its 1:1 phase-locked zero-lag synchrony mode**. **(A)** Sketch of oscillators 1 and 3 with bi-directional pulse-coupling via relay oscillator 2. ε* _{ij}* represents the excitatory synaptic weight from oscillator

*j*to oscillator

*i*, τ

*the conduction delay between relay oscillator 2 and oscillator*

_{k}*k*.

**(B)**An example of 1:1 phase-locked zero-lag synchrony mode. Vertical bars represent the spike times of the oscillators.

*T*

_{0}is the intrinsic period of the oscillators,

*T*<

*T*

_{0}the common period of the oscillators in the phase-locked mode, and θ

*T*the delay between the firing of the relay and the outer oscillators. ${t}_{1,n}^{\text{spike}},$${t}_{2,n}^{\text{spike}},$ and ${t}_{3,n}^{\text{spike}}$ are the

*n*-th spike times of oscillators 1, 2, and 3, respectively.

We have used the Mirollo–Strogatz (MS) phase oscillator (Mirollo and Strogatz, 1990) and the classical Hodgkin and Huxley (HH) neuron (Hodgkin and Huxley, 1952) with parameters given in Vicente et al. (2008). For the MS phase oscillator, *T*_{0} is chosen to be 25 ms and *T*_{0} of the HH neuron is 14.66 ms with the parameters given in Vicente et al. (2008). In Figure 1A ε* _{ij}* represents the coupling strength from pre-synaptic oscillator

*j*to postsynaptic oscillator

*i*. All synaptic couplings in the model are excitatory. The delay time τ

*represents the conduction time for spikes along axons that connect oscillator*

_{k}*k*with the relay oscillator. Here we make the simplifying assumption that delay times are constant and symmetric (i.e., the delay time from oscillator

*i*to

*j*is equal to that from oscillator

*j*to

*i*). These delay times are typically considerably shorter than the period of neuronal oscillations (Fries, 2005). Therefore in the following we will only consider conduction delays shorter than

*T*

_{0}/2, with τ

*expressed as a fraction of the intrinsic period*

_{k}*T*

_{0}and thus in the range between 0 and 0.5.

The MS phase oscillator (Mirollo and Strogatz, 1990) is characterized by a voltage-like state variable *f* ∈ [0, 1], which increases monotonically from 0 toward the threshold value *f* = 1. Within a cycle, the state of the uncoupled neuronal oscillator is defined by a monotonically increasing concave function *f *(φ):[0, 1] → [0, 1]:

with a phase variable φ ∈ [0, 1] and a dissipation parameter *b*. *T*_{0} is the intrinsic firing period of the oscillator. When the threshold is reached, the oscillator fires, the state variable *f* is reset to zero, and the cycle repeats. As in Ernst et al. (1995) and Zeitler et al. (2009), the setting *b* = 3 is used throughout this study. The MS neuron is a so-called type I neuron (Izhikevich, 2007), where excitatory input always gives a phase advance of the neuronal oscillator.

For the classical HH neuron the membrane potential *V* is governed by the differential equation:

with the membrane capacitance *C* and the maximal conductances of sodium *g*_{Na}, potassium *g*_{K}, and leakage *g*_{L}. The corresponding reversal potentials (*E*_{Na}, *E*_{K}, *E*_{L}) and the external current *I*_{ext} are as given in Vicente et al. (2008). The voltage-gated ion channels m, h, and n are described by first order differential equations. Note that the expression for α* _{n}* in Vicente et al. (2008) was not correct and should read α

*(*

_{n}*V*) = [(

*V*+ 55)/100]/{1 − exp[ − 0.1(

*V*+ 55)]}. The excitatory synaptic current

*I*

_{syn}is −ε

*(*

_{ij}S*t*)

*V*, where ε

*is the maximum synaptic conductance and*

_{ij}*S*(

*t*) is the Dirac delta function (in case of an “instantaneous synapse”) or an alpha function (in case of an “alpha synapse”). This classical HH neuron is a so-called type II neuron (Izhikevich, 2007), where excitatory input in early phases of the firing cycle causes a phase delay but a phase advance in later phases of the firing cycle.

### Synaptic Coupling

We investigate two models for pulse-coupling between the oscillators. For an instantaneous synapse with coupling strength ε* _{ij}* the neuronal state

*f*of the MS phase oscillator after arrival of a spike from oscillator

_{i}*j*is incremented instantaneously

with φ* _{i}* the phase of the postsynaptic oscillator

*i*at the time of the spike arrival. The phase φ

*for which the oscillator reaches the threshold after spike input [i.e.,*

_{i}*f*(φ

_{i}*) + ε*

_{i}*= 1] is called the critical phase φ*

_{ij}*(Zeitler et al., 2009) and is given by*

_{c}If a spike arrives at φ* _{i}* < φ

*,*

_{c}*f*will increase instantaneously by ε

_{i}*. The instantaneous change of the state*

_{ij}*f*by ε

_{i}*corresponds to a phase shift $\Delta {\phi}_{i}={{f}_{i}}^{-1}\left[{f}_{i}\left({\phi}_{i}\right)+{\epsilon}_{ij}\right]-{\phi}_{i},$ which yields*

_{ij}with

where φ* _{i}* is the phase of the postsynaptic oscillator just before arrival of the input spike. Note that φ

*(ε*

_{c}*) = [1 − χ*

_{ij}*(ε*

_{b}*)]/[1 + β*

_{ij}*(ε*

_{b}*)].*

_{ij}Tsubo et al. (2007) measured the phase shifts of layer-5 and layer-2/3 pyramidal neurons in rat motor cortex. The maxima of the averaged phase shifts were found to be in the second half of the oscillatory period for these neurons at all frequencies (including the gamma band; see their Figure 4). Hence we require here that the maximum of the phase shift Δφ* _{i}* is in the second part of the intrinsic cycle of the oscillator, and therefore φ

*> 0.5. This restriction imposes an upper bound of ε*

_{c}*< 1 − ln[(1 +*

_{ij}*e*)/2]/

^{b}*b*= 0.21 for

*b*= 3 on the synaptic strength through Eq. 5. In case of ε

*= 0.21, a spike can cause an increase in the state variable*

_{ij}*f*of about 21% of the difference between rest state

*f*= 0 and the threshold value

*f*= 1. For the HH neuron, the upper bound for the maximum synaptic conductance is 3.15 mS/cm

^{2}, which corresponds to an increase of 21% of (

*V*

_{onset}−

*V*

_{rest}), since the resting potential

*V*

_{rest}and the onset of the action potential

*V*

_{onset}are near −65 and −50 mV, respectively.

A more realistic synaptic coupling model is provided by the so-called alpha function. For an “alpha synapse,” the postsynaptic potential after arrival of a spike at time *t*_{0} at the synapse with strength ε* _{ij}* is described by

where τ_{syn} > 0 is the synaptic rise time of the input. Unless stated otherwise, τ_{syn} = 2 ms in this study. The numerical simulations are implemented using an Euler scheme with a time step size equal to 2.5 μs (≈10^{−4}*T*_{0}). In our implementation, there is no current reset after spiking of the postsynaptic neuron, i.e., the “tail” of the alpha function is carried over into the next cycle.

### Spike-Timing Dependent Plasticity

In general, the synaptic coupling strength is not constant, but varies depending on pre- and post-synaptic activity due to STDP (Hebb, 1949; Bi and Poo, 1998). We have implemented the additive STDP rule (Froemke et al., 2006) for both instantaneous and alpha synapses. For a pre-synaptic spike at arrival time ${t}_{k}^{\text{arr}}$ and a postsynaptic spike at ${t}_{l}^{\text{spike}}$ the fractional synaptic modification *W*(Δ*t*) is given by

with $\Delta t\equiv {t}_{l}^{\text{spike}}-{t}_{k}^{\text{arr}}.$

The spike arrival time ${t}_{k}^{\text{arr}}$ is defined as the time of the onset of the postsynaptic potential, as in the experimental protocol by Bi and Poo (1998). Unless stated otherwise, the time constants for potentiation τ_{+} and depression τ_{−} are τ_{+,0} ≡ 16.8 ms and τ_{−,0} ≡ 33.7 ms, respectively. In this study the standard values for the maximum amplitude of potentiation *A*_{+} and depression *A*_{−} are *A*_{+,0} ≡ 0.78 and *A*_{−,0} ≡ −0.27, respectively. These standard parameter values for STDP were fitted to the data from Bi and Poo (1998), who determined the fractional synaptic modification *W* as the relative change in synaptic strength after evoking an input and output spike pair 60 times. We assume that the change of the synaptic weight Δε caused by each input spike is constant and defined by Δε = ε_{0} *W*(Δ*t*)/60 (Song et al., 2000; van Rossum et al., 2000) with initial synaptic weight ε_{0} and *W*(Δ*t*) as in Eq. 9. Therefore, the learning rule used here for a particular pair of pre- and post-synaptic spikes is given by ε_{0} → ε_{0} + ε_{0} *W*(Δ*t*)/60, in agreement with Lee et al. (2009).

The simulation procedure of the network for STDP is as follows: The evolution of the state of the network is studied over 60 consecutive sessions. At the start of each session, the initial phases of the three oscillators are chosen arbitrarily from a uniform random distribution. The first session starts with equal coupling strengths for all synapses, which then change due to STDP. All succeeding sessions start with the coupling strengths that resulted at the end of the previous session, but with re-randomized phases of the oscillators to prevent that the system converges into a locally rather than a globally stable state.

In order to be physiologically relevant for synchronization, convergence should not take too much time. Therefore, we assume that convergence to a stable relative phase between oscillators 1 and 3 in the model takes place within a session consisting of *n*_{sess} = 15 cycles. Thus with STDP we typically run 60 sessions of 15 cycles each to investigate convergence of the network into a stable synchronization state. Without STDP the couplings do not change and a single session of 15 cycles is used. To avoid any spurious dependence on the random initial phases, we repeat each calculation 35^{3} = 45,875 times. Thus in the STDP case almost 40 million cycles are computed for every chosen setting of coupling strengths and delay times.

### Phase-Locking Equation and Stability Analysis

Phase-locking equations are useful to determine the new period of a network and the relative phases of coupled oscillators (van Vreeswijk et al., 1994; Bressloff and Coombes, 1998). For the simple MS phase oscillator, the phase-locking equation can be derived analytically. A full analysis for our model is provided in the Appendix, where we derive the relationships between synaptic weight, conduction delay, and the new period of the oscillators for zero-lag 1:1 phase-locked firing of the outer oscillators. Here we will just consider a simple example (see Figure 1B).

Without loss of generality, we set the time of the (*n*−1)-th firing of the relay neuron ${t}_{2,n-1}^{\text{spike}}=0$ and call the period of the 1:1 phase-locked oscillators *T*. For zero-lag synchrony there is a þeta ∈ [0, 1) such that the outer oscillators fire at ${t}_{1,n-1}^{\text{spike}}={t}_{3,n-1}^{\text{spike}}=\theta T,{t}_{1,n}^{\text{spike}}={t}_{3,n}^{\text{spike}}=\theta T+T,$ etc. For equal delays τ_{1} = τ_{3} = τ, both spikes from the outer oscillators arrive at the relay oscillator at the phase τ + θ*T*/*T*_{0} (if τ + θ*T*/*T*_{0} ≤ *T*/*T*_{0}) with excitatory coupling strength ε. The two synaptic inputs reduce the period of the relay oscillator from the intrinsic period *T*_{0} to

with Δφ as defined in Eq. 6. Since the relay neuron fires at *T*, the outer oscillators will each receive a spike at *T* + τ*T*_{0}. The outer neurons spiked at θ*T* and hence their phase will be τ + (1 − θ)*T*/*T*_{0} (if τ + (1 − θ)*T*/*T*_{0} ≤ *T*/*T*_{0}) at the arrival of the spike from the relay neuron. Therefore the period of the outer oscillators is given by

Equations 10 and 11 are the 1:1 phase-locking equations for the mentioned conditions.

Assume now that both arrival phases exceed the critical phases φ* _{c}*(2ε) and φ

*(ε), respectively. Then rewriting Eqs 10 and 11 using Eq. 6 yields*

_{c}Thus the period of all oscillators is twice the conduction delay and activity switches between inner and outer oscillators every half period. All oscillators immediately spike upon spike input. We will call this mode the “driven synchrony” (DS) mode (see Figure 2B). Two other typical modes are also illustrated in Figure 2: “slave synchrony,” where only the relay oscillator spikes directly after input, and “pacemaker synchrony,” where only the outer oscillators do so.

**Figure 2. Features of synchronization depending on delay time and coupling strength**. **(A)** Synchronization quality (SQ) and **(C)** convergence promptness (CP) for model with instantaneous synapses. Thick black lines in **(A,C)** indicate boundaries of dynamics calculated analytically using the phase-locking equation. **(B)** Illustration of synchrony modes dominant in regions indicated by roman numerals in **(A,C)**. SQ is high in region I, low in regions II and IV, and intermediate in region III. Slave synchrony (regions II and IV) is not stable, whereas driven (region I) and pacemaker (region III) synchrony are asymptotically stable. CP is highest in region I, indicating synchronization within about four cycles.

To investigate the stability of DS we assume small perturbations. Since we use the relay neuron spikes as reference time, the perturbation affects the phase of the outer oscillators to (τ + δφ_{1}, τ + δφ_{3}) at ${t}_{2,n-1}^{\text{spike}}.$ DS is asymptotically stable, if there is a δ > 0 such that the phases of the outer oscillators 1 and 3 will be closer to τ at the next spike at ${t}_{2,n}^{\text{spike}}$ for $\sqrt{\delta {\phi}_{1}^{2}+\delta {\phi}_{3}^{2}}<\delta $. Since φ* _{c}*(ε) < 2τ, we can define a value δ ≡ 2τ − φ

*(ε) > 0, and thus φ*

_{c}*(ε) = 2τ − δ. The spike from the relay arrives when the phase of oscillator 1 is in the range 2τ − δ < φ*

_{c}_{1}< 2τ + δ, which exceeds this critical phase φ

*(ε). Therefore oscillator 1 will spike immediately after receiving input from the relay neuron, and the same is true for oscillator 3. Thus perfect synchronization is re-established as long as the perturbation was sufficiently small, proving that the DS mode is an asymptotically stable mode. To determine the stability of other synchrony modes it is necessary to calculate the eigenvalues λ of the Jacobian of the (phase) return map (Zeitler et al., 2009), see the Appendix for details.*

_{c}We now briefly consider the effect of STDP on slave synchrony. The coupling strengths from the outer oscillators to the relay oscillator remain unchanged for slave synchrony, since the relay oscillator immediately spikes upon input from the outer neurons. In the Appendix we show that the coupling strengths from the relay oscillator to the outer oscillators increase with ε* _{n}* > ε

*for τ ≥ 0.25 using the STDP window*

_{n - 1}*W*(Δ

*t*) of Eq. 9, since the spike from the relay oscillator arrives well before the outer neurons spike. Hence synaptic weights increase to ${\phi}_{c}^{-1}\left(2\tau \right)=1-\text{ln}\left[2\tau \left({e}^{b}-1\right)+1\right]/b$, at which point the outer oscillators also spike immediately after input and STDP stops. Thus STDP turns slave synchrony into DS.

In the driven, pacemaker, and slave synchrony modes oscillators 1 and 3 spike simultaneously. There are other, asynchronous stable modes where this is not the case. Figures 3C,D show the dynamics of the oscillators for a typical asynchronous case. This mode will be referred to as a “pacemaker-slave” because the relay drives oscillator 1 and is driven by oscillator 3.

**Figure 3. Dominance of non-zero relative phases**. **(A)** Normalized histogram of the relative phase φ* _{r}* for stable states with τ = 0.25 and ε = 0.1. For about 90% of possible initial phases, the system converges to a non-zero relative phase.

**(B)**|φ

*| As function of the conduction delay τ and the synaptic weight ε, with averages projected into side panels. |φ*

_{r}*| decreases when ε increases, and has a maximum for τ near 0.2.*

_{r}**(C,D)**Dynamics of “pacemaker-slave” synchrony corresponding to non-zero relative phase. Solid lines represent spike times for parameter values τ = 0.25 and ε = 0.1, dotted lines for a stronger synaptic weight ε = 0.14 in

**(C)**and longer delay time τ = 0.35 in

**(D)**, respectively.

### Synchronization Measures

For a proper analysis of synchronized firing, a quantitative definition of synchrony is required. In experiments perfect synchrony will never be observed as noise causes small variations in the timing of action potentials. Instead we will define synchrony as firing of two oscillators within a small time window, where the time window should be sufficiently large to eliminate the effect of noise, and sufficiently small to provide an accurate estimate of synchrony. While we could measure synchrony here to the limits of the numerical accuracy of our simulations, we define synchrony as the condition when the spikes of the outer oscillators occur within $\left|\phantom{\rule{2.77695pt}{0ex}}{t}_{1,n}^{\text{spike}}-{t}_{3,n}^{\text{spike}}\right|\le 0.02{T}_{\text{0}}$ (Engel and Singer, 2001), which is 0.5 ms for simulations of the model with MS phase oscillators and 0.2932 ms for the model with the classical HH neurons.

Synchronized firing of the outer oscillators does not only depend on the synaptic weights and delay times, but also on the initial phase of the three oscillators. Hence we have to consider the robustness of convergence to synchronous firing for variations in the initial relative phase of the neuronal oscillators. Therefore, we define “*synchronization quality*” (SQ) as the fraction of the number of initial phase combinations that leads to stable synchronous firing of the outer oscillators. SQ has a value between 0 and 100%, where 100% means that the outer oscillators will always converge to synchronous spiking within the simulation period, independent of their initial phases. This provides a measure of the attraction domain of the initial phases for reaching synchronization of the outer oscillators. To determine SQ, we repeat our simulations 35^{3} = 42,875 times with randomly chosen initial phases.

In order to investigate the impact of the various parameters of the STDP learning window, we wish to derive a single value for the SQ, rather than one value for each combination of synaptic weights and conduction delay. For this purpose we consider various synaptic weights in 100 evenly spaced steps in the range from 0.01 to 0.21 and conduction delays in 100 evenly spaced steps in the range from 0.01 to 0.49. Then we average the SQ over these 10,000 pairs to obtain an “average SQ” for each parameter set of the learning window. Note that for an STDP run with 60 sessions, this means that we compute almost 400 billion cycles of the model for each change of the STDP parameters.

For some combinations of the synaptic weights, delay times, and initial phases the state of zero-lag synchronization is reached faster than for other combinations. Therefore, we define a “*convergence promptness*” (CP) for the network to achieve zero-lag synchronization of the outer oscillators. This is calculated as CP = SQ · (1 − 〈*n*_{sync}〉/*n*_{sess}), where 〈*n*_{sync}〉 is the average number of intrinsic periods *T*_{0} needed to achieve zero-lag synchrony. When there is no zero-lag synchrony established within *n*_{sess} = 15 cycles, CP equals to 0. Note that 〈*n*_{sync}〉 ≤ *n*_{sess} with *n*_{sess} = 15 for the simulations in this study, and that the measure is multiplied with the SQ to account for the readiness of the system to achieve zero-lag synchrony at all.

Finally, not only are we interested in synchrony of the outer oscillators, but also in the relative phase φ* _{r}* for stable modes in which the oscillators 1 and 3 are not in synchrony. The relative phase between the outer oscillators is defined by φ

*≡ $\left({t}_{3,n}^{\text{spike}}-{t}_{1,n}^{\text{spike}}\right)/{T}_{0}$ with a value between −1 and 1. Since values for φ*

_{r}*of −1, 0, and 1 all refer to the same state, we report φ*

_{r}*rescaled to the range −1/2 to 1/2 by periodically mapping [−1, −0.5] to [0, 0.5] and [0.5, 1] to [−0.5, 0], respectively.*

_{r}## Results

### Equal Delay Times and Coupling Strengths

Figures 2A,C show simulation results for the SQ and the CP, respectively, as a function of the conduction delay τ and synaptic weight ε for the MS model with instantaneous synapses, identical delay times (τ_{1} = τ_{3} = τ) and identical coupling strengths (ε_{12} = ε_{21} = ε_{23} = ε_{32} = ε). The three solid black lines show the analytically calculated boundaries between different modes of synchrony using the phase-locking equations, for details see Section “Synchrony Modes Derived from the Phase-Locking Equation” in the Appendix. The upper right area, indicated by IV, is bounded by the line

where χ* _{b}* is defined in Eq. 7. For a detailed derivation see the discussion below Eq. A12 in Appendix. Our simulations show that in region IV zero-lag synchronization occurs mainly when the relay oscillator spikes immediately upon arrival of synaptic input from the outer neurons, but the outer neurons do not after input from the relay neuron, see the middle panel of Figure 2B. The period

*T*

_{IV}of this “slave synchrony” mode is given by Eq. A12 “SS2” in Appendix.

The line separating regions I and II is given by τ = φ* _{c}*(ε)/2. In region I, zero-lag synchronization occurs mainly when both the relay and outer oscillators spike immediately when synaptic input arrives, see the upper panel in Figure 2B. This represents “DS” with period

*T*

_{I}for all neurons, see Eq. A7 in Appendix. The line separating regions II and III is given by

where φ* _{c}*, χ

*, and β*

_{b}*defined in Eqs 5 and 7, respectively. For a detailed derivation see the discussion below Eq. A10 in Appendix. Region II shows mainly slave synchrony like region IV, but with a different period*

_{b}*T*

_{II}, see Eq. A12 “SS1” in Appendix. In region III, zero-lag synchronization occurs mainly when the outer oscillators spike immediately upon arrival of synaptic input while the relay oscillator does not. This “pacemaker synchrony,” cf. the lower panel in Figure 2B, has period

*T*

_{III}according to Eq. A10 “PS1” in Appendix.

Figure 2A shows that the SQ changes from high (region I) to low (region II) and back to intermediate values (region III) when the conduction delay decreases from 0.5 to 0. Stability analysis indicates that regions II and IV show poor SQ because slave synchrony, which dominates in these regions, is unstable, whereas driven and pacemaker synchrony, which dominate in region I and III, are asymptotically stable and robust against changes in the initial phases. Most combinations of delay and synaptic weight yield a SQ that is below 50%. Figure 2C shows that the delay and the synaptic weight in region I yield fast convergence to synchronization within about four cycles. In the other regions, zero-lag synchrony is established much more slowly or not at all.

Figure 3A shows a histogram of the relative phase φ* _{r}* for the stable pacemaker-slave (φ

*≠ 0) and unstable slave synchrony (φ*

_{r}*= 0) modes of the model with a conduction delay τ = 0.25 and a synaptic weight ε = 0.1 (parameters in region II, cf. Figure 2A). Two non-zero relative phases, corresponding to the stable pacemaker-slave mode, and one zero-lag phase, corresponding to the unstable slave synchrony mode, occur. The non-zero relative phases have the same absolute value because of the symmetry of the network. For about 90% of the initial phases of the oscillators the system converges to the two non-zero relative phases and for 10% to a state with zero-lag synchrony. For other values of the delay and the weight in region II, III, and IV, the histograms are qualitatively similar to that shown in Figure 3A, i.e., two non-zero and one zero-lag relative phase. Thus in general zero-lag synchrony only occurs for a very limited set of initial phases in parameter regions II, III, and IV.*

_{r}To illustrate how the non-zero relative phase varies for different parameter values in the network, Figure 3B shows the absolute value |φ* _{r}*| as a function of the conduction delay τ and the synaptic weight ε. The three black lines in Figure 3B indicate the boundaries between the regions I, II, III, and IV as in Figure 2 calculated analytically using the phase-locking equation. To understand the results in Figure 3B, we consider the “pacemaker-slave” mode that corresponds to the positive relative phase illustrated in Figures 3C,D. Vertical solid lines represent the spiking times for ε = 0.1 and τ = 0.25, and the dotted lines correspond to a stronger weight ε = 0.14 in Figure 3C and a longer delay τ = 0.35 in Figure 3D, respectively. When the synaptic weight increases in Figure 3C, oscillators 2 and 3 will spike sooner after input from oscillator 1, whereas oscillator 1 always spikes immediately after input from oscillator 2. Hence the relative phase between oscillators 3 and 1 decreases when the coupling strength increases. At longer delay times in Figure 3D, the synaptic inputs from oscillator 2 arrive later in time at oscillators 1 and 3 as indicated by the dotted arrows. The input causes oscillator 1 to spike immediately, but not so for oscillator 3. For the longer delay time, oscillator 3 will be further in its natural cycle and therefore it will spike sooner after the input. Thus oscillators 1 and 3 will both spike later, but their relative phase difference is reduced for longer delays. The left panel of Figure 3B shows the relative phase as a function of the delay, averaged over all values of synaptic coupling. The lower panel shows the relative phase as a function of coupling strength, averaged over all delays.

Next, we study synchronization for a model with HH neurons. Figures 4A,E show simulation results for the model with HH neurons for the SQ and the CP, respectively, as a function of the conduction delay τ and synaptic weight ε for instantaneous synapses. Comparing Figures 4A,E with the same results for the MS oscillator in Figures 2A,C shows that the model with HH neurons yields a high SQ and large CP for a larger range of τ and ε values than the MS neuron.

**Figure 4. Excellent synchronization for alpha synapses with short rise times**. (Top row) Synchronization quality (SQ) and (bottom row) convergence promptness (CP) for a Hodgkin–Huxley based model with instantaneous synapses **(A,E)** and alpha synapses for various rise times **(B–D,F–H)**. Comparing **(A,E)** with the results for the MS neurons (see Figures 2A,C) shows that SQ and CP clearly improve for HH neurons. Moreover, SQ and CP for short rise times of the alpha synapses **(B,F)** are better than for instantaneous synapses **(A,E)**, but decrease for longer rise times.

Figure 4A also shows that perfect SQ is not possible when the delay time is shorter than approximately 0.05 of the intrinsic period, which roughly corresponds to half the refractory period of the HH neuron. If the time delay is 0.05*T*_{0} or less, the time interval from spiking of the outer neurons and spike input from the relay neurons to the outer neurons is less or equal to the refractory period. In that case, input from the relay neuron to the outer neurons arrives in the refractory period, which effectively reduces the coupling strength when the outer neurons tend to synchronize and thereby disables the zero-lag DS mode when the delay time is short.

Figures 4B–D,F–H show the SQ and CP, respectively, as a function of the conduction delay τ and the synaptic weight ε for the HH neuron model with alpha synapses of various synaptic rise times (left, middle, and right columns for τ_{syn} = 1, 2, and 3 ms, respectively). Quite surprisingly, in the context of the results in Figures 4A,E, the SQ is high for very small time delays for τ_{syn} = 1 ms (Figure 4B). This can be understood from the fact that the alpha synapse adds an effective delay such that spike input from the relay neuron to the outer neurons arrives after the refractory period when the outer neurons fire in perfect synchrony. When the synaptic rise time increases, the range of time delays and synaptic strengths with high SQ decreases, cf. Figures 4B–D. To understand this, assume that synaptic inputs from the relay oscillator arrive when the phase of oscillators 1 and 3 is φ_{1} and φ_{1} + Δ, respectively. If Δ differs from zero, then oscillators 1 and 3 will spike at different times, unless the input is strong enough to immediately initiate a spike in both oscillators. For larger rise times and the same synaptic strength ε, less input per unit of time is received, since the synaptic input is spread out over a longer time. Therefore, input from an instantaneous synapse (equivalent to “τ_{syn} = 0”) or fast alpha synapse (small τ_{syn}) can more readily synchronize the outer oscillators than a slow alpha synapse (large τ_{syn}).

High SQ is possible for short delay times (τ < 0.05), if the synaptic rise time is short (Figure 4B), but less so when the synaptic rise time increases (Figures 4C,D). This suggests that there is a range of synaptic rise times which favor a high SQ for short time delays. Figure 5 shows the SQ, as a function of the synaptic rise time τ_{syn}, evaluated at τ = 0.02 for ε = 1 (red), 2 (green), and 3 (blue) mS/cm^{2}. Perfect SQ is obtained when the synaptic rise time is moderately fast, i.e., approximately between 0.5 and 1 ms, for a large range of synaptic coupling strengths. Therefore, moderately fast synaptic rise times favor zero-lag synchrony.

**Figure 5. Optimal synchronization for moderately fast synaptic rise time**. Synchronization quality at τ = 0.02 for three synaptic strengths ε = 1 (red), 2 (green), and 3 (blue) mS/cm^{2} in a model with alpha synapses and Hodgkin–Huxley neurons, shown as a function of the synaptic rise time τ_{syn}. Almost perfect synchronization quality can be obtained when the synaptic rise time is moderately fast, i.e., approximately between 0.5 and 1 ms.

### STDP Facilitates Zero-Lag Synchronization

We now investigate the effect of STDP on synchronization of the MS neuron, starting with the simple instantaneous synapses. Figure 6 shows the SQ and CP with short (left column) and long STDP adaptation (right column), i.e., after the first and sixtieth session, respectively. To allow an easy comparison to the results without STDP in Figure 2, we have drawn the same thick black lines in Figure 6 which separate regions with different dynamics as in Figure 2. Note that the synaptic strength ε along the horizontal axes in Figure 6 represents the *initial* synaptic weight, not the final values after adaptation by STDP. The results in Figures 6A,B show that the effect of STDP is a gradual expansion of the range of coupling strengths which give rise to synchronization. With STDP, the weak coupling gradually increases to larger synaptic strengths that allow synchronous firing, corresponding to DS. Likewise, the CP increases, see Figures 6C,D.

**Figure 6. STDP improves long-range synchronization for instantaneous synapses and decreases the fraction of non-zero relative phases**. **(A,B)** Synchronization quality after the first and the 60-th session, respectively. **(C,D)** Corresponding convergence promptness. For conduction delays τ ≥ 0.25, synchrony improves due to an increase of initial weights ε by STDP. **(E)** |φ* _{r}*| as in Figure 3B but after 60 learning sessions.

**(F)**Synaptic weights from the outer oscillators to the relay oscillator;

**(G)**synaptic weights from the relay oscillator to the outer oscillators. Changes in |φ

*| relative to that in Figure 3B are largely caused by an increase of the synaptic weights by STDP. To compare with the results without STDP, we show the solid black lines which separate regions with different dynamics in Figure 2.*

_{r}The higher SQ and the faster CP in region II for large delay times in Figures 6A–D can be understood as follows: by increasing coupling strengths STDP converts slave synchrony, which is unstable and a dominant zero-lag synchronization mode in region II, into DS, which is asymptotically stable. For slave synchrony STDP increases the coupling strength from the relay oscillator to the outer oscillator. The sharp border is related to the maximum coupling strength ε_{max} = 0.21: the slave to drive synchrony conversion can only occur if τ exceeds ${\tau}_{c}^{+}\equiv {\phi}_{c}\left({\text{\epsilon}}_{\text{max}}\right)/2\approx 0.25.$ Notice that the fully improved domain of attraction for zero-lag synchrony is reached only after 60 sessions with STDP. This suggests that STDP can contribute to zero-lag synchrony, but generally only after many cycles of weight adaptation (here up to 900 oscillations).

Figure 7 shows the number of sessions for STDP required to obtain 100% SQ, if the network begins with four equally strong coupling strengths in the range between 0 < ε ≤ 0.21 and delay times ${\tau}_{c}^{+}\le \tau <{\phi}_{c}\left(\text{\epsilon}\right)/2,$ i.e., where STDP changes slave synchrony into DS as just discussed. The solid black lines in Figure 7 separate the regions with different dynamics as in Figure 2 (without STDP). Figures 7A,B are obtained from direct simulations and from analytical calculations, respectively. The analytical results are obtained using Eqs A26 and A27 in Appendix iteratively. The simulated and analytical results are in good agreement. The network requires a slightly smaller number of sessions to reach a high SQ value in the analytical calculations, because these start from the condition of slave synchrony, whereas the direct simulations start from random initial phases and reach the driven state after achieving slave synchrony first. For smaller initial coupling strengths it takes more time to converge to stable zero-lag synchronization with STDP: 20 cycles (about 0.5 s) or more.

**Figure 7. Time required to reach 100% synchronization quality by STDP**. **(A)** Number of sessions to achieve 100% synchronization quality according to direct simulations; **(B)** same according to analytical calculations for converting slave into driven synchrony for the relevant part ${\phi}_{c}^{+}\le \tau <{\phi}_{c}\left(\epsilon \right)\u22152$ of region II. Solid black lines in **(A,B)** separate regions with different dynamics as shown in Figure 2.

Figure 6E shows the absolute value of the non-zero relative phase |φ* _{r}*| as a function of the delay time τ and the synaptic weight ε after convergence of the network to a stable state by STDP. This result should be compared to the case without STDP in Figure 3B. Note that the synaptic strength ε along the horizontal axes in Figure 6 represents the initial synaptic weight at the start of the simulations, which is assumed to be equally strong for the four connections (cf. Figure 1A) and will change during the simulations by STDP. The average values of final synaptic weights are shown as a function of the delay and the weight in Figure 6F (connections from the outer to the relay oscillator) and Figure 6G (connections from the relay to the outer oscillators), respectively. After training for 60 sessions, the values of |φ

*| for the stable states of the network, shown in Figure 6E, significantly decrease in region II while |φ*

_{r}*| changes slightly in regions I, III, and IV as compared to the values before adaptation of the synapses by STDP (Figure 3B). This is largely caused by STDP-induced increase of the coupling strength in region II, as shown in Figures 6F,G. The final weights overall are quite similar in both directions in region II. In region III they differ because pacemaker synchrony is dominant, which implies different timing of pre- and post-synaptic firing for the relay and outer neurons, and therefore different effects of STDP. The final weights in regions I, II, and IV in Figures 6F,G are similar because, after many learning sessions, DS is dominant here. In DS, the firing behavior of the relay and outer oscillator is the same, i.e., a synaptic input arrives at phase 2τ and makes the neuron spike instantaneously. Since all connections start with the same initial weight, the same weight adaptation is applied to the connections from the relay to the outer oscillators and vice versa.*

_{r}In order to investigate the robustness of our results for variations of the STDP parameter values, we return to the instantaneous synapses. We have varied the amplitudes *A*_{+} and *A*_{−} of the fractional synaptic modification *W*, see Eq. 9 and Figure 8A. Figures 8B,C show the SQ and the CP, respectively, for various values of *A*_{+} and *A*_{−} as a function of the number of learning sessions. The black lines in Figure 8 show the results for the asymmetric learning window with our standard parameters of STDP obtained from the Bi and Poo (1998) data, while the red lines show results for smaller (0.5*A*_{+,0}: dashed line) and larger (1.5*A*_{+,0}: thick line) values of potentiation amplitude *A*_{+}. The blue lines are for smaller (0.5*A*_{−,0}: dashed line) and larger (1.5*A*_{−,0}: thick line) values of depression amplitude *A*_{−}.

**Figure 8. Effects of changing the STDP learning window**. **(A)** Shapes of the learning windows for various amplitudes of *A*_{+} and *A*_{−}. Black lines represent the standard parameter set (Bi and Poo, 1998), while red and blue lines indicate changing values of potentiation *A*_{+} and depression *A*_{−}, respectively, by a factor 0.5 (dashed lines) or 1.5 (thick lines). **(B,C)** Dependence of average synchronization quality and convergence promptness, respectively, on the number of learning sessions. After about 25–40 sessions, the standard learning window (black lines) becomes optimal.

When depression dominates over potentiation (dashed red and solid blue lines), the SQ is poor even after many sessions, see Figure 8B. Larger values for potentiation relative to depression give rise to faster and better synchrony (dashed blue and solid red lines). However, after about 25–40 learning sessions, the standard set of parameters (Bi and Poo, 1998) yields better results, i.e., higher values of average SQ and faster convergence to synchronization, see Figures 8B,C. To explain this, we will consider why the average SQ in region III becomes higher for the standard set of STDP parameters than for the potentiation-dominated parameter sets (dashed blue and solid red lines in Figure 8A). For the other regions (I, II, and IV) both parameter sets yield similar values of average SQ, because there the system converges to the DS mode after several learning sessions.

We first consider the effects of STDP on the coupling strengths of the network starting with initial coupling strength ε and delay time τ in region III just below the line which separates region II and III. STDP will increase the weights ε to larger values due to larger potentiation relative to depression. This moves (ε, τ) from region III to region II, where the zero-lag synchronization modes “SS1” of Eq. A12 in Appendix and “NDS2” of Eq. A13 in Appendix are unstable. Hence the average SQ in region III decreases because the initial weights ε will be increased to values that correspond to unstable zero-lag synchronization modes in region II. However, when the initial weights and the delay times are farther away from the line separating region II and III, the effects of STDP on the weights are different. To understand this, we focus on pacemaker synchrony, a stable synchronization mode in region III. The period of oscillations for pacemaker synchrony in region III becomes longer when the weights ε become smaller, since the delay between the arrival of spike input and firing increases, cf. formula “PS1” in Eq. A10 in Appendix. Since the time window for depression is much larger than that for potentiation, longer delays suppress potentiation more, leaving depression dominant. Close to pacemaker synchrony, this concerns mainly the synapses from the outer oscillators onto the relay neuron, since the outer neurons fire quickly after receiving input. A weakening of these synapses reinforces pacemaker synchrony. Overall then, increasing the strength of potentiation beyond the standard rule will shift more (ε, τ) to region II without stable synchronization, leaving fewer (ε, τ) that will achieve stable pacemaker synchrony. Hence the SQ and CP in region III eventually becomes higher for the standard STDP set than for the potentiation-dominated one; and the average over all regions follows this trend, see Figures 8B,C.

Results for different time constants for potentiation τ_{+} and depression τ_{−} (not shown) are qualitatively similar to the results for changing the amplitudes *A*_{+} and *A*_{−}, i.e., a larger value of τ_{+} relative to τ_{−} gives rise to faster and better synchrony, and the standard values for τ_{+} and τ_{−} yield better results for zero-lag synchrony when the number of the learning sessions is large. Our results are robust for decreases of depression (dashed blue line) or increases in potentiation (solid red line), but not vice versa (dashed red and solid blue lines). Figure 9 shows the SQ and CP for the model with STDP with HH neurons with instantaneous synapses (Figures 9A,C, respectively) and with alpha synapses (Figures 9B,D). Comparing Figure 9A with Figure 4A reveals hardly any improvement of the SQ by STDP. However, for alpha synapses STDP significantly improves the average SQ by increasing the synaptic coupling ε (compare Figures 9B and 4C).

**Figure 9. Spike-timing dependent plasticity improves long-range synchronization for the model with HH neurons**. **(A,B)** Synchronization quality after the 60-th session for instantaneous synapses and the τ_{syn} = 2 ms alpha synapse, respectively. **(C,D)** Corresponding convergence promptness. For conduction delays approximately τ ≥ 0.1, synchrony improves due to STDP-mediated increase of synaptic strength ε.

In summary, the conditions that contribute to zero-lag long range synchronization are as follows. First, the delay times relative to the intrinsic frequency *T*_{0} should be long enough, i.e., more than about 0.35 (8.75 ms for *T*_{0} = 25ms) for the model with MS phase oscillators with the instantaneous synapse (cf. Figures 2A,C) and 0.15 (2.2 ms for *T*_{0} = 14.66 ms) for the model with HH neurons with the alpha synapse (cf. Figure 4). Moreover, alpha synapses with short synaptic rise times are required for zero-lag synchronization for short delay times (cf. Figure 5). Third, STDP facilitates the synchronization but generally only after many cycles of weight adaptation (cf. Figures 6–9).

### When Delay Times from the Relay Neuron to Outer Neurons are not Identical

We now investigate the dynamics of the model in Figure 1A, when the delay times τ_{1} and τ_{3} are different. Simulation results (not shown) show that zero-lag synchrony then is lost, in agreement with Fischer et al. (2006) and Vicente et al. (2008). Similar results (not shown) were found for the model with HH neurons. To understand why zero-lag synchronization disappears when the delay times are different, we consider the MS network with instantaneous synapses. Let us assume that oscillators 1 and 3 fire simultaneously. When oscillator 2 generates an action potential, the spikes arrive at oscillators 1 and 3 after delays τ_{1} and τ_{3}, respectively. When the spikes arrive at oscillators 1 and 3 at a phase exceeding the critical phase, oscillators 1 and 3 will spike immediately after arrival of the spike. This implies that oscillators 1 and 3 spike with a time difference $\left|\phantom{\rule{2.77695pt}{0ex}}{t}_{1,n+1}^{\text{spike}}-{t}_{3,n+1}^{\text{spike}}\right|=\left|{\tau}_{\text{1}}-{\tau}_{3}\right|\phantom{\rule{2.77695pt}{0ex}}{T}_{\text{0}}$. This illustrates that zero-lag synchrony is lost when the delay times τ_{1} and τ_{3} are different. Qualitatively similar results (not shown) are obtained for asymmetric coupling strengths.

Figure 10 shows a histogram of the relative phase φ* _{r}* for various stable modes for the synaptic weight ε = 0.1 for five pairs of conduction delays: (τ

_{1}, τ

_{3}) = (0.35, 0.25) in Figure 10A, (0.25, 0.15) in Figure 10B, (0.25, 0.25) in Figure 10C, (0.25, 0.35) in Figure 10D, and (0.15, 0.25) in Figure 10E. Figure 10A shows two peaks of relative phases near −0.19 and −0.14, instead of the three peaks appearing in Figure 10C, which include one at zero. Figures 10B,D,E also yield two non-zero peaks instead of three. Figures 10A,B are mirror images of Figures 10D,E, respectively, because the delay times are exchanged and the outer neuron with the shorter delay time spikes first. Essentially, the symmetric non-zero peaks in Figure 10C collapse into one for asymmetric delays, and the zero-lag peak shifts to a non-zero value, which is equal to the value of the zero-lag peak. Other combinations of the synaptic weight ε and non-equal conduction delays τ

_{1}and τ

_{3}yield qualitatively similar results.

**Figure 10. Effect of unequal delays on relative phase**. **(A–E)** Histograms showing the fraction of initial phases that converge to a particular relative phase ε = 0.1 for different pairs (τ_{1}, τ_{3}) = (0.35, 0.25), (0.25, 0.15), (0.25, 0.25), (0.25, 0.35), and (0.15, 0.25). When delays between the outer oscillator and the relay neuron become asymmetric, one of the non-zero relative phases disappears and the zero relative phase shifts to non-zero values, corresponding to a state where the outer oscillator with the shorter delay spikes earlier.

## Discussion

In this study we have investigated the conditions for zero-lag synchrony between two neuronal oscillators, which interact via a relay oscillator. The main result of our study is that for the model with type II Hodgkin–Huxley neurons, synchronization is easier to achieve than for type I Mirollo–Strogatz neurons. Synapses with short rise times (typically less than 2 ms) are more suitable to achieve zero-lag synchronization than synapses with longer rise times. With STDP the network converges to zero-lag synchronization at a faster rate and for a larger range of synaptic strengths and time delays than without STDP. However, when the delay times between the two synchronizing oscillators and the relay oscillator are different, zero phase lag may easily get lost.

The network used in this study is a simplified model for interacting neuronal populations. This obviously raises the question whether our results about zero-lag synchrony may be biased by the simplifications inherent in our model. We will argue in the next paragraphs below that this is not the case. Our choice of indirect interactions between oscillating neuronal populations, i.e., via a relay oscillator, was inspired by previous studies, which showed that pulse-coupled neuronal oscillators with direct excitatory coupling and signal delays in general do not oscillate at zero phase lag (van Vreeswijk et al., 1994; Ernst et al., 1995, 1998; Knoblauch and Sommer, 2003; Zeitler et al., 2009), unless the neurons are of type II with biphasic PRC’s (Goel and Ermentrout, 2002; Woodman and Canavier, 2011). Inhibitory coupling between directly interacting oscillators can cause near zero-lag synchrony (van Vreeswijk et al., 1994; Zeitler et al., 2009). However, the dominant connectivity between cortical areas, such as V1, V2, V4, and FEF, is excitatory, rather than inhibitory. These considerations led Fischer et al. (2006) and Vicente et al. (2008) to postulate a network model of oscillators interacting via a relay oscillator, which supports zero-lag synchrony. Subcortical structures like the thalamus are good physiological candidates for such a mediating relay (Theyel et al., 2010). Our study elaborates on this relay network model.

The most simplified version of our model assumes that the oscillators used to represent neuronal population activity are of the Mirollo–Strogatz type. The Mirollo–Strogatz oscillator corresponds to the type I neuron class (Izhikevich, 2007). Although there is evidence that cortical pyramidal cells can switch between type I and type II by means of cholinergic modulation (Ermentrout et al., 2001; Jeong and Gutkin, 2007; Prescott et al., 2008; Stiefel et al., 2008, 2009), the majority of pyramidal cells in neocortex seems to be type I neurons (Reyes and Fetz, 1993a,b; Kawaguchi, 1995; Erisir et al., 1999; Tateno et al., 2004). Our results show that synchronization is hard to achieve for the model with type I MS neurons for weak and strong synaptic coupling strengths, unless the delays are relatively long. In order to appreciate this result, one should realize that the literature on this topic is divided in studies assuming weak coupling using infinitesimal Phase Response Curves and others assuming strong synaptic coupling. Our results are in agreement with the results of Ermentrout (1996), who used a perturbation method, which is equivalent to assuming weak coupling, for networks of type I neurons with excitatory coupling. However, Mirollo and Strogatz (1990) reported that for almost all initial conditions, a network with strongly coupled type I neurons (without delays!) evolves to a state with synchronous firing. When delays are involved, zero-lag synchrony is lost (Ernst et al., 1995, 1998). Recently, Wang et al. (2012) showed that synchrony in a network with strongly coupled type I oscillators is possible in the absence of delays or with delays greater than half the network period. Although we did not study the model with delays exceeding 0.5, our results reveal the largest amount of synchrony for long delays, which is in agreement with the results by Wang et al. (2012).

Replacing the Mirollo–Strogatz neurons by more realistic type II Hodgkin–Huxley neurons allows for a broader range of synaptic strengths and time delays that is compatible with zero-lag synchronization. This result suggests a functional role for changes in neuronal properties from type I to type II, in agreement with suggestions by Prescott et al. (2008). However, there are no experimental data available yet which can be used to test the hypothesis that the properties of pyramidal cells change from type I to type II when synchrony arises in neuronal populations, as far as we know. Another possibility might be that properties of pyramidal cells do not change from type I to type II, but that the activity of interneurons causes strong inhibition after firing of the pyramidal cells (see Börgers et al., 2010, for a more extended description of the effect of inhibition). In that case, the inhibition by the interneurons after firing causes effectively a biphasic PRC with phase delays early in the firing cycle and phase advances later in the cycle for the pyramidal cell/interneuron couple, which contributes to synchrony both for direct coupling between two pyramidal cells, as well as for a pyramidal cells interacting via a relay neuron.

Delays in networks of interacting neurons can give a variety of complex behaviors with a wealth of bifurcations and a rich phase diagram, which includes oscillatory bumps, traveling waves, lurching waves, standing waves arising via a period-doubling bifurcation, aperiodic regimes, and regimes of multistability (Roxin et al., 2005). Synchronous firing is just one of these modes, which only occurs for a limited range of model parameters. A neuronal property, which greatly contributes to synchrony is that the phase response curve of the neurons is biphasic with phase delays early in the firing cycle and phase advances later in the cycle, like for type II Hodgkin–Huxley neurons. This applies both to weakly coupled oscillators (see e.g., Hansel and Mato, 1993; Hansel et al., 1993a,b) as well as for strongly coupled ones (Bressloff and Coombes, 2000; Izhikevich, 2007). Their results are in agreement with the results in our study, which show that the model with type II Hodgkin–Huxley neurons more easily leads to zero-lag synchronization than with type I MS neurons.

Another assumption of this study, which requires some more discussion, is that all oscillators in the model have identical intrinsic properties with the same oscillation period. If the intrinsic periods of the outer oscillators differ, zero-lag synchrony may get lost. Whether or not synchrony will be lost depends on the neuron type. If the synaptic input to the outer neurons, which have different natural frequencies, resets their oscillation periods to the same value, zero-lag synchronization is easily obtained. Synchronization of non-linear oscillators with different oscillatory properties is feasible if the interactions between the oscillators (the synaptic strengths in our study) are sufficiently strong (Pikovsky et al., 2001). A special condition is the situation where the intrinsic period of the relay oscillator is different from that of the outer oscillators. This might not be unusual if the relay oscillators are thalamic cells and the outer oscillators cortical cells. In this case, the synchronization properties change quantitatively but not qualitatively. These results have been confirmed by simulations (not shown) but can be understood from the following: If we make the period of the relay oscillator different from that of the outer oscillators, the combinations of synaptic strength and delay where input can elicit spikes immediately after arrival change only slightly with adjustments of the period. Hence the boundaries of regions with driven, pacemaker, and slave synchrony will change quantitatively, but not qualitatively, unless the differences in the period become too large.

If the periods of the two outer oscillators are different from each other and if input from the relay oscillator does not make the period of the outer oscillators the same, input from the relay oscillator to the outer neurons will elicit spikes at different times. In that case, the spike input from the outer oscillators to the relay neuron also arrives at different times. This is essentially equivalent to the situation with different delay times, which we have studied (see When Delay Times from the Relay Neuron to Outer Neurons are not identical), where we have shown that zero-lag synchrony is easily lost if the delay times become different (see also Figure 10). Therefore, synaptic coupling strengths should be sufficiently strong to ensure zero-lag synchrony when the oscillation periods of the outer oscillators differ.

In our study we have introduced “SQ” as a measure for the robustness of synchrony against variations in the initial phases of the oscillators. SQ was used together with the “CP” to assess zero-lag synchrony between the two outer oscillators in our relay network. The time interval for “synchrony” was chosen as spike coincidence within 0.5 ms (Engel and Singer, 2001), which in experimental settings is long enough to take into account typical noise on spike-timings but short enough to speak about “zero-lag.” Increasing the time interval to 1 ms did not affect the results qualitatively, though quantitatively some minor differences were observed.

In agreement with Knoblauch and Sommer (2003), we found that if STDP adapts synaptic coupling, the network state converges more easily to a stable state with zero-lag synchrony (see Figures 6 and 9). However, adaptation of the synapses by STDP often took quite some time [in general more than 500 cycles (one session corresponds to 60 cycles), see Figure 7], which implies that STDP may not always play a dominant role for the rapid development of zero-lag synchronization. Vicente et al. (2008) reported that the mechanism of synchronization rests on the ability of an excitatory postsynaptic potential to modify the firing latencies of a postsynaptic neuron in a consistent manner. We agree with this conclusion, but our results show that the mechanism of STDP may take too much time (considerably more than the observed time range of 200–250 ms in visual perception, see Rodriguez et al. (1999)), to generate zero-lag synchrony for oscillations in the gamma frequency range.

Overall, our results demonstrate that gamma oscillations in various cortical areas can be synchronized at zero-lag in a network model where neuronal oscillators are coupled via a relay oscillator, in agreement with previous studies (Fischer et al., 2006; Vicente et al., 2008). In addition we show that STDP expands the range of parameter values, which allow zero-lag synchrony.

## Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

## Acknowledgments

This project was partly funded by the Netherlands Organization for Scientific Research (NWO 051.02.050).

## References

Bi, G.-Q., and Poo, M.-M. (1998). Synaptic modifications in cultured hippocampal neurons: dependence on spike timing, synaptic strength, and postsynaptic cell type. *J. Neurosci.* 18, 10464–10472.

Börgers, C., Krupa, M., and Gielen, S. (2010). The response of a classical Hodgkin-Huxley neuron to an inhibitory input pulse. *J. Comput. Neurosci.* 28, 509–526.

Bressloff, P. C., and Coombes, S. (1998). Desynchronization, mode locking, and bursting in strongly coupled integrate-and-fire oscillators. *Phys. Rev. Lett.* 81, 2168.

Bressloff, P. C., and Coombes, S. (2000). Dynamics of strongly coupled spiking neurons. *Neural Comput.* 12, 91–129.

Buonomano, D. V., and Maass, W. (2009). State-dependent computations: spatiotemporal processing in cortical networks. *Nat. Rev. Neurosci.* 10, 113–125.

Castelo-Branco, M., Neuenschwander, S., and Singer, W. (1998). Synchronization of visual responses between the cortex, lateral geniculate nucleus, and retina in the anesthetized cat. *J. Neurosci.* 18, 6395–6410.

Engel, A. K., Konig, P., Kreiter, A. K., Schillen, T. B., and Singer, W. (1992). Temporal coding in the visual-cortex – new vistas on integration in the nervous-system. *Trends Neurosci.* 15, 218–226.

Engel, A. K., and Singer, W. (2001). Temporal binding and the neural correlates of sensory awareness. *Trends Cogn. Sci. (Regul. Ed.)* 5, 16–25.

Erisir, A., Lau, D., Rudy, B., and Leonard, C. S. (1999). Function of specific K+ channels in sustained high-frequency firing of fast-spiking neocortical interneurons. *J. Neurophysiol.* 82, 2476–2489.

Ermentrout, B. (1996). Type I membranes, phase resetting curves, and synchrony. *Neural Comput.* 8, 979–1001.

Ermentrout, B., Pascal, M., and Gutkin, B. (2001). The effects of spike frequency adaptation and negative feedback on the synchronization of neural oscillators. *Neural Comput.* 13, 1285–1310.

Ernst, U., Pawelzik, K., and Geisel, T. (1995). Synchronization induced by temporal delays in pulse-coupled oscillators. *Phys. Rev. Lett.* 74, 1570–1573.

Ernst, U., Pawelzik, K., and Geisel, T. (1998). Delay-induced multistable synchronization of biological oscillators. *Phys. Rev. E* 57, 2150–2162.

Fischer, I., Vicente, R., Buldú, J. M., Peil, M., Mirasso, C. R., Torrent, M. C., and García-Ojalvo, J. (2006). Zero-lag long-range synchronization via dynamical relaying. *Phys. Rev. Lett.* 97, 123901.

Frien, A., Eckhorn, R., Bauer, R., Woelbern, T., and Kehr, H. (1994). Stimulus-specific fast oscillations at zero phase between visual areas V1 and V2 of awake monkey. *Neuroreport* 5, 2273–2277.

Fries, P. (2005). A mechanism for cognitive dynamics: neuronal communication through neuronal coherence. *Trends Cogn. Sci. (Regul. Ed.)* 9, 474–480.

Fries, P., Reynolds, J. H., Rorie, A. E., and Desimone, R. (2001). Modulation of oscillatory neuronal synchronization by selective visual attention. *Science* 291, 1560–1563.

Froemke, R. C., Tsay, I. A., Raad, M., Long, J. D., and Dan, Y. (2006). Contribution of individual spikes in burst-induced long-term synaptic modification. *J. Neurophysiol.* 95, 1620–1629.

Goel, P., and Ermentrout, B. (2002). Synchrony, stability, and firing patterns in pulse-coupled oscillators. *Physica D* 163, 191–216.

Goldbeter, A. (1996). *Biochemical Oscillations and Cellular Rhythms: The Molecular Bases of Periodic and Chaotic Behaviour*. Cambridge: Cambridge University Press.

Gollo, L. L., Mirasso, C., and Villa, A. E. P. (2010). Dynamic control for synchronization of separated cortical areas through thalamic relay. *Neuroimage* 52, 947–955.

Gray, C. M., Konig, P., Engel, A. K., and Singer, W. (1989). Oscillatory responses in cat visual-cortex exhibit inter-columnar synchronization which reflects global stimulus properties. *Nature* 338, 334–337.

Gross, J., Schmitz, F., Schnitzler, I., Kessler, K., Shapiro, K., Hommel, B., and Schnitzler, A. (2004). Modulation of long-range neural synchrony reflects temporal limitations of visual attention in humans. *Proc. Natl. Acad. Sci. U.S.A.* 101, 13050–13055.

Hansel, D., and Mato, G. (1993). Patterns of synchrony in a heterogeneous Hodgkin-Huxley neural network with weak coupling. *Physica A* 200, 662–669.

Hansel, D., Mato, G., and Meunier, C. (1993a). Clustering and slow switching in globally coupled phase oscillators. *Phys. Rev. E* 48, 3470–3477.

Hansel, D., Mato, G., and Meunier, C. (1993b). Phase dynamics for weakly coupled Hodgkin-Huxley neurons. *Europhys. Lett.* 23, 367.

Hennequin, G., Gerstner, W., and Pfister, J. P. (2010). STDP in adaptive neurons gives close-to-optimal information transmission. *Front. Comput. Neurosci.* 4:143. doi:10.3389/fncom.2010.00143

Hodgkin, A. L., and Huxley, A. F. (1952). A quantitative description of membrane current and its application to conduction and excitation in nerve. *J. Physiol. (Lond.)* 117, 500–544.

Izhikevich, E. M. (2007). *Dynamical Systems in Neuroscience: the Geometry of Excitability and Bursting*. Cambridge, MA: MIT Press.

Jeong, H. Y., and Gutkin, B. (2007). Synchrony of neuronal oscillations controlled by GABAergic reversal potentials. *Neural Comput.* 19, 706–729.

Kawaguchi, Y. (1995). Physiological subgroups of nonpyramidal cells with specific morphological-characteristics in layer Ii/Iii of rat frontal-cortex. *J. Neurosci.* 15, 2638–2655.

Knoblauch, A., and Sommer, F. T. (2003). Synaptic plasticity, conduction delays, and inter-areal phase relations of spike activity in a model of reciprocally connected areas. *Neurocomputing* 52–54, 301–306.

Koch, C., and Segev, I. (1998). *Methods in Neuronal Modeling: from Ions to Networks*. Cambridge, MA: MIT Press.

Lee, S., Sen, K., and Kopell, N. (2009). Cortical gamma rhythms modulate NMDAR-mediated spike timing dependent plasticity in a biophysical model. *PLoS Comput. Biol.* 5, e1000602. doi:10.1371/journal.pcbi.1000602

Lindner, B., Gangloff, D., Longtin, A., and Lewis, J. E. (2009). Broadband coding with dynamic synapses. *J. Neurosci.* 29, 2076–2088.

Mirollo, R. E., and Strogatz, S. H. (1990). Synchronization of pulse-coupled biological oscillators. *SIAM J. Appl. Math.* 50, 1645–1662.

Pesaran, B., Pezaris, J. S., Sahani, M., Mitra, P. P., and Andersen, R. A. (2002). Temporal structure in neuronal activity during working memory in macaque parietal cortex. *Nat. Neurosci.* 5, 805–811.

Pikovsky, A., Rosenblum, M., and Kurths, J. (2001). *Synchronization: A Universal Concept in Nonlinear Sciences*. Cambridge: Cambridge University Press.

Prescott, S. A., Ratté, S., De Koninck, Y., and Sejnowski, T. J. (2008). Pyramidal neurons switch from integrators in vitro to resonators under in vivo-like conditions. *J. Neurophysiol.* 100, 3030–3042.

Reyes, A. D., and Fetz, E. E. (1993a). Effects of transient depolarizing potentials on the firing rate of cat neocortical neurons. *J. Neurophysiol.* 69, 1673–1683.

Reyes, A. D., and Fetz, E. E. (1993b). Two modes of interspike interval shortening by brief transient depolarizations in cat neocortical neurons. *J. Neurophysiol.* 69, 1661–1672.

Rodriguez, E., George, N., Lachaux, J. P., Martinerie, J., Renault, B., and Varela, F. J. (1999). Perception’s shadow: long-distance synchronization of human brain activity. *Nature* 397, 430–433.

Roelfsema, P. R., Engel, A. K., Konig, P., and Singer, W. (1997). Visuomotor integration is associated with zero time-lag synchronization among cortical areas. *Nature* 385, 157–161.

Roxin, A., Brunel, N., and Hansel, D. (2005). Role of delays in shaping spatiotemporal dynamics of neuronal activity in large networks. *Phys. Rev. Lett.* 94, 1–4.

Schoffelen, J. M., Oostenveld, R., and Fries, P. (2005). Neuronal coherence as a mechanism of effective corticospinal interaction. *Science* 308, 111–113.

Singer, W., and Gray, C. M. (1995). Visual feature integration and the temporal correlation hypothesis. *Annu. Rev. Neurosci.* 18, 555–586.

Song, S., Miller, K. D., and Abbott, L. F. (2000). Competitive Hebbian learning through spike-timing-dependent synaptic plasticity. *Nat. Neurosci.* 3, 919–926.

Stiefel, K. M., Gutkin, B. S., and Sejnowski, T. J. (2008). Cholinergic neuromodulation changes phase response curve shape and type in cortical pyramidal neurons. *PLoS ONE* 3, e3947. doi:10.1371/journal.pone.0003947

Stiefel, K. M., Gutkin, B. S., and Sejnowski, T. J. (2009). The effects of cholinergic neuromodulation on neuronal phase-response curves of modeled cortical neurons. *J. Comput. Neurosci.* 26, 289–301.

Tateno, T., Harsch, A., and Robinson, H. P. C. (2004). Threshold firing frequency-current relationships of neurons in rat somatosensory cortex: Type 1 and type 2 dynamics. *J. Neurophysiol.* 92, 2283–2294.

Theyel, B. B., Llano, D. A., and Sherman, M. (2010). The corticothalamocortical circuit drives higher-order cortex in the mouse. *Nat. Neurosci.* 13, 84–88.

Tsubo, Y., Takada, M., Reyes, A. D., and Fukai, T. (2007). Layer and frequency dependencies of phase response properties of pyramidal neurons in rat motor cortex. *Eur. J. Neurosci.* 25, 3429–3441.

van Rossum, M. C. W., Bi, G. Q., and Turrigiano, G. G. (2000). Stable Hebbian learning from spike timing-dependent plasticity. *J. Neurosci.* 20, 8812–8821.

van Vreeswijk, C., Abbott, L. F., and Ermentrout, G. B. (1994). When inhibition not excitation synchronizes neural firing. *J. Comput. Neurosci.* 1, 313–321.

Vicente, R., Gollo, L. L., Mirasso, C. R., Fischer, I., and Pipa, G. (2008). Dynamical relaying can yield zero time lag neuronal synchrony despite long conduction delays. *Proc. Natl. Acad. Sci. U.S.A.* 105, 17157–17162.

Wang, S., Chandrasekaran, L., Fernandez, F. R., White, J. A., and Canavier, C. C. (2012). Short conduction delays cause inhibition rather than excitation to favor synchrony in hybrid neuronal networks of the entorhinal cortex. *PLoS Comput. Biol.* 8, e1002306. doi:10.1371/journal.pcbi.1002306

Woodman, M. M., and Canavier, C. C. (2011). Effects of conduction delays on the existence and stability of one to one phase locking between two pulse-coupled oscillators. *J. Comput. Neurosci.* 31, 401–418.

Zeitler, M., Daffertshofer, A., and Gielen, C. C. (2009). Asymmetry in pulse-coupled oscillators with delay. *Phys. Rev. E Stat. Nonlin. Soft Matter Phys.* 79, 065203.

## Appendix

### Detailed Analysis of Phase-Locking and Stability

Analytical derivations of the various synchrony modes and their stability are provided. An overview of these results is provided by Figure A1.

#### Synchrony modes derived from the phase-locking equation

To investigate the existence of zero-lag synchrony as a function of τ and ε, we will derive the 1:1 synchronized phase-locking equation of the oscillators, cf. Figure 1B. The firing of the relay neuron ${t}_{2,n-1}^{\text{spike}}$ is chosen as reference time, i.e., ${t}_{2,n-1}^{\text{spike}}=0.$ For zero-lag synchrony there exists some value þeta ∈ [0, 1) such that the firing times ${t}_{1,n-1}^{\text{spike}}$ and ${t}_{3,n-1}^{\text{spike}}$ of both neurons 1 and 3 are equal to θ*T*, and ${t}_{1,n}^{\text{spike}}$ and ${t}_{3,n}^{\text{spike}}$ are equal to θ*T* + *T*, with *T* the period of the oscillations in the network. The spikes from the outer oscillators arrive at the relay oscillator at phase τ + θ*T*/*T*_{0}, if τ ≤ (1 − θ)*T*/*T*_{0}, or at τ − (1 − θ)*T*/*T*_{0}, if τ > (1 − θ)*T*/*T*_{0}. The new period of the inner oscillator 2 is hence given by

where Δφ is defined in Eqs 6 and 7 of the main text, repeated here for convenience:

**Figure A1. Synchronization regions in the parameter plane of delay time and coupling strength**. Different combinations of possible synchrony modes divide the parameter plane of delay time τ and synaptic weight ε into ten distinct regions. These are shown separated by black lines and the corresponding synchrony modes are indicated by symbols referring to results listed in the Appendix: driven synchrony – Eq. A7 “DS” in Appendix; pacemaker synchrony – Eq. A10 “PS1” and “PS2” in Appendix; slave synchrony – Eq. A12 “SS1” and “SS2” in Appendix; and non-driven synchrony – Eq. A13 “NDS1,” “NDS2,” and “NDS3” in Appendix.

The spike from oscillator 2 will arrive at oscillators 1 and 3 at phase τ + (1 − θ)*T*/*T*_{0}, if τ ≤ θ*T*/*T*_{0}, or τ − θ*T*/*T*_{0}, if τ > θ*T*/*T*_{0}. The new period of the outer oscillators 1 and 3 is hence given by

Equations A1 and A4 are the 1:1 phase-locking equations, which give a relation between the synaptic weight ε, the conduction delay τ, and the new period of the oscillators *T*.

Using Eqs A1 and A4, we will derive four synchronization modes: driven, pacemaker, slave, and non-driven synchrony. In driven synchrony spikes from the outer neurons immediately initiate a spike of the relay neuron, and vice versa, cf. top of Figure 2B. In pacemaker synchrony arrival of a spike from the relay oscillator at the outer oscillators immediately elicits a spike, but not vice versa, cf. bottom of Figure 2B. When the relay neuron spikes immediately after spikes from the outer neurons, but not vice versa, we call this slave synchrony, cf. middle of Figure 2B. Finally, non-driven synchrony occurs when neither the relay neuron nor the outer neurons spike immediately after input.

First, we consider driven synchrony. In order for all input spikes to elicit a spike in the receiving oscillators, we must have the following: for the inner oscillator φ* _{c}*(2ε) ≤ τ + θ

*T*/

*T*

_{0}if τ ≤ (1 − θ)

*T*/

*T*

_{0}, and φ

*(2ε) ≤ τ − (1 − θ)*

_{c}*T*/

*T*

_{0}if τ > (1 − θ)

*T*/

*T*

_{0}; and for the outer oscillators φ

*(ε) ≤ τ + (1 − θ)*

_{c}*T*/

*T*

_{0}, if τ ≤ θ

*T*/

*T*

_{0}, and φ

*(ε) ≤ τ − θ*

_{c}*T*/

*T*

_{0}, if τ > θ

*T*/

*T*

_{0}. Using Eqs A1, A2, and A4, we can then rewrite the new period of the inner oscillator as

and the new period of the outer oscillators as

We can now solve for θ and *T* by using one part of Eqs A5 and A6 each, and combining them under the condition that the new periods *T* of the inner and outer oscillators, respectively, must be identical to achieve synchrony. This yields four possible combinations:

For biologically realistic conditions τ ∈ (0, 0.5) and ε ∈ (0, 0.21) and our choice *b* = 3, only the driven synchrony mode labeled “DS” turns out to be valid. The firing dynamics of “DS” is illustrated in the top of Figure 2B. For the other cases the φ* _{c}* conditions are not fulfilled and therefore these solutions are rejected. Moreover, the third solution also has a θ value out of range for 1:1 phase-locking. We will similarly identify the valid modes of other solutions below. Since φ

*(2ε) ≤ φ*

_{c}*(ε), φ*

_{c}*(ε) ≤ 2τ implies φ*

_{c}*(2ε) ≤ 2τ and the φ*

_{c}*(ε) ≤ 2τ constraint is sufficient. This constraint divides the ε and τ parameter space into two regions: driven synchrony is possible in the region where φ*

_{c}*(ε) ≤ 2τ but not in the region where φ*

_{c}*(ε) > 2τ. These two regions are separated by the line φ*

_{c}*(ε) = 2τ illustrated as the line separating regions I and II in Figures 2A,C.*

_{c}Second, we investigate pacemaker synchrony. For pacemaker synchrony only the outer oscillators spike directly upon input, but not the inner one. Thus for the inner oscillator, we must have φ* _{c}*(2ε) > τ + θ

*T*/

*T*

_{0}if τ ≤ (1 − θ)

*T*/

*T*

_{0}, and φ

*(2ε) > τ − (1 − θ)*

_{c}*T*/

*T*

_{0}if τ > (1 − θ)

*T*/

*T*

_{0}. For the outer oscillators, we must have φ

*(ε) ≤ τ + (1 − θ)*

_{c}*T*/

*T*

_{0}if τ ≤ θ

*T*/

*T*

_{0}, and φ

*(ε) ≤ τ − θ*

_{c}*T*/

*T*

_{0}if τ > θ

*T*/

*T*

_{0}. Using Eq. A2 in Appendix, we can rewrite the new period of the inner oscillator given by Eq. A1 in Appendix as

where as the new period of the outer oscillators is given by Eq. A6 in Appendix. Solving for θ and *T* as before, we find the following four possibilities:

The two valid modes for pacemaker synchrony given our constraints are labeled “PS1” and “PS2” in Eq. A10 in Appendix, respectively. Within our ranges of interest for ε and τ, one can check numerically for the “PS1” and “PS2” pacemaker synchrony that φ* _{c}*(ε) ≤

*T*/

*T*

_{0}implies φ

*(2ε) > 2τ and φ*

_{c}*(2ε) > 2τ −*

_{c}*T*/

*T*

_{0}, respectively. Hence, only the φ

*(ε) ≤*

_{c}*T*/

*T*

_{0}constraint is required. “PS1” pacemaker synchrony is possible in the region {(ε,τ) ∈ [0, 0.21] × [0, 0.5]|φ

*(ε) ≤*

_{c}*T*/

*T*

_{0}}; while “PS2” is possible in the region {(ε,τ) ∈ [0, 0.21] × [0, 0.5]|φ

*(ε) ≤*

_{c}*T*/

*T*

_{0}and (1 − θ)

*T*/

*T*

_{0}< τ}, where (1 − θ)

*T*/

*T*

_{0}< τ refers to the condition given in Eq. A1 in Appendix. The equality φ

*(ε) =*

_{c}*T*/

*T*

_{0}of “PS1” gives the line separating regions II and III in Figures 2A,C. The firing dynamics for “PS1” pacemaker synchrony is illustrated in Figure 2B.

Third, in slave synchrony only the inner oscillator spikes upon input, i.e., in case of the inner oscillator, we must have φ* _{c}*(2ε) ≤ τ + θ

*T*/

*T*

_{0}if τ ≤ (1 − θ)

*T*/

*T*

_{0}, and φ

*(2ε) ≤ τ − (1 − θ)*

_{c}*T*/

*T*

_{0}if τ > (1 − θ)

*T*/

*T*

_{0}. For the outer oscillators, we must have φ

*(ε) > τ + (1 − θ)*

_{c}*T*/

*T*

_{0}if τ ≤ θ

*T*/

*T*

_{0}, and φ

*(ε) > τ − θ*

_{c}*T*/

*T*

_{0}if τ > θ

*T*/

*T*

_{0}. The new period of the inner oscillator is given by Eq. A5 in Appendix and from Eqs A2 and A4 we obtain the following expression for the new period of the outer oscillator

The four combinations are then:

The two valid slave synchronization modes will be called “SS1” and “SS2,” respectively. Within our ranges of interest for ε and τ, one can check numerically for the “SS1” slave synchrony mode that φ* _{c}*(ε) > 2τ implies φ

*(2ε) ≤*

_{c}*T*/

*T*

_{0}; whereas the constraints for “SS2” are always fulfilled in our ranges for ε and τ. Hence “SS1” slave synchrony, the firing dynamic of which is illustrated in Figure 2B, is possible in the region {(ε,τ) ∈ [0, 0.21] × [0, 0.5]|φ

*(ε) > 2τ}; while “SS2” is possible in the region {(ε,τ) ∈ [0, 0.21] × [0, 0.5]|τ > θ*

_{c}*T*/

*T*

_{0}}, where τ > θ

*T*/

*T*

_{0}refers to the condition given in Eq. A4 in Appendix. The equality τ = θ

*T*/

*T*

_{0}of “SS2” gives the line separating regions I and IV in Figures 2A,C.

Finally, for non-driven synchrony none of oscillators fires immediately after receiving a spike. This implies that for the inner oscillator φ* _{c}*(2ε) > τ + θ

*T*/

*T*

_{0}if τ ≤ (1 − θ)

*T*/

*T*

_{0}, and φ

*(2ε) > τ − (1 − θ)*

_{c}*T*/

*T*

_{0}if τ > (1 − θ)

*T*/

*T*

_{0}. For the outer oscillators φ

*(ε) > τ + (1 − θ)*

_{c}*T*/

*T*

_{0}if τ ≤ θ

*T*/

*T*

_{0}, and φ

*(ε) > τ − θ*

_{c}*T*/

*T*

_{0}if τ > θ

*T*/

*T*

_{0}. The new period is then given by Eqs A8 and A11 for the inner and outer oscillators, respectively. The resulting combinations are:

with ρ_{mnop} = [*m* + *n*β* _{b}*(ε)]γ

*(2ε) + [*

_{b}*o*+

*p*β

*(2ε)]γ*

_{b}*(ε), where indices can take on the values + 1, 0, −1 as represented by the signs, and likewise ζ*

_{b}*≡ β*

_{m}*(ε) + β*

_{b}*(2ε) +*

_{b}*m*β

*(ε)β*

_{b}*(2ε). The three valid non-driven synchronization modes will be called NDS1, NDS2, and NDS3, respectively. From the constraints of ε and τ for “NDS1” slave synchrony, φ*

_{b}*(ε) > τ − ρ*

_{c}_{+0−−}/ζ

_{+}implies φ

*(2ε) > τ + ρ*

_{c}_{++−0}/ζ

_{+}. For the “NDS2” slave synchrony, φ

*(ε) > τ − ρ*

_{c}_{+0−0}/ζ

_{0}implies φ

*(2ε) > τ + ρ*

_{c}_{+0−0}/ζ

_{0}. But the constraints of “NDS3” fulfill automatically in our range of ε and τ. Hence the “NDS1” slave synchrony is possible in the region {(ε,τ)∈[0, 0.21] × [0, 0.5]|φ

*(ε) > τ − ρ*

_{c}_{+0−−}/ζ

_{+}}; while “NDS2” is possible in the region {(ε,τ) ∈ [0, 0.21] × [0, 0.5]|φ

*(ε) > τ − ρ*

_{c}_{+0−0}/ζ

_{0}and τ > (1 − θ)

*T*/

*T*

_{0}} and “NDS3” in the region {(ε,τ) ∈ [0, 0.21] × [0, 0.5]|τ > (1 − θ)

*T*/

*T*

_{0}}, where τ > (1 − θ)

*T*/

*T*

_{0}refers to the condition given in Eq. A1 in Appendix.

Based on different combinations of the synchronization modes DS from Eq. A7 in Appendix, PS1 and PS2 from Eq. A10 in Appendix, SS1 and SS2 from Eq. A12 in Appendix, and NDS1, NDS2, and NDS3 from Eq. A13 in Appendix, we can separate the parameter plane of time delay τ and synaptic weight ε into ten regions, see Figure A1. In summary, the regions in the ε and τ parameter space, where the modes are possible, are given by

#### Stability analysis of synchrony modes

In order to investigate the stability of driven synchrony, we assume that the fixed point Eq. A7 “DS” in Appendix suffers from small phase perturbations (φ_{1}, φ_{3}) = (τ + δφ_{1}, τ + δφ_{1}) at time ${t}_{2,n-1}^{\text{spike}}.$ Driven synchrony is asymptotically stable, if there is a δ > 0 such that the phases of the outer oscillators will be closer to τ at the next spike time ${t}_{2,n}^{\text{spike}}$ for $\sqrt{\delta {\phi}_{1}^{2}+\delta {\phi}_{3}^{2}}<\delta .$ Since φ* _{c}*(ε) < 2τ, we can define δ ≡ 2τ − φ

*(ε) > 0, and by this definition φ*

_{c}*(ε) ≡ 2τ − δ. The spike from the relay arrives when oscillator 1 has phase 2τ − δ < φ*

_{c}_{1}< 2τ + δ, which exceeds the critical phase φ

*(ε). Therefore, oscillator 1 will spike immediately after input. Since the same argument applies to oscillator 3, the outer oscillators will again spike in synchrony, proving asymptotic stability for small perturbations. Similar arguments can be made to show that the pacemaker synchrony of Eq. A10 “PS1” and “PS2” in Appendix is asymptotically stable.*

_{c}Next we investigate the stability of the slave synchronization mode Eq. A12 “SS1” in Appendix. Here one has to calculate the eigenvalues λ of the Jacobian of the return map near the fixed point of the outer oscillator phases (Zeitler et al., 2009). If all eigenvalues are within the unit circle (|λ| < 1), then the fixed point is asymptotically stable, otherwise it is unstable. Assume at time ${t}_{2,n-1}^{\text{spike}}$ that phases are perturbed away from the fixed point of Eq. A12 “SS1” in Appendix: (φ_{1}, φ_{3}) = (τ + δφ_{1}, τ + δφ_{3}). Thus synaptic inputs from oscillators 1 and 3 will arrive at oscillator 2 at slightly different times. Without loss of generality, we assume that the synaptic input from oscillator 1 arrives earlier. This first synaptic input can make the relay oscillator spike immediately, if φ* _{c}*(ε) < τ + θ

*T*/

*T*

_{0}. Assuming that the perturbation is small enough for the solution in Eq. A12 “SS1” in Appendix to remain applicable, we find φ

*(ε) < τ + θ*

_{c}*T*/

*T*

_{0}= 1 − χ

*(ε) − 2β*

_{b}*(ε)τ = [1 + β*

_{b}*(ε)]φ*

_{b}*(ε) − 2β*

_{c}*(ε)τ. For non-zero coupling this is equivalent to 2τ < φ*

_{b}*(ε), which is biologically reasonable (Tsubo et al., 2007). Therefore we assume φ*

_{c}*(ε) < τ + θ*

_{c}*T*/

*T*

_{0}, and the spike from oscillator 1 will make the relay oscillator spike instantaneously when it arrives. Define ${\stackrel{\u0303}{\phi}}_{1}$ as the phase of oscillator 1 at this relay spike time ${t}_{2,n}^{\text{spike}},$ then obviously ${\stackrel{\u0303}{\phi}}_{1}=\tau ,$ the conduction delay.

To find the phase ${\stackrel{\u0303}{\phi}}_{3}$ of oscillator 3 when oscillator 2 spikes at ${t}_{2,n}^{\text{spike}},$ we first derive the previous spike times of the outer oscillators ${t}_{1,n-1}^{\text{spike}}$ and ${t}_{3,n-1}^{\text{spike}}.$ The phase of oscillator 1 is τ + δφ_{1} when the relay oscillator fires at ${t}_{2,n-1}^{\text{spike}}\equiv 0.$ After a delay τ*T*_{0}, the spike from the relay oscillator arrives at oscillator 1, thus at its phase ${\phi}_{1}^{\text{old}}=2\tau +\delta {\phi}_{1}.$ If δφ_{1} is small enough we can still assume 2τ + δφ_{1} < φ* _{c}*(ε), i.e., this input from the relay will not make oscillator 1 spike instantaneously. Instead its phase will jump to ${\phi}_{1}^{\text{new}}={\phi}_{1}^{\text{old}}+{\chi}_{b}\left(\epsilon \right)+{\beta}_{b}\left(\epsilon \right){\phi}_{1}^{\text{old}},$ cf. Eq. A2 in Appendix. After receiving the spike, it takes an amount of time ${T}_{\text{0}}\left(1-{\phi}_{1}^{\text{new}}\right)$ until oscillator 1 spikes at ${t}_{1,n-1}^{\text{spike}}.$ Hence ${t}_{1,n-1}^{\text{spike}}$is obtained by adding the conduction delay τ

*T*

_{0}between oscillator 2 and 1 to the time${T}_{\text{0}}\left(1-{\phi}_{1}^{\text{new}}\right)$ that oscillator 1 still needs to evolve before spiking. Spike time ${t}_{3,n-1}^{\text{spike}}$ can be obtained in a similar way, and thus one finds:

The next spike of the relay oscillator occurs at ${t}_{2,n}^{\text{spike}}={t}_{1,n-1}^{\text{spike}}+\tau {T}_{\text{0}},$ as explained above. Since ${\tilde{\phi}}_{3}=\left({t}_{2,n}^{\text{spike}}-{t}_{3,n-1}^{\text{spike}}\right)/{T}_{0}=\tau +\left({t}_{1,n-1}^{\text{spike}}-{t}_{3,n-1}^{\text{spike}}\right)/{T}_{0}$, we find the following return map for our local stability analysis

The corresponding Jacobian matrix is

The eigenvalues of this matrix are 0 and 1 + β* _{b}*(ε) > 1. Consequently, this fixed point is not stable. Equation A12 “SS2” in Appendix can be shown to be unstable in a similar manner.

Finally, we investigate the stability of non-driven synchronization Eq. A13 “NDS1” in Appendix. We will discuss the return map near the fixed point $\left({\phi}_{\text{1}}^{*},{\phi}_{\text{3}}^{*}\right)=\left(\left(1-\theta \right)\text{T},\left(1-\theta \right)T\right).$ At time ${t}_{2,n-1}^{\text{spike}},$ we assume a small perturbation $\left({\phi}_{1},{\phi}_{3}\right)=\left({\phi}_{\text{1}}^{*}+\delta {\phi}_{1},{\phi}_{\text{3}}^{*}+\delta {\phi}_{3}\right).$ The arguments used to derive Eqs A15 and A16 hold here, except that the perturbed fixed point is $\tau +\delta {\phi}_{1}\to {\phi}_{\text{1}}^{*}+\delta {\phi}_{1}$ and $\tau +\delta {\phi}_{3}\to {\phi}_{\text{3}}^{*}+\delta {\phi}_{3},$ respectively. Without loss of generality, we assume that the synaptic input from oscillator 1 arrives at oscillator 2 before that from oscillator 3. Thus when oscillator 1 spikes at ${t}_{1,n-1}^{\text{spike}},$ the spike input from oscillator 1 arrives at oscillator 2 at phase ${t}_{1,n-1}^{\text{spike}}/{T}_{\text{0}}+\tau $ and causes the phase of oscillator 2 to change to ${\chi}_{b}\left(\epsilon \right)+\left[\text{1}+{\beta}_{b}\left(\epsilon \right)\right]\left[{t}_{1,n-1}^{\text{spike}}/{T}_{\text{0}}+\tau \right].$ Then the spike generated at ${t}_{3,n-1}^{s\text{pike}}$ by oscillator 3 will arrive when oscillator 2 has phase ${\chi}_{b}\left(\epsilon \right)+\left[\text{1}+{\beta}_{b}\left(\epsilon \right)\right]\left[{t}_{1,n-1}^{\text{spike}}/{T}_{\text{0}}+\tau \right]+\left({t}_{3,n-1}^{\text{spike}}-{t}_{1,n-1}^{\text{spike}}\right)/{T}_{\text{0}}$ and this causes a change in phase of oscillator 2 to

with oscillator 2 spiking after this additional time Δ*t* has passed. Hence, the phases of the outer oscillators at the time ${t}_{2,n}^{\text{spike}},$ when oscillator 2 spikes, are

The corresponding Jacobian matrix is

with resulting eigenvalues

This fixed point is not stable, because these eigenvalues exceed the value one for ε∈(0, 0.21). In a similar manner one can show that Eq. A13 “NDS2” and “NDS3” in Appendix are unstable.

#### Impact of STDP on coupling in slave synchrony

In slave synchrony, the coupling strengths from the outer oscillators to the relay oscillator remain unchanged by STDP since the relay oscillator instantly spikes upon input from the outer oscillators. How about the coupling from the relay oscillator to the outer oscillators? Say their strengths are ε* _{n - 1}* just after the outer oscillators spike at ${t}_{1,n-1}^{\text{spike}}={t}_{3,n-1}^{\text{spike}}.$ The relay spike arrives when the outer oscillators have phase 2τ, since the relay oscillator fired at their phase τ. Considering all prior spikes of the outer oscillators at times ${t}_{1,l}^{\text{spike}},$ this new relay spike arriving at time ${t}_{1,n-1}^{\text{spike}}+2\tau {T}_{\text{0}}$ changes the coupling strength by

Note that this change is a decrease in coupling strength because ${t}_{1,l}^{\text{spike}}<{t}_{1,n-1}^{\text{spike}}+2\tau {T}_{\text{0}}$ for *l* = 0, 1, …, *n* − 1 and therefore $W\left[{t}_{1,l}^{\text{spike}}-\left({t}_{1,n-1}^{\text{spike}}+2\tau {T}_{\text{0}}\right)\right]<0.$ The input from the relay oscillator also causes a phase shift Δφ(2τ|ε* _{n − 1}*) of the outer oscillators such that their next spike occurs at ${t}_{1,n}^{\text{spike}}={t}_{1,n-1}^{\text{spike}}+{T}_{\text{0}}-{T}_{\text{0}}\Delta \phi \left(2\tau \phantom{\rule{2.77695pt}{0ex}}\left|\phantom{\rule{2.77695pt}{0ex}}{\epsilon}_{n-1}\right.\right)$. Hence the new relay spike arrival at the outer oscillators precedes the next outer oscillator spike by (1 − 2τ)

*T*

_{0}−

*T*

_{0}Δφ(2τ|ε

*), and the previous relay spikes arrive at the outer oscillators at ${t}_{1,n-1}^{\text{spike}}-{t}_{1,k}^{\text{spike}}$ for*

_{n - 1}*k*= 0, 1, …,

*n*− 2 earlier than the new relay spike arrival. Therefore, all previous synaptic inputs considered as a whole lead to an increase of the synaptic weight by

Under certain conditions, see main text, this leads to STDP turning slave synchrony into driven synchrony, and the transition time can be determined by numerically iterating the formulas for ${\tilde{\epsilon}}_{n-1}$ and ε* _{n}*.

Keywords: long-range synchronization, spike-timing dependent plasticity, zero-lag synchronization

Citation: Viriyopase A, Bojak I, Zeitler M and Gielen S (2012) When long-range zero-lag synchronization is feasible in cortical networks. *Front. Comput. Neurosci.* **6**:49. doi: 10.3389/fncom.2012.00049

Received: 26 January 2012; Accepted: 27 June 2012;

Published online: 27 July 2012.

Edited by:

Ad Aertsen, University of Freiburg, GermanyReviewed by:

Carmen Canavier, LSU Health Sciences Center, USAChristian Leibold, Ludwig Maximilians University, Germany

Copyright: © 2012 Viriyopase, Bojak, Zeitler and Gielen. This is an open-access article distributed under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in other forums, provided the original authors and source are credited and subject to any copyright notices concerning any third-party graphics etc.

*Correspondence: Atthaphon Viriyopase, Cognitive Neuroscience, UMC St Radboud, Hp 126, Postbox 9101, 6500 HB Nijmegen, Netherlands. e-mail: a.viriyopase@donders.ru.nl