Ramana Dodla and Charles J. Wilson  Neural Computation 10:2545-2610.

We investigate why electrically coupled neuronal oscillators synchronize or fail to synchronize using the theory of weakly coupled oscillators. Stability of synchrony and antisynchrony is predicted analytically and verified using numerical bifurcation diagrams. The shape of the phase response curve (PRC), the shape of the voltage time course, and the fre- quency of spiking are freely varied to map out regions of parameter spaces that hold stable solutions. We find that type 1 and type 2 PRCs can hold both synchronous and antisynchronous solutions, but the shape of the PRC and the voltage determine the extent of their stability. This is achieved by introducing a five-piecewise linear model to the PRC and a three-piecewise linear model to the voltage time course, and then analyz- ing the resultant eigenvalue equations that determine the stability of the phase-locked solutions. A single time parameter defines the skewness of the PRC, and another single time parameter defines the spike width and frequency. Our approach gives a comprehensive picture of the rela- tion of the PRC shape, voltage time course, and stability of the resultant synchronous and antisynchronous solutions.

Figure 1: Parameterizing the voltage and the PRC shapes.

(a) Model piece-wise linear voltage (thick lines) and (b) phase response curve, PRC (thick lines) shapes employed in the study. The voltage time course gives one time parameter (W/T, normalized spike width parameter) and three amplitude parameters, V_p, V_th, and V_m, which can be cast in terms of a₁, a₂, and a₃. The ratio a₃/a₂ will be used as a useful parameter in the bifurcation diagrams. The PRC, Z(t) gives just one time parameter, the PRC skewness, A, and one amplitude parameter B, which is the maximum delay. The maximum advancement C(> 0) can be used to normalize B. The thin overlaid curves are the model curves of the standard Hodgkin-Huxley equations with an applied current of 10 μA/cm₂. The model can be fit with the parameters A/T = 0.567, W/T = 0.075, B/C = −0.5, and a₃/a₂ = 0.2234. These are the only four independent parameters that the stability boundaries of synchrony and antisynchrony depend on. (c) Compari- son of the interaction function and (d) the growth function computed from the HH model (thin line) and the piecewise linear model (thick line). The piecewise model predicted the stability of the synchrony (slope of G at spike time difference 0) and antisynchrony (slope of G at spike time difference T/2) accurately.

Figure 2: Segmenting the parameter space and finding the stability.

(A) Left: (W, A) space for studying synchrony. Right: The arrangement of Z(t) and V(t) for determining the stability of synchrony. Cases a–d depict the parameter regimes that have different expressions for the eigenvalue and thus different stability criteria. The destabilizing segments of Z(t) and V(t) are gray-shaded in the plots for type 1 and type 2 PRCs in each of the four cases. The six eigenvalue segments λ_1a,...,λ_6a (and likewise λ_1b, and so on) together contribute to the total eigenvalue. (B) Left: (W, A) space for studying antisynchrony. Right: Illustration of the arrangement of Z(t) and V(t) for determining the stability of antisynchrony in 4 of 17 regimes. As in panel A, the shaded portions of the segments contribute to instability. Parameters for a few models and experiments are marked in the (W,A) planes: Hodgkin-Huxley (HH) model depicted in Figure 1, Morris-Lecar model (ML) at an applied current of 0.11 μA/cm2, Erisir et al.’s model at low (E-1) and high (E-2) frequencies (discussed in Figure 12), and Mancill et al.’s experimental recordings at low (M-1) and high (M-2) frequencies (see Figure 12).

Figure 3: When the spike width is zero, the large skewness in type1 PRCs could make the antisynchronous state stable and cause bistability between synchronous and antisynchronous states.

(a) Voltage time course from equation 2.7 at W/T = 0. (b) Z(t) from equation 2.14 when B/C = 0. Eigenvalue for the synchronous state (c) and the antisynchronous state (d) as the PRC skewness is increased. (T = 1, C = 1, a₃ = 24 mV, ε = 1.) (e) Numerical bifurcation diagram showing stable (solid lines) and unstable (dashed lines) phase-locked states with skewness. The synchronous state (phase-locked state at 0) is stable for all skewness, whereas the antisynchronous state (phase-locked state at T/2) acquires stability for large skewness. Other phase-locked states exist at large skewness but are unstable. (f) Stable and unstable regions along skewness.

Figure 4: When the spike width is zero, a large skewness in type 2 PRCs can destabilize a synchronous state, still cause bistability between synchronous and antisynchronous states, and destabilize both synchronous and antisynchronous states.

(a) Voltage time course from equation 2.7 at W/T = 0. (b) PRC profiles from equation 2.14 when B/C = −0.5. The eigenvalue and the different compo- nents that make up the eigenvalue, which determines the stability of the syn- chronous (c) and antisynchronous (d) states. The bifurcation diagram (e) shows stable (solid curves) and unstable (dashed curves) phase-locked states as the skewness is increased. Different stable regimes are marked pictorially in panel f.

Figure 5: The stability diagram and its verification in the absence of spike width.

(a) Stability regions in (A/T, B/C) plane for zero spike width. (At zero spike width, a₃/a₂ does not affect the stability boundaries.) The diagram does not include the edge effects at A = 0. Synchrony is stable for all type 1 PRCs and for type 2 PRCs if the type parameter B/C > ρ₁. The antisynchrony is mostly confined to large skewness but is also possible for very large, positive B/C with small skewness. The circled numbers mark parameter values, which are used to compute numerically the growth function [G(φ)] from the voltage and PRC profiles. (b–d) The growth functions computed from the Z(t) and V(t) profiles for the parameters marked in panel a.

Figure 6: Effect of nonzero spike width on type 1 PRC neurons illustrated when B/C = 0 and W/T = 0.15.

For this choice, antisynchrony rather than synchrony could become stable for very small skewness. (a) Voltage time course from equation 2.7. (b) Z(t) from equation 2.14 when B/C = 0. The eigenvalue (thick curves) and their components (thin curves) that determine the stability of synchrony (c) and antisynchrony (d) are plotted as the skewness is increased. λ_3a is the result of spike downstroke and causes the synchrony to become unstable. The sharp upstroke contributes to γ_4a that stabilizes the antisynchrony, but γ_5a and γ_6a, which are the result of the spike downstroke, cause the instability of the antisynchronous state as the skewness is increased. (e) One-parameter bifurcation diagram showing the stable (solid curves) and unstable (dashed curves) phase-locked solutions. (f) Stable synchrony and antisynchrony are pictorially depicted as skewness is increased. The white space holds other stable phase-locked solutions in panel e. a₃/a₂ = 0.2234.

Figure 7: Effect of nonzero spike width on type 2 PRC neurons illustrated when B/C = −0.5 and W/T = 0.15.

For this choice, antisynchrony could become stable in two ranges: at small and large skewness levels. Boundaries of synchrony are moderately sensitive to decreasing B/C. (a) Voltage time course from equation 2.7. (b) Z(t) from equation 2.14 when B/C = −0.5. Total eigenvalue (thick curve) and its components (thin curves) that determine the stability of synchrony (c) and antisynchrony (d) are shown as a function of skewness. (e) One-parameter bifurcation diagram showing stable (solid curves) and unstable (dashed curves) phase-locked solutions as a function of A/T. (f) The stability regions are shown pictorially as a function of skewness. a₃/a₂ = 0.2234.

Figure 8: At moderate skewness (here A/T = 0.2 which is less than that of the HH model) the parameter regime where the synchronous state is stable expands to larger values of W/T .

The antisynchronous state is not accessible for small W/T. The eigenvalue components (thin curves) and the total eigenvalue (thick curve) that determine the stability of the antisynchronous state as W/T is increased are shown in panel a for a type 1 PRC (at B/C = 0.2). The stability diagram in (W/T, a₃/a₂) plane is plotted in panel b. The interaction and growth functions for a parameter point that supports only other phase-locked states are shown in panel c, and a one-parameter bifurcation diagram is shown (d) at a₃/a₂ = 0.2234. Similar to panels b and d, the results obtained for a type 2 PRC (B/C = −0.25) are shown in panels e and f. In panels d and f, solid lines are stable branches, and dashed lines are unstable.

Figure 9: Large skewness (A/T=0.6 here) can make the antisynchronous state accessible for a range of W/T values starting at zero. Such a regime is thinner in type 1 than in type 2 PRCs.

(a) Stability regions in (W/T,a₃/a₂) space for B/C = 0.2, 0, −0.25. The curves ρ_a and ρ_b are the boundaries for synchrony. The curves σ_d and σ_e are the boundaries for the antisynchrony. (b) One-parameter bifurcation diagram showing stable (solid) and unstable (dashed) states as a function of W/T at a₃/a₂ = 0.22 that shows bistability of synchronous (phase- locked state at 0) and antisynchronous state (phase-locked state at T/2) at small W/T. (c) The eigenvalue components (thin curves) and the total eigenvalue (thick curve) that determine the stability of antisynchronous state as a function of W/T at a₃/a₂ = 0.22. (d) Growth functions at a₃/a₂ = 0.22 at a few values of W/T depicting how the antisynchronous state becomes unstable at large W/T.

Figure 10: Exploring the PRC parameter space

Stability regions of both synchrony and antisynchrony in the PRC skewness versus its type parameter space at W/T = 0.02 (a), 0.05 (b), 0.15 (c), and 0.3 (d) when a₃/a₂ = 0.2234. The HH model discussed in Figure 1 (with I_app = 10μA/cm2 resulting in W/T = 0.075, B/C = −0.5 and A/T = 0.567) lies in a stability diagram that is nearly identical to panel b in the bistability region above the curve ρ_c and slightly to the right of σ_d . The unmarked white space to the right of the vertical dashed line is the for- bidden parameter space by the condition in equation 2.15. But the white space regions to the left of the vertical dashed line that are interspersed between synchronous and antisynchronous states hold other nonzero phase-locked states. While a number of models display a wide range of PRC skewness levels, the type parameter (B/C) for many neuronal models is above −1.

Figure 11: Exploring the spike width and the PRC-type parameters.

The stability of the synchronous and antisynchronous states at A/T = 0.25 (a, c, e), and 0.6 (b, d, f) in the (B/C,W/T) plane. The voltage profile for panels in each column is illustrated at the top. Panels a and b depict the effect of skewness at very small spike threshold (i.e., the case of tall spike). Panels c and d are for the HH model parameter of a₃/a₂ = 0.2234 that is depicted in Figure1. Panels e and f are for a very high spike threshold (i.e., the case of a short spike). The white space holds other nonzero phase-locked states. While a synchronous state may occur for small W/T, the occurrence of an antisynchronous state at either small or large W/T depends on the level of PRC skewness.

Figure 12: Demonstrating piecewise linear formulation for Erisir et al. (1999) model (a–f) and Mancilla et al. (2007) experiments g–i.

(a–f) The voltage, PRC, and resultant growth functions computed from the original models are shown as thin lines, and the corresponding computations with PWL approximation are shown in thick lines. (g–i) Mancilla et al.’s neocortical recordings for the voltage and the PRC and the resultant growth functions are shown at two different frequencies: open circles at low frequency and filled circles at high frequency. (Data were kindly provided by Jaime G. Mancilla, and the data displayed are the average of multiple trials.) The PWL approximation for the experimental voltage and PRC traces and the resultant growth functions predicted stability of synchrony and instability of antisynchrony agreeing with those computed experimentally.