When Long-Range Zero-Lag Synchronization is Feasible in Cortical Networks

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., 2001Fries et al., , 2008Pesaran 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(Ernst et al., , 1998Goel 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 .
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 A B 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, τ k the conduction delay between relay oscillator 2 and oscillator 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 spike 1,n , t spike 2,n , and t spike 3,n are the n-th spike times of oscillators 1, 2, and 3, respectively.
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 τ k represents the conduction time for spikes along axons that connect oscillator 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 τ k expressed as a fraction of the intrinsic period 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 The excitatory synaptic current I syn is −ε ij S(t )V, where ε ij is the maximum synaptic conductance and 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 socalled 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 Frontiers in Computational Neuroscience www.frontiersin.org neuronal state f i of the MS phase oscillator after arrival of a spike from oscillator j is incremented instantaneously with ϕ i the phase of the postsynaptic oscillator i at the time of the spike arrival. The phase ϕ i for which the oscillator reaches the threshold after spike input [i.e., f i (ϕ i ) + ε ij = 1] is called the critical phase ϕ c (Zeitler et al., 2009) and is given by If a spike arrives at ϕ i < ϕ c , f i will increase instantaneously by ε ij . The instantaneous change of the state f i by ε ij corresponds to a phase shift with where ϕ i is the phase of the postsynaptic oscillator just before arrival of the input spike. Note that 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 ϕ c > 0.5. This restriction imposes an upper bound of ε ij < 1 − ln[(1 + e b )/2]/b = 0.21 for b = 3 on the synaptic strength through Eq. 5. In case of ε ij = 0.21, a spike can cause an increase in the state variable 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 arr k and a postsynaptic spike at t spike l the fractional synaptic modification W (∆t ) is given by − t arr k . The spike arrival time t arr k 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.

Frontiers in Computational Neuroscience
www.frontiersin.org

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 zerolag 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 spike 2,n−1 = 0 and call the period of the 1:1 phase-locked oscillators T. For zero-lag synchrony there is a θ ∈ [0, 1) such that the outer oscillators fire at t spike 1,n−1 = t spike 3,n−1 = θT , t spike 1,n = t spike 3,n = θ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 ϕ c (ε), respectively. Then rewriting Eqs 10 and 11 using Eq. 6 yields T = τT 0 + θT and T = τT 0 + (1 − θ) T ⇒ T = 2τT 0 , θ = 1/2.
(12) 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.
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 spike 2,n−1 . 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 spike 2,n for δϕ 2 1 + δϕ 2 3 < δ. Since ϕ c (ε) < 2τ, we can define a value δ ≡ 2τ − ϕ c (ε) > 0, and thus ϕ c (ε) = 2τ − δ. The spike from the relay arrives when the phase of oscillator 1 is in the range 2τ − δ < ϕ 1 < 2τ + δ, which exceeds this critical phase ϕ c (ε). 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.
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 > ε n − 1 for τ ≥ 0.25 using the STDP window W (∆t ) of Eq. 9, since the spike from the relay oscillator arrives well before the outer neurons spike. Hence synaptic weights increase to 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.

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 t spike 1,n − t spike 3,n ≤ 0.02T 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 zerolag 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 zerolag 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 ϕ r ≡ t spike 3,n − t spike 1,n 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.

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 For about 90% of possible initial phases, the system converges to a non-zero relative phase. (B) |ϕ r | 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. (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 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 , χ b , 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 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 (ϕ r = 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 Frontiers in Computational Neuroscience www.frontiersin.org 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.
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 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.05T 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 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.
Frontiers in Computational Neuroscience www.frontiersin.org 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. 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.
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 τ + c ≡ ϕ c (ε max ) 2 ≈ 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 τ + c ≤ τ < ϕ c (ε) 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 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  from the relay to the outer oscillators), respectively. After training for 60 sessions, the values of |ϕ r | 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.
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)  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 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).
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 t spike 1,n+1 − t spike 3,n+1 = |τ 1 − τ 3 | T 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.

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 zerolag 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., 1995Ernst et al., , 1998Knoblauch 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., 2008Stiefel et al., , 2009, the majority of pyramidal cells in neocortex seems to be type I neurons (Reyes and Fetz, 1993a,b; Frontiers in Computational Neuroscience www.frontiersin.org 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. 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(Ernst et al., , 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).

Frontiers in Computational Neuroscience www.frontiersin.org
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 "zerolag." 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.

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

Frontiers in Computational Neuroscience
www.frontiersin.org 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.
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 ϕ c (ε) ≤ τ and ϕ c (2ε) ≤ τ 1 2 , 2 3 τ for ϕ c (ε) ≤ 2 3 τ and ϕ c (2ε) ≤ 2 3 τ .

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 spike 2,n−1 . 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 spike 2,n for δϕ 2 1 + δϕ 2 3 < δ. Since ϕ c (ε) < 2τ, we can define δ ≡ 2τ − ϕ c (ε) > 0, and by this definition ϕ c (ε) ≡ 2τ − δ. The spike from the relay arrives when oscillator 1 has phase 2τ − δ < ϕ 1 < 2τ + δ, which exceeds the critical phase ϕ c (ε). 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.
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 spike 2,n−1 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 − χ b (ε) − 2β b (ε)τ = [1 + β b (ε)]ϕ c (ε) − 2β b (ε)τ. For non-zero coupling this is equivalent to 2τ < ϕ c (ε), which is biologically reasonable (Tsubo et al., 2007). Therefore we assume ϕ c (ε) < τ + θT /T 0 , and the spike from oscillator 1 will make the relay oscillator spike instantaneously when it arrives. Defineφ 1 as the phase of oscillator 1 at this relay spike time t spike 2,n , then obviouslyφ 1 = τ, the conduction delay.
To find the phaseφ 3 of oscillator 3 when oscillator 2 spikes at t spike 2,n , we first derive the previous spike times of the outer oscillators t spike 1,n−1 and t spike 3,n−1 . The phase of oscillator 1 is τ + δϕ 1 when the relay oscillator fires at t spike 2,n−1 ≡ 0. After a delay τT 0 , the spike from the relay oscillator arrives at oscillator 1, thus at its phase ϕ old 1 = 2τ + δϕ 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 ϕ new 1 = ϕ old 1 +χ b (ε)+β b (ε) ϕ old 1 , cf. Eq. A2 in Appendix. After receiving the spike, it takes an amount of time T 0 1 − ϕ new 1 until oscillator 1 spikes at t spike 1,n−1 . Hence t spike 1,n−1 is obtained by adding the conduction delay τT 0 between oscillator 2 and 1 to the timeT 0 1 − ϕ new 1 that oscillator 1 still needs to evolve before spiking. Spike time t spike 3,n−1 can be obtained in a similar way, and thus one finds: The next spike of the relay oscillator occurs at t T 0 , we find the following return map for our local stability analysis Frontiers in Computational Neuroscience www.frontiersin.org

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 spike 1,n−1 = t spike 3,n−1 . 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 spike 1,l , this new relay spike arriving at time t spike 1,n−1 + 2τT 0 changes the coupling strength bỹ ε n−1 = ε n−1 + ε n−1 n−1 l=0 W t spike 1,l − t spike 1,n−1 + 2τT 0 60 < ε n−1 .
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ε n−1 and ε n .

Frontiers in Computational Neuroscience
www.frontiersin.org