Effects of conduction delays on the existence and stability of one to one phase locking between two pulse-coupled oscillators

Gamma oscillations can synchronize with near zero phase lag over multiple cortical regions and between hemispheres, and between two distal sites in hippocampal slices. How synchronization can take place over long distances in a stable manner is considered an open question. The phase resetting curve (PRC) keeps track of how much an input advances or delays the next spike, depending upon where in the cycle it is received. We use PRCs under the assumption of pulsatile coupling to derive existence and stability criteria for 1:1 phase-locking that arises via bidirectional pulse coupling of two limit cycle oscillators with a conduction delay of any duration for any 1:1 firing pattern. The coupling can be strong as long as the effect of one input dissipates before the next input is received. We show the form that the generic synchronous and anti-phase solutions take in a system of two identical, identically pulse-coupled oscillators with identical delays. The stability criterion has a simple form that depends only on the slopes of the PRCs at the phases at which inputs are received and on the number of cycles required to complete the delayed feedback loop. The number of cycles required to complete the delayed feedback loop depends upon both the value of the delay and the firing pattern. We successfully tested the predictions of our methods on networks of model neurons. The criteria can easily be extended to include the effect of an input on the cycle after the one in which it is received.

Thanks to Lakshmi Chandrasekaran for a careful reading of an earlier draft and to Shuoguo Wang for assistance with the figures. This work was supported by NIH grant 5R01MH085387-03 to CCC under the CRCNS Collaborative Research in Computational Neuroscience program.

Authors and Affiliations

1. Neuroscience Center of Excellence, Louisiana State University Health Sciences Center, 2020 Gravier St Suite D, New Orleans, LA, 70112, USA

   Michael Marmaduke Woodman & Carmen C. Canavier

2. Université de la Méditerranée, Institut des Sciences du Mouvement, UMR6233 CNRS, 163 Avenue de Luminy, 13288, Marseille, France

   Michael Marmaduke Woodman

Corresponding author

Correspondence to Carmen C. Canavier.

Action Editor: Bard Ermentrout

Appendix 1

The presence of second order resetting, or an effect on the second spike after the perturbation, implies that the trajectory in the state space has not returned to the limit cycle by the time of the first spike after the perturbation. To correct for the effect on the timing of the second spike, the phase immediately after the first spike is not presumed to be reset to zero, but rather to a value that will cause the second spike to occur at the proper time. Therefore the second order resetting must be included in the following stimulus interval (Achuthan and Canavier 2009; Oprisan et al. 2004; Sieling et al. 2009), and be complete by the end of this interval. Third and higher order resetting cannot be accommodated in this formalism and must be negligible in order for the assumptions and methodology to be valid.

The stimulus interval in the closed loop is then defined based on the open loop stimulus interval, but modified to account for the fact that neurons in the closed circuit are perturbed multiple times as opposed to a single time in the open loop circuit. Therefore, \( t{s_i}\left[ n \right] = {P_i}({\phi_i}\left[ n \right]) + {f_{2i}}({\phi_i}\left[ {n - 1} \right])) \), where ϕ _i [n] is the phase at which the input is received by neuron i in cycle n and the second order resetting (indicated by the 2 in the first subscript) term accounts for the possibility that the trajectory has not fully recovered to the limit cycle at the time of the spike initiating the stimulus interval.

In order to derive stability results from Eqs. (7) and (8), the PRC for both first and second ( jth) order resetting is linearized as

Equations (9 and 10) with second order resetting become

where m _ji is the slope of the jth order PRC f _ji '(ϕ _i [∞]) evaluated at the locking point for the ith neuron. Again, the summation terms are only valid for k > 1.

The addition of second order resetting increases the order for k = 1 from 1 to 2, and adds two for all other k for an order of k + 1. When the linear system constructed from the equation above is evaluated with k = 1, j ₁   = 0 and j ₂   = 0 (Scheme A in Fig. 2), a second order linear system is obtained, and the eigenvalues are the roots of the following quadratic polynomial:

which agrees with results for this firing pattern with no delay (Oprisan et al. 2004).

On the other hand, when the linear system is evaluated at k = 2, j ₁  = 0 and j ₂  = 1 (scheme B in Fig. 2), the three eigenvalues that determine stability are the cubic roots of the following polynomial:

The introduction of second order resetting complicates the analysis. However, it gives the advantage that we no longer require that the trajectory return to the limit cycle between the time a neuron receives an input and when it fires next (the end of the next recovery interval) but rather by the time it next receives an input (the end of the next stimulus interval). Therefore although the “observed phase” when a neuron spikes is one, or equivalently zero, if the trajectory has not yet returned to the limit cycle, the “true” phase is determined by the isochron on which the trajectory lies. The correct phase at the time of the first spike observed after an input can be inferred by observing the shift in the timing of the second spike after the input, and even negative phases which are not on the original limit cycle can be accommodated by this method (see also Oh and Matveev 2009). Second order resetting can be important especially when inputs are received at late phases and bleed over into the second cycle.

Second order resetting is negligible for the Morris Lecar model used in Figs. 5 and 8, but this is not always true. Here we give an example in which second order resetting cannot be neglected. We used two identical mutually inhibitory Wang and Buzsaki model neurons with parameters given in Appendix 2 and an intrinsic period of 16.75 ms. The predicted modes at a delay of 16 ms are indicated by symbols on both the first (solid lines Fig. 9) and second order PRC (dashed lines). A stable antiphase mode (open squares) with a time lag of 9.48 and a stable synchronous mode (open circles) with a network period of 16.77 ms were predicted, and each was observed exactly at the predicted values when the system of differential equations was initialized with their respective basins of attraction. For the antiphase mode, second order resetting was negligible. However, the slope of the second order resetting at a phase of 0.962 corresponding to synchrony was 0.765 compared to the -0.52 for the slope of the first order resetting, so second order resetting is clearly not negligible, and in fact exerted a powerful stabilizing influence on synchrony. The k value for this mode was 1, so the full polynomial expression was used to find the eigenvalues in this case: \( {\lambda^2} - \left( {\left( {{m_{11}} - 1} \right)\left( {{m_{12}} - 1} \right) - {m_{21}} - {m_{22}}} \right)\lambda + {m_{21}}{m_{22}} = 0 \) .

In this case, neglecting second order resetting produces incorrect stability results for synchrony. Stability is only predicted if the second order term is considered. Note that the slope of the second order resetting is steeper than that of the first order resetting at synchrony (open circles). This example is for two Wang and Buzsaki neurons coupled by inhibition with a conduction delay of 16 ms and parameters as given in Appendix 2

Full size image

When second order resetting was included, all eigenvalues had absolute value less than one: a real eigenvalue of -0.76 and a complex conjugate pair with an absolute value of 0.87. On the other hand, ignoring second order resetting moves the locking point to a phase of 0.955, changes the predicted network period to 16.96 ms and results in an eigenvalue of -2.6 that incorrectly predicts that synchrony is not stable at a delay of 16 ms.

Appendix 2

1.1 Model equations and parameters

1.1.1 Morris-Lecar model network

In Fig. 5, we simulated a model network composed of two Morris and Lecar (1981) oscillators with Type II excitability (Rinzel and Ermentrout 1998), reciprocally coupled by identical, excitatory synapses. The differential equations for each oscillator are

where

C _m and V are the membrane capacitance and potential, and t is time. Constants g _Ca , g _K , g _L and g _syn and E _Ca , E _K , E _L and E _syn are the maximal conductances and reversal potentials for the calcium, potassium, leak and synaptic currents. Potassium activation is denoted by w, and m is the steady-state activation of the calcium current; the bar indicates the steady-state values as a function of potential. φ is a time-scale factor.

Post-synaptic activation is modeled by the equation: \( ds/dt = \alpha [1/(1 + \exp ( - {V_{pre}}(t - \delta )/2))](1 - s) - s/{\tau_{syn}}, \)

where α is the presynaptic activation rate, τ _w is the synaptic decay constant, and V _pre (t-δ) is the delayed presynaptic voltage.

C _m was set to 20 μF/cm²; I _stim , 100 μA/cm² ; g _Ca , 4.4 mS/cm², g _K , 8 mS/cm²; g _L , 2 mS/cm²; E _Ca , 120 mV; E _K , -84 mV; E _L , -60 mV; E _syn , 0 mV; φ, 0.04; α, 6.25 ms^-1; τ _syn , 1 ms ; m _half , -1.2 mV; m _slope , 18 mV; w _half , 2 mV; w _slope , 30 mV.

1.1.2 Wang and Buzsaki model network

In Fig. 7b, we simulated a model network composed of two Wang and Buzsaki (1996) model neurons, reciprocally coupled by identical, inhibitory synapses. The differential equations for each oscillator are:

where the capacitance C = 1 μF/cm², V is the cell membrane voltage in millivolts and t is time in milliseconds. The leak current is given by Ι _L   = g _L (V–E _L ). The sodium current is given by Ι _Νa   = g _Νa m _∞ ³ h(V-E _Na ). The steady state activation is m _∞   = α _m /(α _m   + β _m ) where \( {a_m}(V) = --0.1\left( {V + 35} \right)/\{ exp[--0.1\left( {V + 35} \right)]--1\} \) and \( {\beta_m}(V) = 4exp[--\left( {V + 60} \right)/18] \).

The scale factor φ was set to 5. (The symbol was changed from ϕ to φ to avoid confusion with the symbol for phase). The rate constants for the inactivation variable h are given by \( {a_h}(V) = 0.07exp[--\left( {V + 58} \right)/20] \) and \( {\beta_h}(V) = { }1/\{ exp[--0.1\left( {V + 28} \right)] + 1\} \). The potassium current is given by \( {I_K} = {g_K}{n^4}(V--{E_K}) \). The rate constants for n are \( {a_n}(V) = --0.01\left( {V + 34} \right)/\{ exp[--0.1\left( {V + 34} \right)]--1\} \) and \( {\beta_n}(V) = 0.125\,exp[--\left( {V + 44} \right)/80] \). The reversal potentials E _Na , E _K and E _L were set to 55, –90 and –65 mV, respectively. The maximal sodium (g_Na), potassium (g_K) and leak (g_L) conductances were set to 35, 9 and 0.1 mS/cm², respectively. I_stim,i is the applied current and was set at 1.0 μA/cm² . The synaptic current is given by Ι _syn   = g _syn s(V–E _syn ), where g _syn is the maximum synaptic conductance equal to 0.15 mS/cm² and E _syn is equal to -75 mV. Synaptic activation was modeled as described above for the Morris Lecar network.

Appendix 3

1.1 Example interval maps

1.1.1 Map for k = 1 (scheme A)

From Eq. (7 and 8) in main text:

Dependence of recovery interval on stimulus interval is given by the PRC under the assumptions of pulsatile coupling:

The above step can only be performed if second order resetting is ignored; otherwise the map must be formulated in terms of phase and not time intervals.

Substituting for the recovery intervals, then eliminating ts ₁ [n], we obtain the following nonlinear map.

The above nonlinear map contains all information required to extract the stimulus and recovery intervals for both neurons on each cycle as long as the firing pattern remains constant. Linearizing the above expression about the fixed point by setting \( t{s_2}\left[ n \right] = t{s_2}[\infty ] + \Delta t{s_2}\left[ n \right] \), linearizing the PRC in the neighborhood of the fixed point as described in the main text, and cancelling steady terms we obtain the 1D linear map:

where \( P_{1} - P_{2} + ts_{2} [n - 1] + P_{2} f_{2} (ts_{2} [n - 1]/P_{2} ) + \delta _{1} + \delta _{2} )/P_{1} = \phi _{1} [\infty]\) and \( t{s_2}[\infty ]/{P_2} = {\phi_2}[\infty ] \).

This is exactly equivalent to the eigenvalue λ _k=1   = (1–m ₁ )(1–m ₂ )derived in the main text.

1.1.2 Map for k = 2 (scheme B)

From Eqs. (7) and (8) in main text:

Additional Constraint

Substituting the additional constraint in the Equation for ts ₂ [n + 1] above results in the equality:

Making the appropriate substitutions for tr ₁ [n] then ts ₁ [n] into the equation for ts ₂ [n + 1] we obtain:

The above nonlinear map again contains all information required to extract the stimulus and recovery intervals for both neurons on each cycle as long as the firing pattern remains constant. As before, we obtain the 1D linear map:

where \( (2{\delta_1} - t{s_2}[\infty ])/{P_1} = {\phi_1}[\infty ] \) and \( ts_{2} {\left[ \infty \right]}/{\text{P}}_{2} = \phi _{2} {\left[ \infty \right]} \) .

This is exactly equivalent to the eigenvalue \( {\lambda_{k = 2}} = 1 - {m_1} - {m_2} \) derived in the main text.

Reprints and permissions