Effect of Sharp Jumps at the Edges of Phase Response Curves on Synchronization of Electrically Coupled Neuronal Oscillators

We study synchronization phenomenon of coupled neuronal oscillators using the theory of weakly coupled oscillators. The role of sudden jumps in the phase response curve profiles found in some experimental recordings and models on the ability of coupled neurons to exhibit synchronous and antisynchronous behavior is investigated, when the coupling between the neurons is electrical. The level of jumps in the phase response curve at either end, spike width and frequency of voltage time course of the coupled neurons are parameterized using piecewise linear functional forms, and the conditions for stable synchrony and stable antisynchrony in terms of those parameters are computed analytically. The role of the peak position of the phase response curve on phase-locking is also investigated.


Introduction
A phase response curve (PRC) quantifies temporal deviations of an oscillator in response to an oncoming stimulus [1][2][3]. Methods based on PRCs have been extensively used to predict when synchronization could occur between biological oscillators [4][5][6][7]. Studies based on weakly coupled oscillator theory applied to coupled neurons have often sought to relate the shape of the PRC to the emergence and the stability of the phase-locked states [8][9][10][11][12][13][14]. However, some PRCs measured experimentally exhibit significant departures from being smooth and continuous [15][16][17][18]. In particular these PRCs may show sudden increments in their level at the beginning and/or end of the oscillation cycle. There are also some neuronal models that display such behavior. Such sudden increments in the PRC level are often due to fast gating dynamics of the neurons, but could also be due to the effect of higher order PRCs that may become significant in the collective dynamics of coupled neurons. If the sudden rise in the PRC level is indeed due to fast gating dynamics, then, for simplicity, it could rather be approximated by a discontinuous jump in the PRC. When higher order PRCs become significant, then their shapes may also have to be incorporated in the relevant analysis.
Sharp jumps in the PRCs were earlier shown to be caused by fast potassium gating dynamics [19], presence of adaptation currents such as calcium-dependent afterhyperpolarization (AHP) current or muscarinic voltage-dependent potassium (M) current [20,21], or abrupt dynamical changes in the modeling equations. AHP current, for example, can cause the neuron become less sensitive at the early phases and thus impart skewness to the PRC. In some models it can also impart sudden jumps at early phases [21].
The PRCs of a leaky integrate-and-fire model [22] and quadratic integrate-and-fire model [23,24] display discontinuities at the beginning and the end of the oscillation cycle. Adapted exponential integrate-and-fire neuron model [21,25] is another example that displays sharp PRC jumps. In all these cases weakly coupled oscillator theory has been used to predict the stability of synchrony and antisynchrony. As we will also illustrate, adaptation is not necessary to realize PRC with sharp jumps. But even when adaptation was present, we expect that the theory would still be applicable [20] because the effect of adaptation on synchrony is via a modification of the shape of the PRC, and/or the voltage time course; In the synchronized state the change of frequency caused by adaptation can be assumed to be negligible.
Here we address comprehensively the role of the discontinuous jumps at the ends of the PRCs in the synchronizability of coupled oscillatory neurons when the coupling between them is weak and electrical. Only the first order PRC is used in the analysis. To model the phase response curve, we employ a piecewise linear approach that allows a detailed study of the dependence of the PRC shape on synchrony, while at the same time being applicable to experimentally determined PRCs. The PRC profile is constructed with only two piecewise linear segments, and the voltage profile with three piecewise linear segments. We predict when synchrony and antisynchrony become stable as the level of the discontinuous PRC jumps, and the spike width and frequency are varied. Our study complements other similar studies on electrically coupled neuronal networks that used leaky integrateand-fire models [12,26], and generalizes the results of those that used quadratic integrate-and-fire models [23,24] by considering a range of PRC shapes and voltage time courses. We also study how the location of the PRC maximum (the skewness) affects synchrony and antisynchrony. A network simulation of Wang-Buzsáki model neurons when each neuron displays a PRC with sudden rise in its level near zero-phase is also presented.

Model
The leaky integrate-and-fire, quadratic integrate-and-fire, and adaptive exponential integrate-and-fire models whose PRCs are depicted in Fig. 1(a-c) are described in the figure caption. The modified Wang-Buzsáki model [19] whose PRC is depicted in Fig. 1(d) is given by the following evolution equations: The parameters are C m~1 mF/cm 2 , E L~{ 65 mV, E Na~5 5 mV, E K~{ 90 mV, G L~0 :1 mS/cm 2 , G Na~3 5 mS/cm 2 , and G K~9 mS/cm 2 . The applied current I app when set above 0:16 mA/cm 2 triggers spontaneous oscillations. Ignoring a very brief downward and negative swing near zero phase, the PRC of the model may be considered a type-1 and may be treated like a PRC with sharp jump at the left edge [ Fig. 1(d)]. The only difference between our formulation and the original formulation of this model is in making the time-scale factor for sodium inactivation and potassium activation independent: w h and w n that are now used to control the PRC shape. This model is also used later in network simulations where the current due to electrical coupling I ele is proportional to the difference of the voltage of the coupled neurons.
The main results of the paper use the formulation laid out below. Extension of this model that incorporates PRC skewness is presented in the later part of the Results section. We formulate the PRC and the voltage time course using piecewise linear (PWL) functions, and then present a method to find the stability of synchrony and antisynchrony. The choice of PWL functions facilitates analytical determination of stability boundaries. The phase response curve is formulated as a function with two piecewise linear profiles [Y 1 (t) and Y 2 (t)] that exhibit finite jumps or discontinuities at the edges (i.e. at t~0 and at t~T, the period of oscillation of the neuron). The assumption of finite jumps at the edges is not necessarily due to a discontinuity in the PRC, but sharp rise or fall at those phases. But for computation of stability of phase-locked solutions, the assumption of a discontinuous jump at the edges does not lead to any artifacts unless another parameter such as spike width also simultaneously becomes zero; In such a case the effect of discontinuity must be explicitly incorporated into the analysis, or the analysis must be carried out in the limit of those parameters going to zero. The PRC is formulated as below: where B is the amount of jump on the left edge (t~0), B 2 is the amount of jump on the right edge (t~T) of the PRC, and C (.0) is the maximum advancement of the PRC. When B~B 2~0 , the PRC becomes symmetric. The PRC profile is depicted in Fig. 2(a) for a few parameters of B and B 2 . A monotonically increasing PRC as in a leaky integrate-and-fire model is obtained by setting C~(BzB 2 )=2. Such PRCs are also treated in the Results section as a special case. No assumption is made on the sign of B and B 2 in deriving the stability regions. The voltage profile is formulated by the following three piecewise linear curves [ Fig. 2(b)] that are modeled after an empirical observation of the voltage profile of the Hodgkin-Huxley (HH) model equations; Spike width parameter W is the only time varying parameter that appears in the model, and the spike amplitude is controlled by spike peak V p , spike threshold V th , and spike minimum, V m : if T{ W 2 ƒtvT: At an external applied current of 10 mA/cm 2 , the HH voltage time course can be approximated with V m~{ 72 mV, V p~3 5:43 mV, V th~{ 48 mV, and W~1:1 ms, T~14:636 ms (i.e. an oscillation frequency of 68.3 Hz). Thus W =T~0:075 for this model. We define a 2~Vp {V m , and a 3~Vth {V m . We term W the spike width for simplicity although the actual spike width could be up to 5W =2. The spike width and threshold to spike height ratio (a 3 =a 2 ) are freely altered to explore the stability boundaries in these parameter spaces. But we assume that Because the spike downstroke and spike upstroke together add up to a width of 5W =2, we insist that T is bigger than that. That is, Pairs of identical nonlinear oscillatory neurons that may be originally described by several state variables, but oscillate with a constant period T and possess a voltage profile V (t), and a PRC, Y (t), when coupled electrically by a coupling function proportional to the difference of their individual voltages can be reduced to pairs of phase evolution equations under the assumption of weak coupling [1][2][3][27][28][29][30][31]. Two such identical neurons can be reduced to two phase evolution equations with phases w 1 and w 2 (that range from 0 and T) described by the following two equations where (w0) is the strength of coupling. H(w) is the interaction function and quantifies the instantaneous increase of an oscillator frequency due to coupling to the other oscillator, and is expressed as follows: Phase-locking occurs when the phase difference w~w 2 {w 1 remains constant in time, and hence it is convenient to study the equations in terms of the phase difference that evolves according to where represents the growth function that quantifies the instantaneous growth of the phase difference. The change of the growth function as a function of the phase difference [G 0 (w)] at a stable phaselocked state (w~w Ã ) should become negative so that any small perturbation to that steady state subsides. Thus the stability of the synchronous state (w Ã~0 ) of the Eq. 5 is determined by the eigenvalue:  [23]: t _ V V~I app zV 2 where V is reset to V r when it crosses a threshold level V th . (c) PRC of adaptive exponential integrate-and-fire model [21,25]: {w such that V and w are reset, respectively, to V r and wzb whenever V reaches a peak level V Ã . The parameters are C m~0 :1 nF, G L~0 :01 mS, V r~{ 60 mV, V Ã~{ 20 mV, E L~{ 70 mV, and t w~1 00 ms. The level of adaptation is controlled by a. (d) PRC of Wang-Buzsáki model that is described in the Model section. Models in (a, b, d) have no adaptation and their second and higher order PRCs are identical to that of the first order. doi:10.1371/journal.pone.0058922.g001 and the stability of the antisynchronous state (w Ã~T =2) of the Eq. 5 is determined by the eigenvalue: These eigenvalues are computed for the Y (t) and V(t) profiles formulated in Eqs. 1 and 2, and the stability of synchrony and antisynchrony is determined in the following sections. We will see that the PRC jumps (B and B 2 ) can alter the stability of both synchrony and antisynchrony. A symmetric PRC [ Fig. 2(a), thin curve] leads to stable synchrony and unstable antisynchrony [ Fig. 2(c), thin curve]. An example of the synchronous state becoming unstable and the emergence of a non-zero phase-locked state that is very close to the synchronous state due to the finite jumps at the PRC ends is shown in Fig. 2(a,c). Antisynchrony and other non-zero phase-locked states can also undergo stability changes. The case of zero spike width (W =T~0) presents an enigmatic situation when the unstable synchronous branch and a stable non-synchronous but phase-locked branch converge at the same point [ Fig. 2(d)]. This situation arises because of the fact that the discontinuity in the voltage due to W~0 and the discontinuous jump in the PRC occur at the same temporal location. The stable phase-locked branch is very close to the unstable synchrony branch, and they quickly get separated as W =T is increased. Hence strictly at W =T~0, the stability criterion for phaselocked states depends on whether it is computed in the limit of W =T?0 or otherwise. We will compute the stability of synchrony and anti-synchrony in the limit of W ?0 that inherently uses the spike effect, and will see a transition from stable synchrony to unstable synchrony above a critical B=C. For non-zero W =T, this situation does not arise because the edge effects get factored into the eigenvalue components computed due to the up and downstrokes that span finite time widths. The antisynchronous state is unstable for W =T~0 unless the sum of the PRC jumps is bigger than twice the PRC maximum (i.e. BzB 2 w2C). But as W =T increases, more parameter region is filled with antisynchrony, some with synchrony, and some with bistability. The bistability could occur between phase-locked states, synchronous, and antisynchronous states.
The stability boundaries do not depend on the time period T, but in numerical simulations we used the period corresponding to

Results
Except the leaky integrate-and-fire (LIF) model which has an exponential form of PRC, the other three PRCs [quadratic integrate-and-fire (QIF), adaptive exponential integrate-and-fire (aEIF), and modified Wang-Buzsáki model] presented in Fig. 1 can be approximated by PWL formulations in our models. The results will also be applicable to the LIF model in a qualitative manner. The LIF and QIF models both have sharp jumps in the PRC at both early and late phases, and no adaptation currents are present in either model. The aEIF model has an adaptation current variable and shows a left PRC jump. The modified Wang-Buzsáki model that has only sodium and delayed rectifier currents, and no adaptation currents, displays a sudden downward swing followed by a sudden rise in the PRC level. We first study symmetric PRCs depicted in Fig. 2(a). Rigorous analytical arguments are presented for the boundaries of both stable synchrony and stable antisynchrony (Figs. 3 and 4) when the PRC is symmetric (Eq. 1, Fig. 2). Later we introduce PRC skewness that parametrizes the position of the peak PRC level, and analytical results are discussed ( Fig. 5) when the PRC in addition acquires a skewness (Eq. 16) such that the peak of the PRC is either moved to the left or right. Finally we present simulation results using networks of Wang-Buzsáki model.

Synchrony and Antisynchrony when Spike Width is Zero
The boundaries of synchrony (r 1 and r 2 ) and antisynchrony (s 1 ) under the assumption of zero spike width are derived in this section, and are illustrated in Fig. 3(a,b) in the parameter spaces of (B 2 =C, B=C) and (B 2 =C, a 3 =a 2 ). In the absence of any jumps in the PRC, the coupled system synchronizes because the entirely positive PRC encounters a positively sloped voltage segment and thus the eigenvalue [Eq. 6] becomes negative. In other words, if the second neuron leads the first by a small phase, the first neuron speeds up in response because the voltage time course of the second neuron (say, v 2 ) is higher than the first (say, v 1 ) for most of the spike interval (coupling to the first neuron is proportional to v 2 {v 1 ). And if the second neuron lags the first, the converse effect occurs.
A positive jump in the PRC at zero phase (left side edge) helps destabilize synchrony, and a positive jump at T (right side edge) helps stabilize synchrony more. We can visualize these effects by imagining the case of non-zero spike width. The spike downstroke has large negative slope, and the spike upstroke has large positive slope. Since the PRC is non-zero, the convolution (Eq. 6) with the downstroke results in a positive eigenvalue that helps destabilize synchrony, and that with the upstroke results in a negative eigenvalue that helps stabilize synchrony. In other words, if v 2 leads v 1 , the first neuron slows down in the region from spike peak to (nearly) the spike downstroke because the spike downstroke of v 2 falls below that of v 1 , whereas it speeds up during the spike upstroke region because there v 2 is bigger than v 1 . In the limit of zero spike width, these effects remain because the PRC is non-zero due to finite jumps. Thus the destabilizing effect of the left edge (B) can be countered by appropriately increasing the right side jump (B 2 ) leading to a diagonal line of criticality in the parameter space of (B 2 ,B 1 ) [ Fig. 3(a)]. Increasing the spike threshold such that a 3 =a 2 is also increased will cause the voltage time course acquire more positive slope that enhances the negative eigenvalue component, and consequently we see the boundary of synchronous region being pushed into the earlier unstable region [ Fig. 3 or making a transition to synchrony [ Fig. 3 The instability of synchrony occurring for W~0 is more subtle than that occurring for non-zero spike width. Imagine again the case of non-zero spike width. When v 2 leads v 1 , the slow down of the first neuron occurs within the duration of spike downstroke, and then it speeds up for the rest of its cycle. When the spike width is zero, the slow down regime is really confined to zero width, but it still exists because of the non-zero value of the PRC at the edge. Thus the turn around from speed up to slow down occurs right at the zero phase, thus creating an equilibrium point (w Ã~0z ) that becomes stable when the synchrony becomes unstable. The growth function displaying a stable non-zero equilibrium (due to the negative slope) merging with unstable synchronous state (due to positive slope) as the spike width becomes zero is illustrated in Fig. 3(c). Similar phenomenon occurs even in the stable synchrony regime, leading to a stable synchrony coexisting with unstable equilibrium point (w Ã~0z ) [ Fig. 3(d)]. Thus the boundaries of synchrony derived below are in the limit of W going to zero. And for very small non-zero W these boundaries begin to change slightly in the parameter spaces, and the non-zero equilibrium and the synchronous states become more separated.
If the PRC is maximum at half period, the time-shifted (by half period) voltage not only has a large negative discontinuity, but also encounters a large PRC level at the discontinuity. Together this contributes to large positive integral in the eigenvalue of the antisynchrony (Eq. 7). In other words, if v 2 leads v 1 , v 2 {v 1 that drives the first neuron becomes large negative at half period, and thus the first neuron slows further. If this effect is not countered by the other segments of the voltage (which occurs when condition in Eq. 10 is satisfied), the antisynchrony remains unstable.
The eigenvalue for the synchronous state (w Ã~0 ) is obtained by using the formula in Eq. 6, and then adding the edge contributions from both the jumps at the end. Using the formula in Eq. 6 we get Adding all these four eigenvalue components results in the total eigenvalue that determines the stability of the synchronous state: We directly see that a positive left side jump makes the eigenvalue more positive and hence making it less likely to synchronize, and a positive right side jump makes the eigenvalue more negative and hence making it more likely to synchronize. The maximum PRC advancement works in favor of synchrony. The synchronous state is stable when the eigenvalue is negative, i.e. when where a 0 3~a 3 =a 2 , B 0~B =C, and B 0 2~B 2 =C. The region of synchrony bounded by the curve r 1 is illustrated in Fig. 3(a) for three levels of a 3 =a 2 . The above critical condition can also be written in terms of the ratio a 3 =a 2 that gives more convenient way to visualize the stability region as a function of the level of spike threshold: The stable synchrony region bounded by r 2 is illustrated in Fig. 3(b) as a function of the jump at the right edge at three levels of the left edge jumps. Next we find the stability conditions for antisynchrony. Since V (t{T=2) is discontinuous at t~T=2, the discontinuity computed as before results in an eigenvalue component that is equal to For most PRCs that may have jumps, this condition may be difficult to satisfy because it requires large drops at the edges that cumulatively exceed the maximum PRC advancement found at t~T=2. The region falls outside the depicted range of B=C in Figure 3. Edge effects when spike width is zero. (a) Stable synchrony (shaded region) and the unstable synchrony (white region) for a 3 =a 2~0 :2234 in the plane of B 2 =C and B=C. Boundary curves for two other levels of a 3 =a 2 are also shown. (b) Same as in (a) but in the plane of B 2 =C and a 3 =a 2 for B=C~0. Boundary curves for two other level of B=C are also displayed. These curves are obtained by inverting equation r 1 for a 3 =a 2 . Antisynchrony is unstable in the displayed parameter ranges in (a) and (b). (c) Growth function G(w), in the limit of zero spike width, displaying unstable synchrony but a stable phase-locked state that is very close to the synchronous state for a parameter value that is in the unstable synchrony region. (d) One-parameter bifurcation diagram as a function of B 2 =C. doi:10.1371/journal.pone.0058922.g003 Fig. 3(a), and thus the depicted range holds an unstable antisynchronous state. The stability is independent of a 3 =a 2 .

Non-zero Spike Width and Effect of Frequency
Boundaries of synchrony (r 3 and r 4 ) and antisynchrony (s 2 and s 3 ) in the presence of non-zero spike width are derived in this section. For non-zero W , the non-zero equilibrium and the synchronous state that were found merging get separated as seen in Fig. 3(c). The slopes of the spike downstroke and upstroke are proportional to the spike width, and hence their contributions to the growth of the phase difference also is proportional to W for small W . Consequently, the resultant boundaries for small nonzero W are near those obtained for zero spike width. As in the case of zero spike width, the PRC jump on the left diminishes the chances of synchrony, and that on the right promotes the chances of synchrony. Correspondingly, large PRC jump on the right and small jump on the left signifies stable synchronous region, and large jump on the left side and small jump on the right destabilizes the synchronous state [ Fig. 4(a)].
The antisynchronous state was unstable in most of the parameter space at W =T~0 due to the destabilizing effect of the voltage discontinuity of V (t{T=2) at t~T=2. With increasing W =T, the upstroke broadens at a lower rate than the downstroke, and consequently the stabilizing effect of the upstroke dominates resulting in boundary of the stable antisynchrony becoming sensitive for W =T. At large W =T, most of the parameter space is filled with antisynchronous state [ Fig. 4(b)]. Unlike the case of zero spike width, the stable near-zero phaselocked state and the synchronous state are well separated [ Fig. 4(c)], and a bistability between a non-zero phase-locked state and the antisynchronous state is found for large left PRC jumps. Increasing the frequency causes bistability between synchronous and antisynchronous states [ Fig. 4(d)], before the synchrony loses stability [ Fig. 4(e)].
Synchrony. As the spike width becomes bigger than zero, the edge effects that we had to include earlier (when W~0) are naturally contained in the contributions of the up and downstrokes of the spike. Apart from these two components, the regime from the spike minimum to the PRC maximum, and PRC maximum to the spike threshold provide the other two components to the integral in Eq. 6. But if the spike minimum occurs after T=2, then the downstroke contribution extends all the way up to T=2, and the other three components are contained in the regime twT=2.
When the spike width is small such that 0ƒW v T 4 , we obtain the region of stable synchrony as (see Methods) where , and The parameters T and a 2 can be used to further normalize W and a 3 respectively. The region of stable synchrony bounded by r 3 is illustrated in Fig. 4(a) for W =T~0:05 and a 3 =a 2~0 :2234.
And when the spike width is large (or the frequency is high) such that the spike minimum occurs after the PRC peak ( the region of stability of synchrony is given by (see Methods).
The parameters W and a 3 can be normalized, respectively with T and a 2 . The stable synchrony region bounded by r 4 is illustrated in Fig. 4(b) for W =T~0:3 and a 3 =a 2~0 :2234.
Antisynchrony. The effect due to discontinuity of V (t{T=2) at t~T=2 that existed for zero spike width is now contained in the V 3 and V 1 segments that fall on either side of t~T=2. Along with these two segments, the segment preceding the spike peak and the segment that follows the spike threshold contribute to the integral in Eq. 7. When the spike width is small such that 0ƒW v T 4 , we obtain the region of stability for antisynchrony as (see Methods) where The stable antisynchronous region bounded by s 2 is illustrated in Fig. 4(a) at W =T~0:05 and a 3 =a 2~0 :2234 (that corresponds to the HH model).
When the spike width is large such that T 4 ƒW v 2T 5 , the spike downstroke effect lasts from t~T=2 to t~T, and the other three segments fall in the first half of the period. The corresponding boundary for stable antisynchrony is more complex than before, and is given in two regimes of spike width (see Methods). The stable antisynchrony exists when When W wx 19 , we arrive at the region of stable antisynchrony as At a 3 =a 2~0 :2234 the limits for the curves s 2 and s 3 are given respectively by 1=4ƒW =Tv0:365 and 0:365ƒW =Tv0:4. At W =T~0:3 only s 2 exists, and the region bounded by it is illustrated in Fig. 4(b).

Effect of PRC Skewness
The location of the PRC peak could indeed affect the onset of synchrony and antisynchrony. We introduce the skewness by redefining the two PRC segments of Eq. 1 as below.
where the parameter A is the skewness parameter and could range from {T when the maximum PRC advancement occurs at zerophase to T when the same occurs at phase T. The shapes of the PRCs are illustrated in two special cases when only the left side jump is present [ Fig. 5(a) left top] or when only the right side jump is present [ Fig. 5(b) left top]. The other panels in the figure present corresponding transitions of stability as the level of these jumps and the skewness are altered. The one-parameter bifurcation diagrams are computed numerically, and the stability boundaries presented are derived from the expressions obtained by solving the eigenvalue equation as was done in the previous sections. For PRCs that have jumps only on the left side an example of transition of stability between phase-locked states is shown in Fig. 5(a, left middle) for small spike width. In this illustration the level of the jump is set at half of the peak PRC level. Synchrony is unstable for most of the skewness, except in a small window at very large skewness. For negative skewness the level of Z 1 during the spike downstroke is bigger than that at zero skewness, and thus the positive eigenvalue component contributed by the downstroke increases resulting in a loss of stable synchrony. However, at very large skewness the destabilizing effect of the downstroke can be countered by the enormously negative eigenvalue contributed by the spike upstroke in association with the increasing level of Z 2 segment, thus giving a thin regime of stable synchrony [ Fig. 5 ]. An example of transition of stability is shown for negative B=C in Fig. 5(a, left bottom). There is no sudden transition as B=C goes from positive to negative regime. The stability curves are continuous. For negative B=C a portion of the spike downstroke contributes to stable synchrony due to negative Z 1 segment, and thus stable synchrony can be sustained even for large negative skewness. In the previous section we saw that it is possible to make the antisynchronous state stable by increasing the spike width or frequency. This advantage is now countered by a positive skewness; Now the downstroke that is responsible for instability dominates because the entire (half-period shifted) spike profile encounters an increasing segment Z 1 . Consequently the PRC level during the downstroke is higher than that during upstroke. Thus it is difficult to attain stable antisynchrony for large skewness [ Figs. 5(a, right panels)].
The corresponding results for a right side jump of the PRC are illustrated in Fig. 5(b). For synchronous and antisynchronous states, the parameter B 2 =C behaves qualitatively like {B=C, i.e. the right side jump is qualitatively equivalent to a left side jump with the opposite sign, and vice versa [panels in Fig. 5(b)]. This equivalency seems better at low frequencies or spike widths, i.e. when W =T is very small. We can indeed determine when this equivalency is a valid approximation by examining the eigenvalues Effect of PRC Discontinuities on Synchrony PLOS ONE | www.plosone.org that define the stability of synchrony and antisynchrony. For ease of analysis, we examine the case when the skewness is absent (A~0). The eigenvalue under study is l 1 (see Methods), and for showing the exact equivalency we must show that the dependence on B (i.e. the coefficient of B) is the same as the dependence on {B 2 (i.e. the negative of the coefficient of B 2 ). In other words, the sum of the coefficients of these terms must be zero. The sum of the numerators of these two coefficients is proportional to: Thus the equivalency is valid when x eq~0 , i.e. when where a 0 3~a 3 =a 2 . This expression is valid for A=T~0, and similar expressions could be obtained for non-zero A=T but they become cumbersome. For a 3 =a 2~0 :2234 employed in Fig. 5, (W =T) eq~0 :082. A similar analysis carried out for the antisynchronous state would yield the same value. The equivalency is less accurate at large value of W =T as is also evident from Fig. 5(b, panel for W =T~0:2). We conclude that negative skewness is in general favorable for antisynchrony and positive skewness is not. Synchrony is stable predominantly when the left PRC jump is negative, or when the right PRC jump is positive. But for sufficiently positive skewness synchrony can be achieved for any B and B 2 , and for sufficiently positive skewness antisynchrony can be made unstable.

Network Simulations
We demonstrate here numerically that a PRC with left jump causes synchrony lose its stability using the network of identical modified Wang-Buzsáki model neurons coupled electrically with all-to-all connectivity such that the current I ele received by i th neuron is given by where N is the number of neurons in the network, and g ele is the strength of electrical coupling. The coupling can be termed weak if in the synchronized state the frequency of the neurons is not significantly altered. Each neuron in the network is made to oscillate by employing an appropriate steady external current such that it's spontaneous firing rate is about 5:33 Hz. We test here two parameter sets: w n~2 , I app~0 :17791 mA/cm 2 , and w n~9 , I app~0 :17 mA/cm 2 . By making the dynamics of the potassium activation faster, we made the nearly symmetric PRC [curve marked 1 in Fig. 6(a)] acquire a sharp rise at zero phase [curve marked 2 in Fig. 6(a)]. The PRC does acquire negative values at early phases below 0:23 ms, and then rises with a finite slope. But even if this regime was reset such that the PRC gradually approached a positive value at zero phase (corresponding to a left side jump), the results would not differ much. In addition to the sharp rise near zero phase, making the potassium dynamics faster also shifted the peak position of the PRC toward left. From the theory presented thus far, we expect that both these attributes help destabilize synchrony. The corresponding growth functions are shown in Fig. 6(b) that clearly reveal that left side jump together with the leftward tilt of the PRC caused a destabilization of the synchrony. The antisynchrony remains unstable. A non-zero phase-locked state has acquired stability. Actual solving for the time course using two coupled neurons (N = 2) [ Fig. 6(c)] for these two parameter sets verifies this observation.
A network of 100 neurons are coupled electrically starting from initial conditions distributed uniformly on the periodic orbit of the uncoupled neuron, and numerically integrated. Their spike times are shown in Fig. 6(e, top two panels) corresponding to the two parameter sets discussed above. The initial conditions constitute the most inhomogeneous set possible in the network, and the transients decay slowly but gradually and yield either synchronous state for the first parameter set, or an asynchronous state for the second parameter set.
The non-zero phase-locked state is a more sensitive function of the voltage time course and the PRC shape than the synchronous and antisynchronous states. This is because the exact level of the steady state separation between two coupled neurons corresponding to a non-zero phase-locked state is not identical for all parameter changes, unlike synchrony and antisynchrony. Thus their stability is also difficult to compute. But we noted earlier that the case of zero spike width leads to a stable phase-locked state merging with an unstable synchronous state. It is indeed difficult to carry out simulations in such a regime, but a stable phase-locked state can be created near (but not merging with) an unstable synchronous state. In Fig. 6(d) we have used a third parameter set to find that both synchrony and antisynchrony are unstable, but a non-zero phase-locked state very near synchronous state is stable. Carrying out a similar simulation as before using a 100 neuron network [ Fig. 6(e, bottom)] reveals that an almost-synchrony may be achieved where the spike times continue to show a small amount of jitter (proportional to the level of phase-locked state), but the network exhibits long transients. The dynamics of such a state depends on a number of factors including how fast a two neuron subnetwork repels a synchronous state (slope of the growth function at zero phase), and how fast such a network approaches the non-zero phase-locked state, and whether there are any other locked states in the system and their stability.

Special PRCs
Flat PRCs. A special case of Y (t) is when the PRC is a flat horizontal line: B~B 2~C . All the eigenvalues corresponding to synchronous and antisynchronous states become zero under this assumption. We can also see the emergence of zero eigenvalues directly from Eqs. 6 and 7. After effecting a transformation of variables in the definition of G(w), the integral equations expressing l and c can be written with derivatives of Y (t) which would become zero under the above assumption. Thus a flat PRC results in neutrally stable synchronous as well as antisynchronous states.
Linearly increasing or decreasing PRCs. A special PRC type that is linearly increasing or decreasing from left to right is obtained by setting the PRC maximum to the average of the two jumps: in Eq. 1. First consider stability of synchrony. The two regimes that were treated earlier (0ƒW v T 4 and into one because the earlier two PRC segments will now become just one due to the special value C. Some of the components r 1 , r 2 , r 3 , and r 4 (see Methods) will have different values than s 1 , s 2 , s 3 , and s 4 because of altered limits of integration, but their sum will become identical to one another. Substituting for C from Eq. 17 in the eigenvalue components r 1 , r 2 , r 3 , and r 4 , and computing the eigenvalue, we get l 1 as E(B{B 2 )(T{2W )(a 2 x 20 {a 3 )=T 2 where . Solving l 1~0 clearly yields the critical curve as: The coefficient of B is positive if a 3 =a 2 vx 20 which clearly is the case because it is sufficient for x 20 to be bigger than 1 for this relation to hold, but in fact x 20 §2 in the range 0ƒW =Tv2=5 (because x 20 acquires a value of 2 at W =T~0 and 5 at W =T~2=5, and has a positive slope of 3T=½2(T{2W ) 2 in between). Thus large negative B results in l 1 v0 and thus stable synchrony, and large positive B results in unstable synchrony, and the curve B~B 2 separates the two parameter regions. In summary, the stable region for synchrony is given by The region of stability did not change with finite spike width although the eigenvalue expression does acquire a dependence on W =T because eigenvalue in either case depends on the difference of the PRC jumps, and the spike width modulates that difference. Next consider stability of antisynchrony. Consider the case 0vW v T 4 . Substituting for C from Eq. 17 in the expression for c 1 and simplifying, we get c 1 as (B 2 {B)W ½(5a 2 {a 3 )(2T{5W )z 15a 3 W =½2T 2 (2T{5W ) which has clearly a negative coefficient for B. Thus at large negative B, the eigenvalue becomes positive leading to unstable antisynchrony, and at large positive B the eigenvalue becomes negative leading to stable antisynchrony. The critical curve that separates these two regimes is obtained by solving c 1~0 which yields B~B 2 . Next in the regime , substituting for C again from Eq. 17 in the expression for c 2 , we get the eigenvalue as (B 2 {B) 2a 3 (T{2W )W za 2 T(T{4W )z5a 2 W 2 Â Ã =(2T 2 W ) which again has a negative coefficient for B because within the bracketed expression each term is positive due to the current range of W . Thus the stability region in the same as above. In summary the region of stability of antisynchrony is given by We can see that at W~0, c 1 becomes zero independent of B or B 2 leading to neutrally stable antisynchrony at W~0. Note that the eigenvalue c 1 when expanded in Taylor series has a linear dependence on W at small W : Hence as W approaches zero, there is no sudden transition from stability to neutral stability. Thus for numerical purposes its stability at W~0 may be considered to follow the condition in Eq. 19. In summary, stable synchrony (Eq. 18) and stable antisynchrony (Eq. 19) exist in complementary parameter regions in (B,B 2 ) space, and the line B~B 2 represents critical state where both the locked states merge. It is also clear that the stability regions are independent of spike width and frequency.

Discussion and Conclusions
The PRC and the voltage are not completely independent. The PRCs are often computed numerically since analytical forms of the PRCs can only be obtained for very few simple models [3]. The standard approach is computing the so called adjoint using averaging technique. In the popular XPPAUT software program [32] this is computed numerically by using a method that uses backward integration technique. If the dynamical system under study exhibits sharp discontinuities such as in adaptive exponential integrate-and-fire model [25], a similar method was recently introduced by Ladenbauer et al. [21]. However, for these two techniques to be employed we must have full knowledge of the differential equations of all the variables in the system. In other words the PRC cannot be completely derived from the voltage time course alone, but is actually derived from the inverse of the derivative of the voltage, and is thus dependent on not only the voltage but all other variables and their derivatives computed on the trajectory.
To compute synchronization properties from experimentally measured PRCs, the above methods are not applicable because differential equation models are not known in general. But the theory of weakly coupled oscillators becomes advantageous if we note that the only two components that are responsible for determining the stability of synchrony and antisynchrony are the shape of the voltage and the shape of the PRC irrespective of all other variables. This gives us an opportunity to comprehensively study the role of each of the shape parameters of the voltage and the PRC on the network behavior without invoking any specific model system. We parameterized the shapes of both the voltage and the PRC. The PRC shape is parameterized by discontinuities and the skewness. And the voltage shape is parametrized by the spike peak, its minimum, threshold level, spike width and period. The Hodgkin-Huxley model voltage time course is only used for the limited purpose of extracting the relationship between the parameters. The information that is extracted is that the voltage may be divided into three segments, and a single spike width parameter W could be used in quantifying both the spike rise and fall profiles. Obviously, not all the PRC and voltage shapes spanned by our study are relevant for a particular given model. But we have presented analytical boundary expressions that can be used in applying the theory to many experimentally determined PRCs and voltages (within the purview of our model), and thus our theory has predictive power. Our study is not the first to incorporate variable parameters in theoretical models. Chow and Kopell [26] and Lewis and Rinzel [12] introduced free parameters of voltage spike width and shape to study the role of them within the context of mostly leaky integrate-and-fire (LIF) like models. But not many experimental results resemble PRCs of the LIF models. Thus we need to parametrize the PRC shapes as well to extend such theories to wider applicability. Our work is an attempt in this direction.
We have investigated the role of sharp jumps in the PRC at its edges on synchronization of electrically coupled neuronal oscillators. The PRC is modeled using a two-piecewise linear function with discontinuous jumps at the edges, i.e. phases corresponding to time 0 and the period T. The temporal relationships between the voltage segments are empirically derived from the Hodgkin-Huxley model. But the spike frequency, the spike period, the spike height parameters, and the level of PRC jumps are all freely altered. Analytical boundaries for stability of synchrony and antisynchrony are determined. The main results for small spike widths and frequencies are given by r 3 for synchrony and s 2 for antisynchrony (relations in Eq. 11 and 13). At large levels of W =T, the boundaries are given by r 4 for synchrony, and s 3 and s 4 for antisynchrony. We have also studied some special cases of the PRC where its level is either constant, linearly increasing or decreasing.
A positive jump at the left edge contributes to destabilizing the synchronous state, whereas a positive jump at the right edge contributes to its stability. When the depolarizing phase of V (t) has zero slope (i.e. when a 3~0 ), then only the spike upstroke and downstroke play role in deciding the stability, and the eigenvalue for stability is proportional to the difference of the jumps, and more specifically B{B 2 ; Thus if both edges have equal jumps, then the synchrony is neutrally stable (but the adjacent unstable region has stable non-zero phase-locked states near synchrony), and stability is achieved when the right jump is bigger than the left jump. This relationship holds even for negative jumps. If the depolarizing phase of the voltage has finite slope (a 3 w0), which usually is the case, then it contributes more favorably to synchrony. This leads to the curves r 1 that have decreasing slopes with increasing a 3 =a 2 as shown in Fig. 3(a). Large spike width/ frequency (W =Tw0) makes it harder to synchronize, and the stable region is moved to lower right as in Fig. 4(a,b).
The jumps at the edges do not directly contribute to altering the stability of antisynchronous state, but they affect the slope of the PRC segments that are connected to them, and those segments alter the stability. But the parameter that affects the stability of antisynchronous state in a sensitive manner is W =T. In the absence of spike width antisynchrony must be very rare because it occurs only when the sum of the left and right jumps exceeds twice the PRC level at t~T=2 (Eq. 10). It is unstable in the parameter region shown in Fig. 3(a) because the discontinuous jump associated with the voltage at t~T=2 imparts a positive eigenvalue component to the stability and it dominates the magnitude of the negative components from the depolarizing segment of the voltage. When the spike width is finite, the positive eigenvalue contribution from the downstroke (due to its less sharp slope) becomes smaller than the negative eigenvalue contribution from the upstroke (due to its sharper slope) leading to lessening of the destabilizing effect. At appropriate level of frequency (or W =T), the antisynchrony becomes stable. The overlap between stable synchrony and stable antisynchrony can become more at higher frequencies leading to bistability.
Coupled quadratic integrate-and-fire neuron models at different oscillation frequencies and different spike threshold levels were studied earlier by Pfeuty et al. [23,24]. Altering the ratio of the spike reset level and the threshold level, (for example from 0.1 to 10 in their Fig. 8 of [23]), transformed the PRC from a state of large jump on the left to a large jump on the right resulting in a transition from stable antisynchrony to unstable antisynchrony. Since we have special parameters to quantify the PRC jumps, in our notation, this means that increasing B=B 2 from above 1 toward zero should make the antisynchrony unstable. This is exactly the case in Fig. 4(a) (refer to the first quadrant). Our boundaries r 3 and s 2 (Eqs. 11 and 13) also quantify the effect of spike width and frequency more explicitly.
The classic leaky integrate-and-fire (LIF) neurons display a PRC (which is e t=T =I where T is the period, and I is the applied current) that is indeed discontinuous at either end, but also exponential [12,26]. Using the PWL methodology we see that BvB 2 for such a model which ensures stable synchrony (relation in Eq. 18), but relation in Eq. 19 predicts instability for antisynchrony. Our results on the role of skewness (Fig. 5) also indicated that for PRCs that have large positive skewness, antisynchrony is unstable. If the PRC had a diminished, or no response, at early phases, then large PRC skewness could reduce the effectiveness of the spike downstroke (when the spike times are shifted by half period) and cause the antisynchrony become stable.
The piecewise linear approach has been very successful in predicting the stability of synchronous and antisynchronous states. We also note that similar piecewise linear approach when synaptic excitation or inhibition is employed yields quite successful predictions of the role of PRC and voltage shapes on synchrony and antisynchrony [33]. But the non-zero phase-locked states and their stability may not be accurately predicted by a simple twopiecewise linear PRC model because of the fact that the electrical coupling invokes the entire voltage time course. A more elaborate PRC model is needed to address those states. Also, only the first order PRCs are included in the analysis. When higher order PRCs play a significant role in determining the discontinuous jumps, then the approach to and the stability of the phase-locked states may also be affected by their shapes. Finally, strong synaptic input or strong coupling can cause other dynamical solutions that are beyond the scope of the weakly coupled oscillators, such as n : m phase-locking that is different from 1 : 1, and dynamic spike order switches [34]. Near-synchronous or other non-zero phase-locked solutions could be affected by including the effects of higher-order PRCs [34]. However, weak coupling may still be accurate in predicting the stability of synchronous and antisynchronous states.

Eigenvalue for Synchronous State
Case (a). 0ƒW v T 4 . Consider the case of the spike downstroke occurring before the PRC peak, i.e. 0ƒW v T 4 . The integral in Eq. 6 can be split into four regimes each of which contains linear segments of Y (t) and V (t). These four regimes contribute to the following four eigenvalue components: a 3 (T{4W )½B(T{4W )zC(Tz4W ) T 2 (2T{5W ) , a 3 (T{W )½B 2 (T{W )zC(TzW ) T 2 (2T{5W ) , The combined eigenvalue is the sum of the eigenvalue components, and is obtained as.  W w0 leading to x 1 w0. And utilizing the condition in Eq. 3, we conclude that the coefficient of B in the eigenvalue expression is positive, and hence, at any fixed value of B and other parameters, at large positive B the eigenvalue becomes positive, and at large negative B, the eigenvalue becomes negative. Hence the synchronous state is unstable at large positive B, and stable at large negative B. Also since the dependence on B is linear, the derivative of l with respect to B is simply the coefficient of B which is positive. Thus the eigenvalue transitions from a negative to a positive value as B is increased across any critical curve B Ã . Such critical curve is obtained by solving the equation l 1~0 for B, and the stable synchronous region falls below this critical curve. That is the region of stable synchrony is given by BvB Ã and is expressed in Eq.11.
Case (b). . The limits of the four components computed in the previous case are altered because now the spike minimum occurs after the PRC peak position. The four eigenvalue components that contribute to the combined eigenvalue are given below.
Using the formula in Eq. 6 for these four components, we get the total eigenvalue (we term it l 2 ) as. Since the coefficient of B is positive, we directly see that at large negative B, the synchronous state becomes stable because the eigenvalue would become negative, and at large positive B it comes unstable. And the transition across any critical curve B Ã in between is governed by the coefficient of B which is positive. And hence the stable synchronous region lies below the critical curve that is obtained by solving the equation l 2~0 : BvB Ã and is as written in Eq. 12.

Eigenvalue for Antisynchronous State
Case (a). 0ƒW v T 4 . Finally we find conditions for stability of antisynchronous state. As before, first consider the case 0ƒW v T 4 . The eigenvalue integral in Eq. 7 can be split into four regimes that have PWL segments, and the four integrals give rise to the following four components: Using the formula in Eq. 7 for these four components, we obtain the total eigenvalue as. We clearly see (utilizing the condition in Eq. 4) that the coefficient of B in the eigenvalue expression is negative, and thus for large negative B, the eigenvalue becomes positive resulting in an unstable antisynchronous state, and for large positive B it becomes negative resulting in a stable antisynchronous state. This is exactly opposite to what we obtained above for the synchronous state. And the transition of the eigenvalue as B is increased is such that it goes from positive to negative (the derivative of c 1 with respect to B is negative). Thus stable antisynchronous state lies above a critical curve B Ã which is obtained by solving the equation c 1~0 for B: BwB Ã which is expressed in Eq. 13.
. Next consider the case of T 4 ƒW v 2T 5 . The four eigenvalue components are given by w 1~{ 2 T